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

Last change on this file since 30954 was 30954, checked in by Ivan Raikov, 6 years ago

flsim: update to SML RK library

File size: 4.9 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
96fun gen_soln3 (integrator,h,t,st) =
97  let
98      val (stn,en,inp) = integrator (t,st)
99      val tn       = Real.+(t,h)
100  in
101      if t >= 5.0
102      then (putStr (showst (tn,stn));
103            putStrLn ("\t" ^ (showReal (inp 1.0))))
104      else gen_soln3 (integrator,h,tn,stn)
105  end
106
107fun do_case3 integrator n =
108  let
109      val h = if n < 0 then Real.Math.pow (2.0,Real.fromInt (~n))
110              else Real.Math.pow (0.5,Real.fromInt n)
111      val sep = if n <= 4 then "\t\t" else "\t"
112  in
113      putStr (String.concat [(showReal h), sep]);
114      gen_soln3 (integrator h,h,t0,y0)
115  end
116
117fun solver3 (integrator,stats) =
118  (putStrLn stats;
119   putStrLn "# step yf uf err";
120   List.app (do_case3 (integrator (scaler,summer,deriv)))
121            (List.tabulate (15, fn x => x - 2));
122   putStrLn "# All done!\n")
123
124
125val rkfe: real stepper1 = make_rkfe()
126val rk3:  real stepper1 = make_rk3()
127val rk4a: real stepper1 = make_rk4a()
128val rk4b: real stepper1 = make_rk4b()
129
130val rkhe:  real stepper2  = make_rkhe()
131val rkbs:  real stepper2  = make_rkbs()
132val rkf45: real stepper2 = make_rkf45()
133val rkck:  real stepper2 = make_rkck()
134val rkdp:  real stepper2 = make_rkdp()
135val rkdpb:  real stepper2 = make_rkdpb()
136val rkf78: real stepper2 = make_rkf78()
137val rkv65: real stepper2 = make_rkv65()
138
139val rkhe_aux:  real stepper1  = make_rkhe_aux()
140val rkbs_aux:  real stepper1  = make_rkbs_aux()
141val rkf45_aux: real stepper1 = make_rkf45_aux()
142val rkck_aux:  real stepper1 = make_rkck_aux()
143val rkdp_aux:  real stepper1 = make_rkdp_aux()
144val rkdpb_aux:  real stepper1 = make_rkdpb_aux()
145val rkf78_aux: real stepper1 = make_rkf78_aux()
146val rkv65_aux: real stepper1 = make_rkv65_aux()
147
148val cerkdp:  real stepper3 = make_cerkdp()
149
150
151fun run() =
153  List.app solver1 [(rkfe, show_rkfe),
154                    (rk3,  show_rk3),
155                    (rk4a, show_rk4a),
156                    (rk4b, show_rk4b)];
158  List.app solver2 [(rkhe, show_rkhe),
159                    (rkbs, show_rkbs),
160                    (rkf45, show_rkf45),
161                    (rkck, show_rkck),
162                    (rkdp, show_rkdp),
163                    (rkdpb, show_rkdpb),
164                    (rkf78, show_rkf78),
165                    (rkv65, show_rkv65)];
166  putStrLn "#### Continuous Solvers";
167  List.app solver3 [(cerkdp, show_cerkdp)];
168  putStrLn "#### Auxiliary Solvers: Error Estimators from Adaptives";
169  List.app solver1 [(rkhe_aux, show_rkhe_aux),
170                    (rkbs_aux, show_rkbs_aux),
171                    (rkf45_aux, show_rkf45_aux),
172                    (rkck_aux, show_rkck_aux),
173                    (rkdp_aux, show_rkdp_aux),
174                    (rkdpb_aux, show_rkdpb_aux),
175                    (rkf78_aux, show_rkf78_aux),
176                    (rkv65_aux, show_rkv65_aux)]
177  )
178
179val _ = run()
180end
Note: See TracBrowser for help on using the repository browser.