1 | (* Module for testing Runge-Kutta integration of ODEs |
---|
2 | Based on code by Uwe Hollerbach <uh@alumni.caltech.edu> |
---|
3 | *) |
---|
4 | |
---|
5 | structure RungeKuttaTest = |
---|
6 | struct |
---|
7 | |
---|
8 | open RungeKutta |
---|
9 | |
---|
10 | val summer = Real.+ |
---|
11 | val scaler = Real.* |
---|
12 | |
---|
13 | infix 7 */ |
---|
14 | infix 6 +/ |
---|
15 | infix // |
---|
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 | |
---|
25 | val con = ~0.4 |
---|
26 | fun deriv (t,y) = con*y |
---|
27 | val t0 = 0.0 |
---|
28 | val y0 = 1.75 |
---|
29 | fun exact t = y0*Real.Math.exp(con*(t - t0)) |
---|
30 | |
---|
31 | fun putStr str = |
---|
32 | (TextIO.output (TextIO.stdOut, str)) |
---|
33 | |
---|
34 | fun putStrLn str = |
---|
35 | (TextIO.output (TextIO.stdOut, str); |
---|
36 | TextIO.output (TextIO.stdOut, "\n")) |
---|
37 | |
---|
38 | fun showReal n = Real.toString n |
---|
39 | |
---|
40 | fun showst (t, y) = String.concat [(showReal y), "\t", (showReal (y - (exact t)))] |
---|
41 | |
---|
42 | fun 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 | |
---|
52 | fun 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 | |
---|
62 | fun 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 | |
---|
72 | fun 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 | |
---|
79 | fun 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 | |
---|
89 | fun 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 | |
---|
97 | val rkfe: real stepper1 = make_rkfe() |
---|
98 | val rk3: real stepper1 = make_rk3() |
---|
99 | val rk4a: real stepper1 = make_rk4a() |
---|
100 | val rk4b: real stepper1 = make_rk4b() |
---|
101 | |
---|
102 | val rkhe: real stepper2 = make_rkhe() |
---|
103 | val rkbs: real stepper2 = make_rkbs() |
---|
104 | val rkf45: real stepper2 = make_rkf45() |
---|
105 | val rkck: real stepper2 = make_rkck() |
---|
106 | val rkdp: real stepper2 = make_rkdp() |
---|
107 | val rkf78: real stepper2 = make_rkf78() |
---|
108 | |
---|
109 | val rkhe_aux: real stepper1 = make_rkhe_aux() |
---|
110 | val rkbs_aux: real stepper1 = make_rkbs_aux() |
---|
111 | val rkf45_aux: real stepper1 = make_rkf45_aux() |
---|
112 | val rkck_aux: real stepper1 = make_rkck_aux() |
---|
113 | val rkdp_aux: real stepper1 = make_rkdp_aux() |
---|
114 | val rkf78_aux: real stepper1 = make_rkf78_aux() |
---|
115 | |
---|
116 | |
---|
117 | |
---|
118 | fun run() = |
---|
119 | (putStrLn "#### Non-Adaptive Solvers"; |
---|
120 | List.app solver1 [(rkfe, show_rkfe), |
---|
121 | (rk3, show_rk3), |
---|
122 | (rk4a, show_rk4a), |
---|
123 | (rk4b, show_rk4b)]; |
---|
124 | putStrLn "#### Adaptive Solvers"; |
---|
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 | |
---|
140 | val _ = run() |
---|
141 | end |
---|