#] #] ********************* #] loaddefs link d_Qndfs 'Cash-Karp-Runge-Kutta integrate ODEs with accuracy monitoring.ndf - #] adapted rkck Cash-Karp-Runge-Kutta step used by rkqs, #] rkqs integrate one step of ODEs with accuracy monitoring # 18May07 www.BillHowell.ca # Warning: Use Courier constant-width font and no-line-wrap to view this file! # Otherwise the code doesn't line up and it's hard to read. f_MindSym := 'Cash-Karp-Runge-Kutta integrate ODEs with accuracy monitoring.ndf' ; loaddefs_start f_MindSym ; #**************************** # List of operators, generated with : # $ cat "$d_Qndfs""Cash-Karp-Runge-Kutta integrate ODEs with accuracy monitoring.ndf" | grep "^#]" | sed 's/^#\]/ /' # #*********************************** # Original code from: ?Editor? "Numerical Recipes for C: The Art of Scientific Computing" Cambridge University Press, http://www.nrbook.com/a/ #********************************* # globals & identifier sources for tracing program elements ; f_Cash_Karp_Runge_Kutta := 'Cash-Karp-Runge-Kutta integrate ODEs with accuracy monitoring.ndf' ; # variables defined in this file a1 := 0.0 ; c1 := 37.0/378.0 ; dc1 := c1 - (2825.0/27648.0) ; a2 := 0.2 ; c2 := 0.0 ; dc2 := 0.0 ; a3 := 0.3 ; c3 := 250.0/621.0 ; dc3 := c3 - (18575.0/48384.0) ; a4 := 0.6 ; c4 := 125.0/594.0 ; dc4 := c4 - (13525.0/55296.0) ; a5 := 1.0 ; c5 := 0.0 ; dc5 := -277.00/14336.0 ; a6 := 0.875 ; c6 := 512.0/1771.0 ; dc6 := c6 - 0.25 ; b21 := 0.2 ; b31 := 3.0/40.0 ; b32 := 9.0/40.0 ; b41 := 0.3 ; b42 := -0.9 ; b43 := 1.2 ; b51 := -11.0/54.0 ; b52 := 2.5 ; b53 := -70.0/27.0 ; b54 := 35.0/27.0 ; b61 := 1631.0/55296.0; b62 := 175.0/512.0 ; b63 := 575.0/13824.0 ; b64 := 44275.0/110592.0 ; b65 := 253.0/4096.0 ; a1_v := c1_v := dc1_v := f_Cash_Karp_Runge_Kutta ; a2_v := c2_v := dc2_v := f_Cash_Karp_Runge_Kutta ; a3_v := c3_v := dc3_v := f_Cash_Karp_Runge_Kutta ; a4_v := c4_v := dc4_v := f_Cash_Karp_Runge_Kutta ; a5_v := c5_v := dc5_v := f_Cash_Karp_Runge_Kutta ; a6_v := c6_v := dc6_v := f_Cash_Karp_Runge_Kutta ; b21_v := f_Cash_Karp_Runge_Kutta ; b31_v := b32_v := f_Cash_Karp_Runge_Kutta ; b41_v := b42_v := b43_v := f_Cash_Karp_Runge_Kutta ; b51_v := b52_v := b53_v := b54_v := f_Cash_Karp_Runge_Kutta ; b61_v := b62_v := b63_v := b64_v := b65_v := f_Cash_Karp_Runge_Kutta ; ERRCON := 1.89e-4 ; ERRCON_v := f_Cash_Karp_Runge_Kutta ; MAXSTP := 10000 ; MAXSTP_v := f_Cash_Karp_Runge_Kutta ; PGROW := -0.2 ; PGROW_v := f_Cash_Karp_Runge_Kutta ; PSHRNK := -0.25 ; PSHRNK_v := f_Cash_Karp_Runge_Kutta ; SAFETY := 0.9 ; SAFETY_v := f_Cash_Karp_Runge_Kutta ; TINY := 1.0e-30 ; TINY_v := f_Cash_Karp_Runge_Kutta ; % kmax don't need as is argument to an operation? ; # operators defined in this file sign_o := f_Cash_Karp_Runge_Kutta ; rkck_o := f_Cash_Karp_Runge_Kutta ; rkqs_o := f_Cash_Karp_Runge_Kutta ; odeint_o := f_Cash_Karp_Runge_Kutta ; # externals IF flag_debug THEN write 'loading derivs' ; ENDIF ; #] derivs IS EXTERNAL operation ; derivs_e := f_Cash_Karp_Runge_Kutta ; - derivs IS EXTERNAL operation ; derivs_e := f_Cash_Karp_Runge_Kutta ; # NONLOCAL grab bag, leftovers, changes NONLOCAL ERRCON MAXSTP PGROW PSHRNK SAFETY TINY ; NONLOCAL a1_v c1_v dc1_v a2_v c2_v dc2_v a3_v c3_v dc3_v a4_v c4_v dc4_v a5_v c5_v dc5_v a6_v c6_v dc6_v b21_v b31_v b32_v b41_v b42_v b43_v b51_v b52_v b53_v b54_v b61_v b62_v b63_v b64_v b65_v ; #******************************** General routines # sign - used for odeint IF flag_debug THEN write 'loading sign' ; ENDIF ; #] sign IS OP a b { IF b >= 0.0 THEN abs a ELSE opposite abs a ENDIF } - sign IS OP a b { IF b >= 0.0 THEN abs a ELSE opposite abs a ENDIF } #******************************** rkck - Cash-Karp-Runge-Kutta step used by rkqs calling: Outputs := Operator Inputs yout yerr := rkck x y dx dydx deriv_name ; IF flag_debug THEN write 'loading rkck' ; ENDIF ; #] rkck IS OP x y dx dydx deriv_name - rkck IS OP x y dx dydx deriv_name { NONLOCAL a1_v c1_v dc1_v a2_v c2_v dc2_v a3_v c3_v dc3_v a4_v c4_v dc4_v a5_v c5_v dc5_v a6_v c6_v dc6_v b21_v b31_v b32_v b41_v b42_v b43_v b51_v b52_v b53_v b54_v b61_v b62_v b63_v b64_v b65_v ; % ; ytemp := y + (dx * b21 * dydx ) ; ak2 := APPLY deriv_name ( (a2 * dx + x) ytemp ) ; ytemp := y + (dx * ((b31 * dydx) + (b32 * ak2) )) ; ak3 := APPLY deriv_name ( (a3 * dx + x) ytemp ) ; ytemp := y + (dx * ((b41 * dydx) + (b42 * ak2) + (b43 * ak3) )) ; ak4 := APPLY deriv_name ( (a4 * dx + x) ytemp ) ; ytemp := y + (dx * ((b51 * dydx) + (b52 * ak2) + (b53 * ak3) + (b54 * ak4) )) ; ak5 := APPLY deriv_name ( (a5 * dx + x) ytemp ) ; ytemp := y + (dx * ((b61 * dydx) + (b62 * ak2) + (b63 * ak3) + (b64 * ak4) + (b65 * ak5) )) ; ak6 := APPLY deriv_name ( (a6 * dx + x) ytemp ) ; yout := y + (dx * ((c1 * dydx) + (c3 * ak3) + (c4 * ak4) + (c6 * ak6) )) ; yerr := dx * ((dc1 * dydx) + (dc3 * ak3) + (dc4 * ak4) + (dc5 * ak5) + (dc6 * ak6) ) ; % ; yout yerr } #******************************** rkqs - integrate one step of ODEs with accuracy monitoring calling: Outputs := Operator Inputs x y yerr dx_next stop_msg := rkqs x y dx dydx eps yscal deriv_name IF flag_debug THEN write 'loading rkqs' ; ENDIF ; #] rkqs IS OP x y dx dydx eps yscal deriv_name - rkqs IS OP x y dx dydx eps yscal deriv_name { NONLOCAL ERRCON MAXSTP PGROW PSHRNK SAFETY TINY ; stop_flag := o ; % ; REPEAT y yerr := rkck x y dx dydx deriv_name ; errmax := (max abs (yerr / yscal)) / eps ; IF errmax <= 1.0 THEN stop_flag := l ; stop_msg := 'Next step, errmax <= 1.0, in routine rkqs' ; ELSE dx_temp := SAFETY * dx * (power errmax PSHRNK) ; IF dx >= 0.0 THEN dx := max dx_temp (0.1 * dx) ELSE dx := min dx_temp (0.1 * dx) ENDIF ; xnew := x + dx ; IF xnew = x THEN stop_flag := l ; stop_msg := fault '?stepsize underflow in routine rkqs' ; ENDIF ; ENDIF ; UNTIL stop_flag ENDREPEAT ; % ; IF errmax > ERRCON THEN dx_next := SAFETY * dx * (power errmax PGROW) ; % original *dx_next ; ELSE dx_next := 5.0 * dx ; ENDIF ; % ; (x + dx) y yerr dx dx_next stop_msg } #****************************************************** # odeint.ndf - Cash-Karp-Runge-Kutta integrate ODEs with accuracy monitoring, adapted ; calling: Outputs := Operator Inputs stop_msg kount nok nbad xp yp yerrp := odeint x1 x2 ystart dx1 dx_min dxsav eps kmax deriv_name IF flag_debug THEN write 'loading odeint' ; ENDIF ; #] odeint IS OP x1 x2 ystart dx1 dx_min dxsav eps kmax deriv_name - odeint IS OP x1 x2 ystart dx1 dx_min dxsav eps kmax deriv_name { NONLOCAL ERRCON MAXSTP PGROW PSHRNK SAFETY TINY ; % ; kount := nok := nbad := ntsp := yerr := 0 ; yp := yerrp := kmax (first shape ystart) reshape 0 ; xp := kmax reshape 0.0 ; xp@0 := x := x1 ; yp|[0,]:= y := ystart ; dx := sign dx1 (x2 - x1) ; xsav := x - (dxsav * 2.0) ; stop_flag := o ; % ; REPEAT dydx := APPLY deriv_name (x y) ; yscal := (abs y) + (abs (dydx * dx)) + TINY ; IF AND (kount < kmax-1) (abs (x - xsav) > (abs dxsav)) THEN kount := kount + 1 ; xp@kount := xsav := x ; yp|[kount,] := y ; yerrp|[kount,] := yerr ; ENDIF ; IF ((x + dx - x2) * (x + dx - x1)) > 0.0 THEN dx := x2 - x ENDIF ; x y yerr dx_did dx_next rkqs_msg := rkqs x y dx dydx eps yscal deriv_name ; IF isfault rkqs_msg THEN stop_flag := l ; stop_msg := rkqs_msg ; ELSE % inserted for tracing progression ; kount := kount + 1 ; xp@kount := xsav := x ; yp|[kount,] := y ; yerrp|[kount,] := yerr ; % end of tracing progression ; IF dx_did = dx THEN nok := nok + 1 ; ELSE nbad := nbad + 1 ; ENDIF ; ntsp := ntsp + 1 ; IF ntsp = MAXSTP THEN stop_flag := l ; stop_msg := fault '?Too many steps, nstp > MAXSTP, in routine odeint' ; ENDIF ; IF (abs dx) <= dx_min THEN stop_flag := l ; stop_msg := fault '?Stepsize too small, dx < dx_min, in routine odeint' ; ENDIF ; IF ((x - x2) * (x2 - x1)) >= 0.0 THEN stop_flag := l ; stop_msg := 'Answer was found in routine odeint' ; ENDIF ; dx := dx_next ; ENDIF ; UNTIL stop_flag ENDREPEAT ; % ; i := tell (kount + 1) ; stop_msg kount nok nbad xp|[i] yp|[i,] yerrp|[i,] } #*********************************************** Description of variables # arguments for rkqs: y dydx n x dx_try eps yscal dx_did dx_next deriv_name # arguments to odeint x1 starting point x2 end of interval eps used by rkqs to indicate when error criteria is met? dx1 used for initial stepsize, either positive or negative value taken depending on DIRECTION of start & end points of interval dx_min minimum step size allowable deriv_name external function to calculate derivatives? rkqs not used as an argument in my version, I simply implement a direct call loaddefs_ended f_MindSym ; # enddoc