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
