source: project/release/4/flsim/trunk/sml-lib/rk/rk.sml @ 30946

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

flsim: added rkv65 method to rk library

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