[Toss-devel-svn] SF.net SVN: toss:[1721] trunk/Toss

[<luk...@us...>]

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.