Changeset 30954 in project
 Timestamp:
 06/03/14 08:24:47 (6 years ago)
 Location:
 release/4/flsim/trunk/smllib/rk
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

release/4/flsim/trunk/smllib/rk/rk.sml
r30946 r30954 30 30 * rkhe, rkbs, rkf45, rkck, rkdp, rkf78, rkv65 31 31 * 32 * adaptive solvers with interpolation (CERK): 33 * cerkdp 34 * 32 35 * auxiliary nonadaptive solvers (error estimators from the adaptive ones): 33 36 * rkhe_aux, rkbs_aux, rkf45_aux, rkck_aux, rkdp_aux, rkf78_aux, rkv65_aux … … 47 50 struct 48 51 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 52 62 53 exception InsufficientArguments … … 294 285 end 295 286 287 288 289 (* Helper function to sum a list of b_i (theta) K_i *) 290 291 fun 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 *) 317 fun 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 330 type 'a stepper3 = ((real * 'a > 'a) * 331 ('a * 'a > 'a) * 332 (real * 'a > 'a)) > 333 (real > (real * 'a) > ('a * 'a * (real > 'a))) 334 335 fun 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 296 351 (* Helper routines to show the internal tables *) 297 298 352 299 353 fun rcl_show (d,ns) = … … 308 362 fun rk_show2 (title,cs,ar: RCL list,bs,ds) = 309 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 368 fun 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) ^ 310 371 "\ncs:\t" ^ ((def_list_show Real.toString) cs) ^ 311 372 "\nbs:\t" ^ (rcl_show bs) ^ … … 465 526 val bs_dp = ratToRCL r1_dp 466 527 val ds_dp = ratToRCL (diffs (r1_dp, r2_dp)) 528 (* interpolation coeffs for continuous method *) 529 val 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 537 fun make_cerkdp (): 'a stepper3 = core3 (cs_dp, as_dp, bs_dp, ds_dp, ws_dp) 538 val show_cerkdp = rk_show3 ("DormandPrince 5(4)", cs_dp, as_dp, bs_dp, ds_dp, ws_dp) 539 467 540 fun make_rkdp (): 'a stepper2 = core2 (cs_dp, as_dp, bs_dp, ds_dp) 468 541 val show_rkdp = rk_show2 ("DormandPrince 5(4) \"DOPRI5\"", cs_dp, as_dp, bs_dp, ds_dp) … … 471 544 fun make_rkdp_aux (): 'a stepper1 = core1 (cs_dp, as_dp, bs_dp_aux) 472 545 val show_rkdp_aux = rk_show1 ("DormandPrince (4)", cs_dp, as_dp, bs_dp_aux) 546 547 548 (* Alternative coefficients for DormandPrince, 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 553 val cs_dpb = ratToReals [RAT 0, 2//9, 1//3, 5//9, 2//3, RAT 1, RAT 1] 554 val 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 (* fifthorder coeffs *) 562 val r1_dpb = [19//200, RAT 0, 3//5, ~243//400, 33//40, 7//80] 563 (* fourthorder coeffs *) 564 val r2_dpb = [431//5000, RAT 0, 333//500, ~7857//10000, 957//1000, 193//2000, ~1//50] 565 val bs_dpb = ratToRCL r1_dpb 566 val ds_dpb = ratToRCL (diffs (r1_dpb, r2_dpb)) 567 fun make_rkdpb (): 'a stepper2 = core2 (cs_dpb, as_dpb, bs_dpb, ds_dpb) 568 val show_rkdpb = rk_show2 ("DormandPrince 5(4) B", cs_dpb, as_dpb, bs_dpb, ds_dpb) 569 570 val bs_dpb_aux = ratToRCL r2_dpb 571 fun make_rkdpb_aux (): 'a stepper1 = core1 (cs_dpb, as_dpb, bs_dpb_aux) 572 val show_rkdpb_aux = rk_show1 ("DormandPrince (4) B", cs_dpb, as_dpb, bs_dpb_aux) 473 573 474 574 (* 
release/4/flsim/trunk/smllib/rk/rktest.sml
r22665 r30954 94 94 putStrLn "# All done!\n") 95 95 96 fun 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 107 fun 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 117 fun 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 96 124 97 125 val rkfe: real stepper1 = make_rkfe() … … 105 133 val rkck: real stepper2 = make_rkck() 106 134 val rkdp: real stepper2 = make_rkdp() 135 val rkdpb: real stepper2 = make_rkdpb() 107 136 val rkf78: real stepper2 = make_rkf78() 137 val rkv65: real stepper2 = make_rkv65() 108 138 109 139 val rkhe_aux: real stepper1 = make_rkhe_aux() … … 112 142 val rkck_aux: real stepper1 = make_rkck_aux() 113 143 val rkdp_aux: real stepper1 = make_rkdp_aux() 144 val rkdpb_aux: real stepper1 = make_rkdpb_aux() 114 145 val rkf78_aux: real stepper1 = make_rkf78_aux() 146 val rkv65_aux: real stepper1 = make_rkv65_aux() 115 147 148 val cerkdp: real stepper3 = make_cerkdp() 116 149 117 150 … … 128 161 (rkck, show_rkck), 129 162 (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)]; 131 168 putStrLn "#### Auxiliary Solvers: Error Estimators from Adaptives"; 132 169 List.app solver1 [(rkhe_aux, show_rkhe_aux), … … 135 172 (rkck_aux, show_rkck_aux), 136 173 (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)] 138 177 ) 139 178
Note: See TracChangeset
for help on using the changeset viewer.