# source:project/release/4/flsim/trunk/sml-lib/rk/rktest.sml@22665

Last change on this file since 22665 was 22665, checked in by Ivan Raikov, 10 years ago

flsim: changes to Runge-Kutta solver interfaces

File size: 3.6 KB
Line
1(* Module for testing Runge-Kutta integration of ODEs
2   Based on code by  Uwe Hollerbach <uh@alumni.caltech.edu>
3*)
4
5structure RungeKuttaTest =
6struct
7
8open RungeKutta
9
10val summer = Real.+
11val scaler = Real.*
12
13infix 7 */
14infix 6 +/
15infix //
16
17
18
19(* Solve the test problem dy/dt = -t^3 y^3
20
21   which has the exact solution y = 1/sqrt(C + t^4/2)
22   where C = 1/y_0^2 - t_0^4/2
23*)
24
25val con = ~0.4
26fun deriv (t,y) = con*y
27val t0 = 0.0
28val y0 = 1.75
29fun exact t = y0*Real.Math.exp(con*(t - t0))
30
31fun putStr str =
32    (TextIO.output (TextIO.stdOut, str))
33
34fun putStrLn str =
35    (TextIO.output (TextIO.stdOut, str);
36     TextIO.output (TextIO.stdOut, "\n"))
37
38fun showReal n = Real.toString n
39
40fun showst (t, y) = String.concat [(showReal y), "\t", (showReal (y - (exact t)))]
41
42fun gen_soln1 (integrator,h,t,st) =
43  let
44      val stn = integrator (t,st)
45      val tn  = Real.+(t,h)
46  in
47      if t >= 5.0
48      then putStrLn (showst (tn,stn))
49      else gen_soln1 (integrator,h,tn,stn)
50  end
51
52fun gen_soln2 (integrator,h,t,st) =
53  let
54      val (stn,en) = integrator (t,st)
55      val tn       = Real.+(t,h)
56  in
57      if t >= 5.0
58      then putStrLn (showst (tn,stn))
59      else gen_soln2 (integrator,h,tn,stn)
60  end
61
62fun do_case1 integrator n =
63  let
64      val h = if n < 0 then Real.Math.pow (2.0,Real.fromInt (~n))
65              else Real.Math.pow (0.5,Real.fromInt n)
66      val sep = if n <= 4 then "\t\t" else "\t"
67  in
68      putStr (String.concat [(showReal h), sep]);
69      gen_soln1 (integrator h,h,t0,y0)
70  end
71
72fun solver1 (integrator,stats) =
73  (putStrLn stats;
74   putStrLn "# step yf err";
75   List.app (do_case1 (integrator (scaler,summer,deriv)))
76            (List.tabulate (15, fn x => x - 2));
77   putStrLn "# All done!\n")
78
79fun do_case2 integrator n =
80  let
81      val h = if n < 0 then Real.Math.pow (2.0,Real.fromInt (~n))
82              else Real.Math.pow (0.5,Real.fromInt n)
83      val sep = if n <= 4 then "\t\t" else "\t"
84  in
85      putStr (String.concat [(showReal h), sep]);
86      gen_soln2 (integrator h,h,t0,y0)
87  end
88
89fun solver2 (integrator,stats) =
90  (putStrLn stats;
91   putStrLn "# step yf err";
92   List.app (do_case2 (integrator (scaler,summer,deriv)))
93            (List.tabulate (15, fn x => x - 2));
94   putStrLn "# All done!\n")
95
96
97val rkfe: real stepper1 = make_rkfe()
98val rk3:  real stepper1 = make_rk3()
99val rk4a: real stepper1 = make_rk4a()
100val rk4b: real stepper1 = make_rk4b()
101
102val rkhe:  real stepper2  = make_rkhe()
103val rkbs:  real stepper2  = make_rkbs()
104val rkf45: real stepper2 = make_rkf45()
105val rkck:  real stepper2 = make_rkck()
106val rkdp:  real stepper2 = make_rkdp()
107val rkf78: real stepper2 = make_rkf78()
108
109val rkhe_aux:  real stepper1  = make_rkhe_aux()
110val rkbs_aux:  real stepper1  = make_rkbs_aux()
111val rkf45_aux: real stepper1 = make_rkf45_aux()
112val rkck_aux:  real stepper1 = make_rkck_aux()
113val rkdp_aux:  real stepper1 = make_rkdp_aux()
114val rkf78_aux: real stepper1 = make_rkf78_aux()
115
116
117
118fun run() =
120  List.app solver1 [(rkfe, show_rkfe),
121                    (rk3,  show_rk3),
122                    (rk4a, show_rk4a),
123                    (rk4b, show_rk4b)];
125  List.app solver2 [(rkhe, show_rkhe),
126                    (rkbs, show_rkbs),
127                    (rkf45, show_rkf45),
128                    (rkck, show_rkck),
129                    (rkdp, show_rkdp),
130                    (rkf78, show_rkf78)];
131  putStrLn "#### Auxiliary Solvers: Error Estimators from Adaptives";
132  List.app solver1 [(rkhe_aux, show_rkhe_aux),
133                    (rkbs_aux, show_rkbs_aux),
134                    (rkf45_aux, show_rkf45_aux),
135                    (rkck_aux, show_rkck_aux),
136                    (rkdp_aux, show_rkdp_aux),
137                    (rkf78_aux, show_rkf78_aux)]
138  )
139
140val _ = run()
141end
Note: See TracBrowser for help on using the repository browser.