Changeset 30954 in project


Ignore:
Timestamp:
06/03/14 08:24:47 (6 years ago)
Author:
Ivan Raikov
Message:

flsim: update to SML RK library

Location:
release/4/flsim/trunk/sml-lib/rk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • release/4/flsim/trunk/sml-lib/rk/rk.sml

    r30946 r30954  
    3030 *      rkhe, rkbs, rkf45, rkck, rkdp, rkf78, rkv65
    3131 *
     32 * adaptive solvers with interpolation (CERK):
     33 *      cerkdp
     34 *
    3235 * auxiliary non-adaptive solvers (error estimators from the adaptive ones):
    3336 *      rkhe_aux, rkbs_aux, rkf45_aux, rkck_aux, rkdp_aux, rkf78_aux, rkv65_aux
     
    4750struct
    4851
    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 *)
    6152
    6253exception InsufficientArguments
     
    294285  end
    295286
     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
    296351(* Helper routines to show the internal tables *)
    297 
    298352
    299353fun rcl_show (d,ns) =
     
    308362fun rk_show2 (title,cs,ar: RCL list,bs,ds) =
    309363    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) ^
    310371    "\ncs:\t" ^ ((def_list_show Real.toString) cs) ^
    311372    "\nbs:\t" ^ (rcl_show bs) ^
     
    465526val bs_dp = ratToRCL r1_dp
    466527val 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
    467540fun make_rkdp (): 'a stepper2  = core2 (cs_dp, as_dp, bs_dp, ds_dp)
    468541val show_rkdp = rk_show2 ("Dormand-Prince 5(4) \"DOPRI5\"", cs_dp, as_dp, bs_dp, ds_dp)
     
    471544fun make_rkdp_aux (): 'a stepper1 = core1 (cs_dp, as_dp, bs_dp_aux)
    472545val 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)
    473573
    474574(*
  • release/4/flsim/trunk/sml-lib/rk/rktest.sml

    r22665 r30954  
    9494   putStrLn "# All done!\n")
    9595
     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
    96124
    97125val rkfe: real stepper1 = make_rkfe()
     
    105133val rkck:  real stepper2 = make_rkck()
    106134val rkdp:  real stepper2 = make_rkdp()
     135val rkdpb:  real stepper2 = make_rkdpb()
    107136val rkf78: real stepper2 = make_rkf78()
     137val rkv65: real stepper2 = make_rkv65()
    108138
    109139val rkhe_aux:  real stepper1  = make_rkhe_aux()
     
    112142val rkck_aux:  real stepper1 = make_rkck_aux()
    113143val rkdp_aux:  real stepper1 = make_rkdp_aux()
     144val rkdpb_aux:  real stepper1 = make_rkdpb_aux()
    114145val rkf78_aux: real stepper1 = make_rkf78_aux()
     146val rkv65_aux: real stepper1 = make_rkv65_aux()
    115147
     148val cerkdp:  real stepper3 = make_cerkdp()
    116149
    117150
     
    128161                    (rkck, show_rkck),
    129162                    (rkdp, show_rkdp),
    130                     (rkf78, show_rkf78)];
     163                    (rkdpb, show_rkdpb),
     164                    (rkf78, show_rkf78),
     165                    (rkv65, show_rkv65)];
     166  putStrLn "#### Continuous Solvers";
     167  List.app solver3 [(cerkdp, show_cerkdp)];
    131168  putStrLn "#### Auxiliary Solvers: Error Estimators from Adaptives";
    132169  List.app solver1 [(rkhe_aux, show_rkhe_aux),
     
    135172                    (rkck_aux, show_rkck_aux),
    136173                    (rkdp_aux, show_rkdp_aux),
    137                     (rkf78_aux, show_rkf78_aux)]
     174                    (rkdpb_aux, show_rkdpb_aux),
     175                    (rkf78_aux, show_rkf78_aux),
     176                    (rkv65_aux, show_rkv65_aux)]
    138177  )
    139178
Note: See TracChangeset for help on using the changeset viewer.