Here is a Runge Kutta integrator (using the Cash-Karp 5th order
method) that I know works. But it is very un-J-like: I just
transcribed the Fortran from Numerical Recipes, as I was in a
hurry to get something working :-)
Patrick
On Tue, 6 May 2014, David Lambert 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
