source: project/release/4/flsim/trunk/sml-lib/rk/rk.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: 21.4 KB
Line 
1(*
2 * Runge-Kutta integration of ODEs
3 * Based on Haskell code by Uwe Hollerbach <uh@alumni.caltech.edu>
4 *
5 * dy
6 * -- = f(t, y)
7 * dt
8 *
9 * y_{n+1} = y_n + h sum_{i=1}^s b_i k_i
10 * k_i = f(t_n + c_i h, y_n + h sum_{j=1}^s a_{ij} k_j)
11 * * "Butcher Tableau" is
12 *
13 * c_1  a_11 a_12 ... a_1s
14 * c_2  a_21 a_22 ... a_2s
15 * ...  ...
16 * c_s  a_s1 a_s2 ... a_ss
17 *      b_1  b_2  ... b_s
18 *
19 * This module implements a method that can do a generic tableau, then
20 * specializes with different tableaux to let the user pick a specific
21 * method. Adaptive step-size methods, see below, add a row of d_j
22 * coefficients and use that to report the error:
23 *
24 * e_{n+1} = h sum_{i=1}^s d_i k_i
25 *
26 * non-adaptive solvers:
27 *      rkfe, rk3, rk4a, rk4b
28 *
29 * adaptive solvers:
30 *      rkhe, rkbs, rkf45, rkck, rkdp, rkf78, rkv65
31 *
32 * adaptive solvers with interpolation (CERK):
33 *      cerkdp
34 *
35 * auxiliary non-adaptive solvers (error estimators from the adaptive ones):
36 *      rkhe_aux, rkbs_aux, rkf45_aux, rkck_aux, rkdp_aux, rkf78_aux, rkv65_aux
37 *
38 * use rk4[ab] if you don't need an adaptive solver, rkdp or rkv65 if you do;
39 * or use what you need if you're an expert.
40 *
41 * DO NOT USE rkfe EXCEPT FOR DEMONSTRATIONS OF INSTABILITY!
42 * (Or if you're an expert.)
43 *
44 * Reference: E. Hairer, S. P. Norsett, G. Wanner,
45 * Solving Ordinary Differential Equations I: Nonstiff Problems
46 * (second revised edition, 1993).
47*)
48
49structure RungeKutta =
50struct
51
52
53exception InsufficientArguments
54exception InvalidCoefficients
55
56fun putStr str =
57    (TextIO.output (TextIO.stdOut, str))
58
59fun putStrLn str =
60    (TextIO.output (TextIO.stdOut, str);
61     TextIO.output (TextIO.stdOut, "\n"))
62
63
64fun foldl1 f (a::b::lst) = List.foldl f (f(a,b)) lst
65  | foldl1 f (a::[]) = a
66  | foldl1 f _ = raise InsufficientArguments
67
68fun list_show (toString,sep,lb,rb) xs =
69    let 
70        fun loop (x::xs,str) = loop (xs, str ^ sep ^ (toString x))
71          | loop ([],str) = str
72
73        val ls = case xs of x::xs => loop (xs, toString x) | [] => ""
74    in
75        lb ^ ls ^ rb
76    end
77
78fun def_list_show toString xs =
79    list_show (toString,",","[","]") xs
80
81(* Rational number implementation.
82 * Copyright (C) 2007 Vesa Karvonen
83 *
84 * This code is released under the MLton license, a BSD-style license.
85 * See the LICENSE file or http://mlton.org/License for details.
86 *)
87
88
89(* Silly implementation of rational numbers
90 *
91 * HINT: Someone should really implement a nice lib for rational numbers!
92 *)
93datatype rational =
94   RAT of LargeInt.int
95 | //  of LargeInt.int * LargeInt.int
96
97infix 7 */
98infix 6 +/
99infix //
100
101fun rat_show (RAT i) = LargeInt.toString i
102  | rat_show (n // d) = String.concat [(LargeInt.toString n), "//", (LargeInt.toString d)]
103
104val fromRational = 
105    fn (RAT n) => Real.fromLargeInt n
106     | (n // d) => Real./ (Real.fromLargeInt n, Real.fromLargeInt d)
107   
108val numerator =
109    fn (RAT n) => n
110     | (n // d) => n
111
112val denominator =
113    fn (RAT n) => 1
114     | (n // d) => d
115
116fun gcd (a, b) : LargeInt.int = if 0 = b then a else gcd (b, a mod b)
117
118val normalize =
119 fn 0 // _     => RAT 0
120  | r as RAT _ => r
121  | n // d     => let
122       val c = gcd (n, d)
123    in
124       if c = d
125       then RAT (n div c)
126       else n div c // d div c
127    end
128
129val op +/ = let
130   fun sym i n d = n + i * d // d
131in
132   fn (RAT l,   RAT r) => RAT (l + r)
133    | (RAT i,  n // d) => sym i n d
134    | (n // d,  RAT i) => sym i n d
135    | (n // d, m // e) =>
136      normalize (if d = e then n + m // d else n*e + m*d // d*e)
137end
138
139val op */ = let
140   fun sym i n d = normalize (i*n // d)
141in
142   fn (RAT l,   RAT r) => RAT (l * r)
143    | (RAT i,  n // d) => sym i n d
144    | (n // d,  RAT i) => sym i n d
145    | (n // d, m // e) => normalize (n*m // d*e)
146end
147
148
149(* Store a list of rational numbers as a common denominator, then a list
150   of numerators, all stored as doubles: this lets us write the values in
151   the source as exact rational numbers, yet we can take advantage of not
152   having to do too many conversions and also of not having to multiply
153   by 0 or 1 in some cases.
154*)
155
156type RCL = real * real list
157
158(* ratToRCL :: [rational] -> RCL *)
159fun ratToRCL [] = (1.0, [])
160  | ratToRCL rs =
161    let 
162        val ds = map denominator rs
163        val dp = foldl1 ( op * ) ds
164        val ns = map (numerator o (fn (x) => (RAT dp) */ x)) rs
165        val g  = foldl gcd dp ns
166    in 
167        (Real.fromLargeInt (LargeInt.quot (dp, g)), 
168         map (Real.fromLargeInt o (fn x => LargeInt.quot (x,g))) ns)
169    end
170
171(* ratToRCLs :: [[rational]] -> [RCL] *)
172val ratToRCLs = fn x => map ratToRCL x
173
174(* ratToReals :: [rational] -> [real] *)
175val ratToReals = fn x => map fromRational x
176
177(* Helper function to sum a list of K_i, skipping
178   un-necessary multiplications and additions *)
179
180fun k_sum (sc_fn: real * 'a -> 'a, 
181           sum_fn: 'a * 'a -> 'a, 
182           h: real) 
183          ((d,ns), ks) =
184    let fun m_scale (s,v) =
185             (if Real.compare(s,0.0) = EQUAL then NONE
186              else (if Real.compare(s,1.0) = EQUAL then SOME v
187                    else SOME (sc_fn (s,v))))
188        val ns_ks = ListPair.zip (ns,ks)
189    in
190        sc_fn (Real./ (h,d), foldl1 sum_fn (List.mapPartial m_scale ns_ks  ))
191    end
192
193
194(* Helper function to generate a list of K_i *)
195
196fun gen_ks (ksum_fn,sum_fn: 'a * 'a -> 'a,der_fn: real * 'a -> 'a,
197            h,(tn,yn),ks,[],[]) = ks
198
199  | gen_ks (ksum_fn,sum_fn,der_fn,h,old as (tn,yn),ks,(c::cs),(a::ar)) =
200    let
201        val yn1 = if (List.null ks) then yn else sum_fn (yn, ksum_fn (a,ks))
202    in
203        gen_ks (ksum_fn, sum_fn, der_fn, h, old, (ks @ [der_fn (tn + c*h, yn1)]), cs, ar)
204    end
205
206  | gen_ks (ksum_fn,sum_fn,der_fn,h,(tn,yn),ks,_,_) =
207    raise InvalidCoefficients
208   
209
210(*
211
212   This is the first core routine: it does not get exported,
213   only partial applications of it; see below.
214
215   Its arguments are:
216
217     c table (specified internally)
218     a table (specified internally)
219     b table (specified internally)
220
221   user-specified arguments:
222
223     scale function to scale a Y state vector ::
224        (real * a -> a)
225
226     sum function to add two Y state vectors ::
227        (a * a -> a)
228
229     derivative function F ::
230        (real * a -> a)
231
232     step size H ::
233        real
234
235     current state (T,Y) ::
236        (real, a)
237
238     and the return value is the new state Y_new
239*)
240
241type 'a stepper1 =  ((real * 'a -> 'a) *
242                     ('a * 'a -> 'a)   *
243                     (real * 'a -> 'a)) ->
244                     (real -> (real * 'a) -> 'a)
245
246fun core1
247        (cl: real list, al: RCL list, bl: RCL)
248        (sc_fn:  real * 'a -> 'a, 
249         sum_fn: 'a * 'a -> 'a,
250         der_fn: real * 'a -> 'a)
251        (h: real) 
252        (old as (tn,yn: 'a)) =
253    (let
254         val ksum = k_sum (sc_fn,sum_fn,h)
255         val ks   = gen_ks (ksum, sum_fn, der_fn, h, old, [], cl, al)
256     in
257         sum_fn (yn, ksum (bl, ks))
258     end)
259
260
261(*
262   This is the second core routine, analogous to the previous one.
263   The difference is that this gets an additional internal table arg,
264   and it returns a tuple instead of a single value: (ynew,enew),
265   where enew is the error state vector
266   e_{n+1} = h sum_{i=1}^s (b_i - b'_i) k_i
267*)
268
269type 'a stepper2 =  ((real * 'a -> 'a) *
270                     ('a * 'a -> 'a)   *
271                     (real * 'a -> 'a)) ->
272                    (real -> (real * 'a) -> ('a * 'a))
273
274fun core2 (cl: real list, al: RCL list, bl: RCL, dl: RCL) 
275          (sc_fn: real * 'a -> 'a, 
276           sum_fn: 'a * 'a -> 'a,
277           der_fn: real * 'a -> 'a)
278           (h: real)
279           (old as (tn,yn: 'a)) =
280  let
281      val ksum = k_sum (sc_fn,sum_fn,h)
282      val ks   = gen_ks (ksum, sum_fn, der_fn, h, old, [], cl, al)
283  in
284      (sum_fn (yn, ksum (bl, ks)), ksum (dl, ks))
285  end
286
287
288
289(* Helper function to sum a list of b_i (theta) K_i *)
290
291fun bk_sum (bs: RCL list)
292           (sc_fn: real * 'a -> 'a, sum_fn: 'a * 'a -> 'a) 
293           (ks: 'a list, h: real)
294           (theta: real) = 
295    let fun m_scale (s,v) =
296            (if Real.compare(s,0.0) = EQUAL then NONE
297             else (if Real.compare(s,1.0) = EQUAL then SOME v
298                   else SOME (sc_fn (s,v))))
299
300        fun recur ((d,ns)::bs, k::ks, fs) =
301            let
302                val (bsum,_) = foldl (fn (n,(sum,theta)) => 
303                                         ((n*theta)+sum,theta*theta)) (0.0,theta) ns
304            in
305                case m_scale (bsum, k) of 
306                    SOME bk => recur (bs, ks, (sc_fn (h/d, bk))::fs)
307                  | NONE => recur (bs, ks, fs)
308            end
309          | recur ([], [k], fs) = foldl1 sum_fn fs
310          | recur (_, _, _) = raise InvalidCoefficients
311    in
312        recur (bs, ks, [])
313    end
314   
315
316(* Interpolation routine for continuous explicit RK (CERK) methods *)
317fun interp (ws: RCL list)
318           (sc_fn: real * 'a -> 'a, 
319            sum_fn: 'a * 'a -> 'a)
320           (ks: 'a list, h: real)
321           (tn,yn: 'a) (theta: real) = 
322    sum_fn (yn, bk_sum ws (sc_fn,sum_fn) (ks,h) theta)
323
324
325(* Core routine for constructing continuous methods.  It returns a
326   triple (ynew,enew,interp), where interp is the interpolation
327   function for this timestep.
328*)
329
330type 'a stepper3 =  ((real * 'a -> 'a) *
331                     ('a * 'a -> 'a)   *
332                     (real * 'a -> 'a)) ->
333                    (real -> (real * 'a) -> ('a * 'a * (real -> 'a)))
334
335fun core3 (cl: real list, al: RCL list, bl: RCL, dl: RCL, wl: RCL list) 
336          (sc_fn: real * 'a -> 'a, 
337           sum_fn: 'a * 'a -> 'a,
338           der_fn: real * 'a -> 'a)
339           (h: real)
340           (old as (tn,yn: 'a)) =
341  let
342      val interp'   = interp wl (sc_fn,sum_fn)
343      val ksum      = k_sum (sc_fn,sum_fn,h)
344      val ks        = gen_ks (ksum, sum_fn, der_fn, h, old, [], cl, al)
345  in
346      (sum_fn (yn, ksum (bl, ks)),
347       ksum (dl, ks),
348       interp' (ks, h) (tn,yn))
349  end
350
351(* Helper routines to show the internal tables *)
352
353fun rcl_show (d,ns) =
354    "<" ^ (Real.toString d) ^ ", " ^ (def_list_show Real.toString ns) ^ ">"
355
356
357fun rk_show1 (title,cs,ar: RCL list,bs) =
358    title ^ ":\ncs:\t" ^ ((def_list_show Real.toString) cs) ^
359    "\nas:\t" ^ (list_show (rcl_show,"\n\t","","") ar) ^
360    "\nbs:\t" ^ (rcl_show bs) 
361
362fun rk_show2 (title,cs,ar: RCL list,bs,ds) =
363    title ^ ":\nds:\t" ^ (rcl_show ds) ^
364    "\ncs:\t" ^ ((def_list_show Real.toString) cs) ^
365    "\nbs:\t" ^ (rcl_show bs) ^
366    "\nas:\t" ^ (list_show (rcl_show,"\n\t","","") ar) 
367
368fun rk_show3 (title,cs,ar: RCL list,bs,ds,ws) =
369    title ^ ":ws:\t" ^ (list_show (rcl_show,"\n\t","","") ws) ^
370    "\nds:\t" ^ (rcl_show ds) ^
371    "\ncs:\t" ^ ((def_list_show Real.toString) cs) ^
372    "\nbs:\t" ^ (rcl_show bs) ^
373    "\nas:\t" ^ (list_show (rcl_show,"\n\t","","") ar) 
374
375
376(*
377   Some specific explicit methods, taken from
378   "List of Runge-Kutta methods" at Wikipedia
379*)
380
381(* forward Euler: unconditionally unstable: don't use this! *)
382
383val cs_fe = ratToReals [RAT 0]
384val as_fe = ratToRCLs [[]]
385val bs_fe = ratToRCL  [RAT 1]
386fun make_rkfe (): 'a stepper1 = core1 (cs_fe, as_fe, bs_fe)
387val show_rkfe = rk_show1 ("Forward Euler", cs_fe, as_fe, bs_fe)
388
389(* Kutta's third-order method: *)
390
391val cs_rk3 = ratToReals [RAT 0, 1//2, RAT 1]
392val as_rk3 = ratToRCLs [[], [1//2], [RAT ~1, RAT 2]]
393val bs_rk3 = ratToRCL  [1//6, 2//3, 1//6]
394fun make_rk3 (): 'a stepper1 = core1 (cs_rk3, as_rk3, bs_rk3)
395val show_rk3 = rk_show1 ("Kutta's third-order method", cs_rk3, as_rk3, bs_rk3)
396
397(* Classic fourth-order method *)
398
399val cs_rk4a = ratToReals [RAT 0, 1//2, 1//2, RAT 1]
400val as_rk4a = ratToRCLs [[], [1//2], [RAT 0, 1//2], [RAT 0, RAT 0, RAT 1]]
401val bs_rk4a = ratToRCL  [1//6, 1//3, 1//3, 1//6]
402fun make_rk4a (): 'a stepper1 = core1 (cs_rk4a, as_rk4a, bs_rk4a)
403val show_rk4a = rk_show1 ("Classic fourth-order method", cs_rk4a, as_rk4a, bs_rk4a)
404
405(* Kutta's other fourth-order method... "The first [above] is more popular,
406   the second is more precise." (Hairer, Norsett, Wanner) *)
407
408val cs_rk4b = ratToReals [RAT 0, 1//3, 2//3, RAT 1]
409val as_rk4b = ratToRCLs [[], [1//3], [~1//3, RAT 1], [RAT 1, RAT ~1, RAT 1]]
410val bs_rk4b = ratToRCL  [1//8, 3//8, 3//8, 1//8]
411fun make_rk4b (): 'a stepper1 = core1 (cs_rk4b, as_rk4b, bs_rk4b)
412val show_rk4b = rk_show1 ("Kutta's other classic fourth-order method", cs_rk4b, as_rk4b, bs_rk4b)
413
414(*
415   Some adaptive-stepsize methods, also from Wikipedia; more from HNW.
416   These don't auto-adapt, but they do allow the user to make a somewhat
417   intelligent decision about what the step size ought to be at each step.
418*)
419
420(*
421   A helper function to take the difference of two lists of rationals:
422   we don't want to use zipWith (-) because that gives only the head
423   where both lists have entries; we want implicit zeros at the end,
424   as far as is necessary.
425*)
426
427fun negate x = (RAT ~1) */ x
428
429fun diffs ([], []) = []
430  | diffs (xs, []) = xs
431  | diffs ([], xs) = map negate xs
432  | diffs (x::xs, y::ys) = (x +/ (negate y)) :: (diffs (xs,ys))
433
434(* Heun-Euler, order 2/1 *)
435
436val cs_he = ratToReals [RAT 0, RAT 1]
437val as_he = ratToRCLs [[], [RAT 1]]
438val r1_he = [1//2, 1//2]   (* second-order coeffs *)
439val r2_he = [RAT 1]        (* first-order coeffs *)
440val bs_he = ratToRCL r1_he
441val ds_he = ratToRCL (diffs (r1_he, r2_he))
442fun make_rkhe (): 'a stepper2  = core2 (cs_he, as_he, bs_he, ds_he)
443val show_rkhe = rk_show2 ("Heun-Euler 2(1)", cs_he, as_he, bs_he, ds_he)
444
445val bs_he_aux = ratToRCL r2_he
446fun make_rkhe_aux (): 'a stepper1 = core1 (cs_he, as_he, bs_he_aux)
447val show_rkhe_aux = rk_show1 ("Heun-Euler (1)", cs_he, as_he, bs_he_aux)
448
449(* Bogacki-Shampine, order 3/2 *)
450
451val cs_bs = ratToReals [RAT 0, 1//2, 3//4, RAT 1]
452val as_bs = ratToRCLs [[], [1//2], [RAT 0, 3//4], [2//9, 1//3, 4//9]]
453val r1_bs = [2//9, 1//3, 4//9]          (* third-order coeffs *)
454val r2_bs = [7//24, 1//4, 1//3, 1//8]   (* second-order coeffs *)
455val bs_bs = ratToRCL r1_bs
456val ds_bs = ratToRCL (diffs (r1_bs, r2_bs))
457fun make_rkbs (): 'a stepper2 = core2 (cs_bs, as_bs, bs_bs, ds_bs)
458val show_rkbs = rk_show2 ("Bogacki-Shampine 3(2)", cs_bs, as_bs, bs_bs, ds_bs)
459
460val bs_bs_aux = ratToRCL r2_bs
461fun make_rkbs_aux (): 'a stepper1 = core1 (cs_bs, as_bs, bs_bs_aux)
462val show_rkbs_aux = rk_show1 ("Bogacki-Shampine (2)", cs_bs, as_bs, bs_bs_aux)
463
464(* Runge-Kutta-Fehlberg, order 4/5 *)
465
466val cs_rkf = ratToReals [RAT 0, 1//4, 3//8, 12//13, RAT 1, 1//2]
467val as_rkf = ratToRCLs [[],
468                        [1//4],
469                        [3//32, 9//32],
470                        [1932//2197, ~7200//2197, 7296//2197],
471                        [439//216, RAT ~8, 3680//513, ~845//4104],
472                        [~8//27, RAT 2, ~3544//2565, 1859//4104, ~11//40]]
473(* fourth-order coeffs *)
474val r1_rkf = [25//216, RAT 0, 1408//2565, 2197//4104, ~1//5]
475(* fifth-order coeffs *)
476val r2_rkf = [16//135, RAT 0, 6656//12825, 28561//56430, ~9//50, 2//55]
477val bs_rkf = ratToRCL r1_rkf
478val ds_rkf = ratToRCL (diffs (r1_rkf, r2_rkf))
479fun make_rkf45 (): 'a stepper2  = core2 (cs_rkf, as_rkf, bs_rkf, ds_rkf)
480val show_rkf45 = rk_show2 ("Runge-Kutta-Fehlberg 4(5)", cs_rkf, as_rkf, bs_rkf, ds_rkf)
481
482val bs_rkf_aux = ratToRCL r2_rkf
483fun make_rkf45_aux (): 'a stepper1 = core1 (cs_rkf, as_rkf, bs_rkf_aux)
484val show_rkf45_aux = rk_show1 ("Runge-Kutta-Fehlberg (5)", cs_rkf, as_rkf, bs_rkf_aux)
485
486(* Cash-Karp, order 4/5 (use 4th-order sol'n,
487   coeffs chosen to minimize error of 4th-order sol'n) *)
488
489val cs_ck = ratToReals [RAT 0, 1//5, 3//10, 3//5, RAT 1, 7//8]
490val as_ck = ratToRCLs [[],
491                       [1//5],
492                       [3//40, 9//40],
493                       [3//10, ~9//10, 6//5],
494                       [~11//54, 5//2, ~70//27, 35//27],
495                       [1631//55296, 175//512, 575//13824, 44275//110592, 253//4096]]
496(* fourth-order coeffs *)
497val r1_ck = [2825//27648, RAT 0, 18575//48384, 13525//55296, 277//14336, 1//4]
498(* fifth-order coeffs *)
499val r2_ck = [37//378, RAT 0, 250//621, 125//594, RAT 0, 512//1771]
500val bs_ck = ratToRCL r1_ck
501val ds_ck = ratToRCL (diffs (r1_ck, r2_ck))
502fun make_rkck (): 'a stepper2  = core2 (cs_ck, as_ck, bs_ck, ds_ck)
503val show_rkck = rk_show2 ("Cash-Karp 4(5)", cs_ck, as_ck, bs_ck, ds_ck)
504
505val bs_ck_aux = ratToRCL r2_ck
506fun make_rkck_aux (): 'a stepper1 = core1 (cs_ck, as_ck, bs_ck_aux)
507val show_rkck_aux = rk_show1 ("Cash-Karp (5)", cs_ck, as_ck, bs_ck_aux)
508
509(* Dormand-Prince, order 5/4 (use 5th-order sol'n, coeffs chosen to
510   minimize error of 5th-order sol'n) This is DOPRI5 from Hairer,
511   Norsett, Wanner
512*)
513
514val cs_dp = ratToReals [RAT 0, 1//5, 3//10, 4//5, 8//9, RAT 1, RAT 1]
515val as_dp = ratToRCLs [[],
516                       [1//5],
517                       [3//40, 9//40],
518                       [44//45, ~56//15, 32//9],
519                       [19372//6561, ~25360//2187, 64448//6561, ~212//729],
520                       [9017//3168, ~355//33, 46732//5247, 49//176, ~5103//18656],
521                       [35//384, RAT 0, 500//1113, 125//192, ~2187//6784, 11//84]]
522(* fifth-order coeffs *)
523val r1_dp = [35//384, RAT 0, 500//1113, 125//192, ~2187//6784, 11//84]
524(* fourth-order coeffs *)
525val r2_dp = [5179//57600, RAT 0, 7571//16695, 393//640, ~92097//339200, 187//2100, 1//40]
526val bs_dp = ratToRCL r1_dp
527val ds_dp = ratToRCL (diffs (r1_dp, r2_dp))
528(* interpolation coeffs for continuous method *)
529val ws_dp = ratToRCLs [[RAT 1, ~1337//480, 1039//360, ~1163//1152],
530                       [],
531                       [RAT 0, 4216//1113, ~18728//3339, 7580//3339],
532                       [RAT 0, ~27//16, 9//2, ~415//192],
533                       [RAT 0, ~2187//8480, 2673//2120, ~8991//6784],
534                       [RAT 0, 33//35, ~319//105, 187//84]]
535                                                         
536
537fun make_cerkdp (): 'a stepper3  = core3 (cs_dp, as_dp, bs_dp, ds_dp, ws_dp)
538val show_cerkdp = rk_show3 ("Dormand-Prince 5(4)", cs_dp, as_dp, bs_dp, ds_dp, ws_dp)
539
540fun make_rkdp (): 'a stepper2  = core2 (cs_dp, as_dp, bs_dp, ds_dp)
541val show_rkdp = rk_show2 ("Dormand-Prince 5(4) \"DOPRI5\"", cs_dp, as_dp, bs_dp, ds_dp)
542
543val bs_dp_aux = ratToRCL r2_dp
544fun make_rkdp_aux (): 'a stepper1 = core1 (cs_dp, as_dp, bs_dp_aux)
545val show_rkdp_aux = rk_show1 ("Dormand-Prince (4)", cs_dp, as_dp, bs_dp_aux)
546
547
548(* Alternative coefficients for Dormand-Prince, from Butcher p. 195.
549   "Although this has larger error constants overall ... it has the
550   advantage of a longer stability interval [than DOPRI5]."
551*)
552
553val cs_dpb = ratToReals [RAT 0, 2//9, 1//3, 5//9, 2//3, RAT 1, RAT 1]
554val as_dpb = ratToRCLs [[],
555                        [2//9],
556                        [1//12, 1//4],
557                        [55//324, ~25//108, 50//81],
558                        [83//330, ~13//22, 61//66, 9//110],
559                        [~19//28, 9//4, 1//7, ~27//7, 22//7],
560                        [19//200, RAT 0, 3//5, ~243//400, 33//40, 7//80]]
561(* fifth-order coeffs *)
562val r1_dpb = [19//200, RAT 0, 3//5, ~243//400, 33//40, 7//80]
563(* fourth-order coeffs *)
564val r2_dpb = [431//5000, RAT 0, 333//500, ~7857//10000, 957//1000, 193//2000, ~1//50]
565val bs_dpb = ratToRCL r1_dpb
566val ds_dpb = ratToRCL (diffs (r1_dpb, r2_dpb))
567fun make_rkdpb (): 'a stepper2  = core2 (cs_dpb, as_dpb, bs_dpb, ds_dpb)
568val show_rkdpb = rk_show2 ("Dormand-Prince 5(4) B", cs_dpb, as_dpb, bs_dpb, ds_dpb)
569
570val bs_dpb_aux = ratToRCL r2_dpb
571fun make_rkdpb_aux (): 'a stepper1 = core1 (cs_dpb, as_dpb, bs_dpb_aux)
572val show_rkdpb_aux = rk_show1 ("Dormand-Prince (4) B", cs_dpb, as_dpb, bs_dpb_aux)
573
574(*
575   Fehlberg, order 7/8: "... of frequent use in all high precision
576   computations, e.g., in astronomy."  --Hairer, Norsett, Wanner.
577   But caveat: suffers from the drawback that error estimate is
578   identically 0 for quadrature problems y' = f(x)
579
580   NOTE BUG in Hairer Norsett Wanner: the third-last A coefficient in the
581   last row of the tableau is listed as "19/41" in the book. This is WRONG:
582   that row does not then sum to 1, and the convergence of the auxiliary
583   solver is then order 1 or 2, not 8
584*)
585
586val cs_f78 = ratToReals [RAT 0, 2//27, 1//9, 1//6, 5//12, 1//2, 5//6, 1//6, 2//3, 1//3, RAT 1, RAT 0, RAT 1]
587(* TODO re-check this *)
588val as_f78 = ratToRCLs
589                 [[],
590                  [2//27],
591                  [1//36, 1//12],
592                  [1//24, RAT 0, 1//8],
593                  [5//12, RAT 0, ~25//16, 25//16],
594                  [1//20, RAT 0, RAT 0, 1//4, 1//5],
595                  [~25//108, RAT 0, RAT 0, 125//108, ~65//27, 125//54],
596                  [31//300, RAT 0, RAT 0, RAT 0, 61//225, ~2//9, 13//900],
597                  [RAT 2, RAT 0, RAT 0, ~53//6, 704//45, ~107//9, 67//90, RAT 3],
598                  [~91//108, RAT 0, RAT 0, 23//108, ~976//135, 311//54, ~19//60, 17//6, ~1//12],
599                  [2383//4100, RAT 0, RAT 0, ~341//164, 4496//1025, ~301//82, 2133//4100, 45//82, 45//164, 18//41],
600                  [3//205, RAT 0, RAT 0, RAT 0, RAT 0, ~6//41, ~3//205, ~3//41, 3//41, 6//41, RAT 0],
601                  [~1777//4100, RAT 0, RAT 0, ~341//164, 4496//1025, ~289//82, 2193//4100, 51//82, 33//164, 12//41, RAT 0, RAT 1]]
602(* seventh-order coeffs *)
603val r1_f78 = [41//840, RAT 0, RAT 0, RAT 0, RAT 0, 34//105, 9//35, 9//35, 9//280, 9//280, 41//840]
604(* eighth-order coeffs *)
605val r2_f78 = [RAT 0, RAT 0, RAT 0, RAT 0, RAT 0, 34//105, 9//35, 9//35, 9//280, 9//280, RAT 0, 41//840, 41//840]
606val bs_f78 = ratToRCL r1_f78
607val ds_f78 = ratToRCL (diffs (r1_f78, r2_f78))
608fun make_rkf78 (): 'a stepper2 = core2 (cs_f78, as_f78, bs_f78, ds_f78)
609val show_rkf78 = rk_show2 ("Fehlberg 7(8)", cs_f78, as_f78, bs_f78, ds_f78)
610
611val bs_f78_aux = ratToRCL r2_f78
612fun make_rkf78_aux (): 'a stepper1 = core1 (cs_f78, as_f78, bs_f78_aux)
613val show_rkf78_aux = rk_show1 ("Fehlberg (8)", cs_f78, as_f78, bs_f78_aux)
614
615
616(* Verner, order 6/5 DVERK *)
617
618val cs_v65 = ratToReals [RAT 0, 1//6, 4//15, 2//3, 5//6, RAT 1, RAT 115, RAT 1]
619val as_v65 = ratToRCLs
620                 [[],
621                  [1//6],
622                  [4//75, 16//75],
623                  [5//6, ~8//3, 5//2],
624                  [~165//64, 55//6, ~425//64, 85//96],
625                  [12//5, RAT ~8, 4015//612, ~11//36, 88//255],
626                  [~8263//15000, 124//75, ~643//680, ~81//250, 2484//10625],
627                  [3501//1720, ~300//43, 297275//52632, ~319//2322, 24068//84065, RAT 0, 3850//26703]]
628(* sixth-order coeffs *)
629val r1_v65 = [3//40, RAT 0, 875//2244, 23//72, 264//1955, RAT 0, 125//11592, 43//616]
630(* fifth-order coeffs *)
631val r2_v65 = [13//160, RAT 0, 2375//5984, 5//16, 12//85, 3//44]
632val bs_v65 = ratToRCL r1_v65
633val ds_v65 = ratToRCL (diffs (r1_v65, r2_v65))
634fun make_rkv65 (): 'a stepper2 = core2 (cs_v65, as_v65, bs_v65, ds_v65)
635val show_rkv65 = rk_show2 ("Verner 6(5) \"DVERK\"", cs_v65, as_v65, bs_v65, ds_v65)
636
637val bs_v65_aux = ratToRCL r2_v65
638fun make_rkv65_aux (): 'a stepper1 = core1 (cs_v65, as_v65, bs_v65_aux)
639val show_rkv65_aux = rk_show1 ("Verner (5)", cs_v65, as_v65, bs_v65_aux)
640
641end
Note: See TracBrowser for help on using the repository browser.