Changeset 30939 in project


Ignore:
Timestamp:
05/29/14 11:14:27 (6 years ago)
Author:
Ivan Raikov
Message:

flsim: support for variable timestep integration

File:
1 edited

Legend:

Unmodified
Added
Removed
  • release/4/flsim/trunk/flsim.scm

    r30106 r30939  
    575575       ("val " ,solver ": (real list) stepper1 = make_" ,solver "()" ,nl)
    576576       ("fun make_stepper (deriv) = " ,solver " (scaler,summer,deriv)" ,nl)
    577        ("fun integrate (f,x: real,y: real list,h,i) = ((make_stepper f) h) (x,y)" ,nl)
    578        ))
    579 
     577       . ,(case solver 
     578            ;; adaptive solvers
     579            ((rkhe rkbs rkf45 rkck rkdp rkf78 rkv65)
     580             `(
     581#<<EOF
     582
     583val tol = Real.Math.pow (10.0, ~6.0)
     584val lb = 0.5 * tol
     585val ub = 1.0 * tol
     586
     587
     588datatype ('a, 'b) either = Left of 'a | Right of 'b
     589
     590
     591fun predictor tol (h,ys) =
     592  let open Real
     593      val e = foldl (fn (y,ax) => Real.+ ((abs y),ax)) 0.0 ys
     594  in
     595      if e < lb
     596      then Right (1.414*h)              (* step too small, accept but grow *)
     597      else (if e < ub
     598            then Right h                (* step just right *)
     599            else Left (0.5*h))          (* step too large, reject and shrink *)
     600  end
     601
     602
     603fun integrate (f,x,ys,h,i) =
     604  let open Real
     605      val stepper = make_stepper f
     606      fun recur (hs) =
     607          let
     608              val (stn, etn) = (stepper hs) (x,ys)
     609          in
     610              if hs < 1.0 andalso Real.== (hs, hf)
     611              then (hs, xn, stn)
     612              else
     613                  (case predictor tol (h,etn) of
     614                       Left bad => recur bad
     615                     | Right good => good::stn)
     616          end
     617  in
     618      recur h
     619  end
     620
     621EOF
     622
     623               ))
     624            (else
     625             `(
     626               ("fun integrate (f,x: real,y: real list,h,i) = ((make_stepper f) h) (x,y)" ,nl)
     627               ))
     628            ))
     629       )
    580630))
    581631
Note: See TracChangeset for help on using the changeset viewer.