[Toss-devel-svn] SF.net SVN: toss:[1721] trunk/Toss
Status: Beta
Brought to you by:
lukaszkaiser
From: <luk...@us...> - 2012-06-05 20:14:09
|
Revision: 1721 http://toss.svn.sourceforge.net/toss/?rev=1721&view=rev Author: lukaszkaiser Date: 2012-06-05 20:13:57 +0000 (Tue, 05 Jun 2012) Log Message: ----------- The Cash-Karp non-stiff adaptive RK code. Modified Paths: -------------- trunk/Toss/Formula/Formula.ml trunk/Toss/Formula/Formula.mli trunk/Toss/Formula/FormulaTest.ml trunk/Toss/examples/Parsing.toss Modified: trunk/Toss/Formula/Formula.ml =================================================================== --- trunk/Toss/Formula/Formula.ml 2012-06-05 10:39:06 UTC (rev 1720) +++ trunk/Toss/Formula/Formula.ml 2012-06-05 20:13:57 UTC (rev 1721) @@ -637,3 +637,107 @@ let k4 = eq_f (tinit +. tstep) (aadd vals_arr (amul tstep k3)) in let (dk2, dk3) = (amul (2.) k2, amul (2.) k3) in aadd vals_arr (amul tstepdiv6 (aadd k1 (aadd dk2 (aadd dk3 k4)))) + +let dist a1 a2 = + let d = Aux.array_fold_left2 (fun d x y -> d+.(x-.y)*.(x-.y)) 0. a1 a2 in + d ** 0.5 + +let sf = 0.9 + +let rkCK_default_start () = (ref 1.5, ref 1.1, ref 100., ref 100.) + +(* Implements the Cash-Karp algorithm with varying order and adaptive step, + precisely as described in the paper "A Variable Order Runge-Kutta Method + for Initial Value Problems with Varying Right-Hand Sides" by J. R. Cash and + Alan H. Karp, ACM Trans. Math. Software, 1990. + The interface is as in rk4_step above, except now the tolerance is used + (given by [epsilon]) and some initial references must be set (use the ones + from rkCK_default_start). We also return time passed and next step size.*) +let rkCK_step ?(epsilon=0.0001) (twiddle1, twiddle2, quit1, quit2) tn h f yn = + let k1 = f tn yn in + let k2 = f (tn +. (h/.5.)) (aadd yn (amul (h/.5.) k1)) in + let y1 = aadd yn (amul h k1) in + let y2 = aadd yn (amul h (aadd (amul (-3./.2.) k1) (amul (5./.2.) k2))) in + let e1 = ((dist y2 y1) /. epsilon) ** (1./.2.) in + if e1 > !twiddle1 *. !quit1 then (* abandon the step *) + let esttol = e1 /. !quit1 in + (yn, 0., h *. (max (1. /. 5.) (sf /. esttol))) + else + let k3 = f (tn +. (3.*.h/.10.)) + (aadd yn (aadd (amul (3.*.h/.40.) k1) (amul (9.*.h/.40.) k2))) in + let k4 = f (tn +. (3.*.h/.5.)) + (aadd yn (aadd (amul (3.*.h/.10.) k1) + (aadd (amul (-9.*.h/.10.) k2) (amul (6.*.h/.5.) k3)))) in + let y3 = aadd yn (amul h + (aadd (amul (19./.54.) k1) + (aadd (amul (-10./.27.) k3) (amul(55./.54.) k4)))) in + let e2 = ((dist y3 y2) /. epsilon) ** (1./.3.) in + if e2 > !twiddle2 *. !quit2 then (* try a lower order solution *) + if e1 < 1. then (* check the error of the second-order solution *) + let res15 = aadd yn (amul (h /. 10.) (aadd k1 k2)) in + let res15bar = aadd yn (amul (h /. 5.) k1) in + if dist res15 res15bar < epsilon then (* accept second-order solution *) + (res15, h /. 5., h /. 5.) + else (* abandon the step *) + (yn, 0., h /. 5.) + else (* abandon the step *) + let esttol = e2 /. !quit2 in + (yn, 0., h *. (max (1. /. 5.) (sf /. esttol))) + else + let k5 = f (tn +. h) + (aadd yn (aadd (amul (-11.*.h/.54.) k1) + (aadd (amul (5.*.h/.2.) k2) + (aadd (amul (-70.*.h/.27.) k3) + (amul (35.*.h/.27.) k4))))) in + let k6 = f (tn +. (7.*.h/.8.)) + (aadd yn (aadd (amul (1631.*.h/.55296.) k1) + (aadd (amul (175.*.h/.512.) k2) + (aadd (amul (575.*.h/.13824.) k3) + (aadd (amul (44275.*.h/.110592.) k4) + (amul (253.*.h/.4096.) k5)))))) in + let y4 = aadd yn (amul h + (aadd (amul (2825./.27648.) k1) + (aadd (amul (18575./.48384.) k3) + (aadd (amul (13525./.55296.) k4) + (aadd (amul (277./.14336.) k5) + (amul (1./.4.) k6)))))) in + let y5 = aadd yn (amul h + (aadd (amul (37./.378.) k1) + (aadd (amul (250./.621.) k3) + (aadd (amul (125./.594.) k4) + (amul (512./.1771.) k6))))) in + let e4 = ((dist y5 y4) /. epsilon) ** (1./.5.) in + if e4 > 1. then ( (* try a lower order solution *) + (* readjust twiddle factors *) + if e1 /. !quit1 < !twiddle1 then twiddle1 := max 1.1 (e1 /. !quit1); + if e2 /. !quit2 < !twiddle2 then twiddle2 := max 1.1 (e2 /. !quit2); + if e2 < 1. then ( (* try the order-3 solution *) + let res35 = aadd yn (amul h + (aadd (amul (1./.10.) k1) + (aadd (amul (2./.5.) k3) + (amul (1./.10.) k4)))) in + let res35bar = aadd yn (amul (3. *. h /. 5.) k3) in + if dist res35 res35bar < epsilon then (* accept order-3 solution *) + (res35, 3. *. h /. 5., 3. *. h /. 5.) + else (* try an even lower order solution *) + if e1 < 1. then (* check the error of the order-2 solution *) + let res15 = aadd yn (amul (h /. 10.) (aadd k1 k2)) in + let res15bar = aadd yn (amul (h /. 5.) k1) in + if dist res15 res15bar < epsilon then (* accept order-2 solution*) + (res15, h /. 5., h /. 5.) + else (* abandon the step *) + (yn, 0., h /. 5.) + else (* abandon the current step *) + (yn, 0., h *. (max (1. /. 5.) (sf /. e4))) + ) else (* abandon the current step *) + (yn, 0., h *. (max (1. /. 5.) (sf /. e4))) + ) else ( (* accept the full order-5 solution *) + let q1, q2 = e1 /. e4, e2 /. e4 in + let q1 = if q1 > !quit1 then min q1 (10. *. !quit1) else + max q1 (2. *. !quit1 /. 3.) in + let q2 = if q2 > !quit2 then min q2 (10. *. !quit2) else + max q2 (2. *. !quit2 /. 3.) in + quit1 := max 1. (min 10000. q1); + quit2 := max 1. (min 10000. q2); + (y5, h, h *. (min 5. (sf /. e4))) + ) Modified: trunk/Toss/Formula/Formula.mli =================================================================== --- trunk/Toss/Formula/Formula.mli 2012-06-05 10:39:06 UTC (rev 1720) +++ trunk/Toss/Formula/Formula.mli 2012-06-05 20:13:57 UTC (rev 1721) @@ -145,3 +145,15 @@ side [eq_terms]. Time variable [tvar] starts at [tinit] and moves [tstep]. *) val rk4_step : float -> float -> (float -> float array -> float array) -> float array -> float array + +(** Implements the Cash-Karp algorithm with varying order and adaptive step, + precisely as described in the paper "A Variable Order Runge-Kutta Method + for Initial Value Problems with Varying Right-Hand Sides" by J. R. Cash and + Alan H. Karp, ACM Trans. Math. Software, 1990. + The interface is as in rk4_step above, except now the tolerance is used + (given by [epsilon]) and some initial references must be set (use the ones + from rkCK_default_start). We also return time passed and next step size.*) +val rkCK_default_start : unit -> float ref * float ref * float ref * float ref +val rkCK_step: ?epsilon: float -> float ref* float ref* float ref* float ref -> + float -> float -> (float -> float array -> float array) -> + float array -> float array * float * float Modified: trunk/Toss/Formula/FormulaTest.ml =================================================================== --- trunk/Toss/Formula/FormulaTest.ml 2012-06-05 10:39:06 UTC (rev 1720) +++ trunk/Toss/Formula/FormulaTest.ml 2012-06-05 20:13:57 UTC (rev 1721) @@ -81,5 +81,35 @@ "0.00110, 3.6e-10, 2.99999, -1.9991, 7.16443, -0.0008" (float_arr_str 7 (rk4_step 0. 0.02 ceqs init)); ); + + "rkCK" >:: + (fun () -> + let float_arr_str n fa = String.concat ", " (List.map (fun f -> + String.sub (string_of_float f) 0 n) (Array.to_list fa)) in + let eqs = eqs_of_string ":f(a)' = :f(a) + :t" in + let ceqs = compile ":t" eqs in + let (res, _, _) = rkCK_step (rkCK_default_start()) 0. 0.1 ceqs [|0.|] in + assert_equal ~printer:(fun x -> x) "0.00517" (float_arr_str 7 res); + + let tyson_model_eqs = eqs_of_string + (" :y(e1)' = 0.015 + -200. * :y(e1) * :y(e4); " ^ + ":y(e2)' = -0.6 * :y(e2) + 1. * :y(e5); " ^ + ":y(e3)' = -100. * :y(e3) + 1. * :y(e5) + 100. * :y(e4); " ^ + ":y(e4)' = -100. * :y(e4) + 100.*:y(e3) + -200.*:y(e1) * :y(e4); "^ + ":y(e5)' = -1. * :y(e5) + 0.018 * :y(e6) + " ^ + " 180. * :y(e6) * :y(e5) * :y(e5); " ^ + ":y(e6)' = -0.018 * :y(e6) + 200. * :y(e1) * :y(e4) + " ^ + " -180. * :y(e6) * :y(e5) * :y(e5)") in + let ceqs = compile ":t" tyson_model_eqs in + let init = [| 0.; 0.; 1.; 0.; 0.; 0. |] in + (* Step in this try is too big - should abandon, i.e. time passed = 0.*) + let (_, t, _) = rkCK_step (rkCK_default_start()) 0. 0.02 ceqs init in + assert_equal ~printer:(fun x -> string_of_float x) 0. t; + (* Now it is ok, should report results. *) + let (res, _, _) = rkCK_step (rkCK_default_start()) 0. 0.003 ceqs init in + assert_equal ~printer:(fun x -> x) + "4.28845, 1.90309, 0.77440, 0.22559, 3.00456, 2.11544" + (float_arr_str 7 res); + ); ] Modified: trunk/Toss/examples/Parsing.toss =================================================================== --- trunk/Toss/examples/Parsing.toss 2012-06-05 10:39:06 UTC (rev 1720) +++ trunk/Toss/examples/Parsing.toss 2012-06-05 20:13:57 UTC (rev 1721) @@ -13,4 +13,4 @@ } START [ | Tp:2 {}; Tlist:1 {}; Pone (one); Ptrue (t); Pnil (nil); - S { (one, t); (t, nil) } | ] + S { (one, t); (t, nil) } | - ] This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |