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
