Yes, those too. Thank you.
On Wed, May 7, 2014 at 6:00 PM, Raul Miller <[email protected]> wrote: > Also http://www.eelvex.net/programming/the-choice-2-j-and-lisp/ for > your J specific comments and > http://www.eelvex.net/programming/the-choice-1-intro/ for your > introduction. > > Thanks, > > -- > Raul > > > On Wed, May 7, 2014 at 10:33 AM, EelVex <[email protected]> wrote: > > You can see here: > https://github.com/EelVex/thechoice/blob/master/sol.ijs > > for a simple Runge-Kutta-4 solver of an N-dimensional system of ODEs. > > > > I think converting to RKf45 is straightforward but haven't tried it so I > > may be wrong. > > > > > > On Wed, May 7, 2014 at 12:30 AM, David Lambert <[email protected] > >wrote: > > > >> (Unless I've overlooked it) I propose we include a new adverb into > >> addons/math/misc to solve systems of first order initial value > differential > >> equations. This initial attempt I believe meets the "stage 1 > requirements" > >> of product development, summarized simply as "make one that works". In > >> > >> u Rkf45 > >> > >> u is a dyad. x is the grid position at which to evaluate the function > >> (time in the pendulum example), y are the estimated function values at > that > >> grid (angle and angular speed in the pendulum example). u computes and > >> returns the derivatives (angular speed and angular acceleration in the > >> pendulum example). Since the example hasn't got a forcing function, the > >> pendulum is not explicitly dependent on time. > >> > >> u Rkf45 is also dyadic. x is the rank 1 grid (times at which to advance > >> the solution in the pendulum example), y are the initial values (angle > and > >> angular speed in the pendulum example). The grid defaults to 20 steps > from > >> 0 to 1. > >> > >> u Rkf45 returns a vector of 2 boxes. The head is a matrix GRID ,. > >> SOLUTION . (In the pendulum example each row is time, angle, angular > >> speed.) The tail is a matrix of estimated 5th order error. In the > example > >> each row is angular_error, angular_speed_error. > >> > >> The complete grid left input to Rkf45 is unusual. A more standard input > >> would be start, stop, step_size. Better still would be the ends of the > >> interval, an error tolerance, maximum number of function evaluations > with > >> an adaptive algorithm. Let's make this a great code and enter it into > >> math/misc, or tell me where it's already been done. The output may > also be > >> non-standard and we can fix that too. > >> > >> Thanks, Dave. > >> > >> > >> NB. Rkf45.ijs > >> > >> Rkf45 =: 1 : 0 > >> NB. T u Rkf45 U0 > >> NB. dyad u computes derivatives > >> NB. x is the grid > >> NB. y are the initail function values. > >> y u Rkf45~ (%~ i.@:>:)20 > >> : > >> float =. _1&x: > >> mp =. +/ .* > >> A =. ,:,0 > >> A =. A,1r4 0 > >> A =. A,3r32 9r32 > >> A =. A,1932r2197 _7200r2197 7296r2197 > >> A =. A,439r216 _8 3680r513 _845r4104 > >> A =. A,_8r27 2 _3544r2565 1859r4104 _11r40 > >> C =. float +/"1 A > >> B4=. float 25r216 0 1408r2565 2197r4104 _1r5 0 > >> B5=. float 16r135 0 6656r12825 28561r56430 _9r50 2r55 > >> A =. float A > >> GRID =. x > >> ERROR =. RESULT =. ,: y > >> a =. A&({~ [: < (; i.)) > >> for_I. i. <: # GRID do. > >> H =. -/ GRID {~ 1 0 + I > >> T =. (I { GRID) + H * C > >> Y =. {: RESULT > >> K =. ,: H * (0 { T) u Y > >> K =. K , (H * {&T u Y + K mp~ a) 1 > >> K =. K , (H * {&T u Y + K mp~ a) 2 > >> K =. K , (H * {&T u Y + K mp~ a) 3 > >> K =. K , (H * {&T u Y + K mp~ a) 4 > >> K =. K , (H * {&T u Y + K mp~ a) 5 > >> RESULT =. RESULT , Y + B4 mp K > >> ERROR =. ERROR , Y + B5 mp K > >> end. > >> (GRID ,. RESULT) ; ERROR - RESULT > >> ) > >> > >> > >> NB. test case > >> > >> TAU =: 2p1 NB. tauday.com > >> > >> NB. small amplitude period, ACCELERATION %:@:% LENGTH > >> PERIOD =: 1 > >> ANGULAR_FREQUENCY =: TAU % PERIOD > >> > >> pendulum =: 3 : 0 :([: $: ]) > >> NB. A second order differential equation > >> NB. converted to two first order equations. > >> 'PHI U' =: y > >> DPHI_DT =: U > >> DU_DT =: - (*: ANGULAR_FREQUENCY) * 1 o. PHI > >> DPHI_DT , DU_DT > >> ) > >> > >> NB. 3.14 0 means "start near the top, stationary" > >> BIG_SWING =: (100%~i.1001) pendulum Rkf45 3.14 0 > >> > >> NB. a small amplitude swing > >> ONE_PERIOD =: pendulum Rkf45 0.01 0 > >> > >> require 'plot' > >> > >> plot;/2{.|:>{.BIG_SWING > >> plot;/2{.|:>{.ONE_PERIOD > >> > >> ---------------------------------------------------------------------- > >> For information about J forums see http://www.jsoftware.com/forums.htm > >> > > ---------------------------------------------------------------------- > > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
