Changeset 5959 in project
 Timestamp:
 09/09/07 15:34:12 (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

odelmm/trunk/rkg4.scm
r5955 r5959 34 34 (declare (export ode:rkg4)) 35 35 36 (define Ct32 2)37 (define Ct21 (/ 3 2))38 36 39 37 (define Ch32 (/ 3 2)) … … 44 42 (define Ch33 (/ 9 8)) 45 43 44 (define C2t .33) 45 (define C3t .66) 46 46 47 47 48 ;; machine epsilon … … 51 52 (if (= 1.0 (+ 1.0 x1)) x (loop x1))))) 52 53 53 ;; unit roundoff error for IEEE doubleprecision floating point 54 (define u 4.4409e16) 55 56 ;; scaling tolerance when computing lambda 57 (define stol 1e6) 58 59 (define (rkg4 startit t tstep gdval) 60 61 (let loop ((it startit) (t t) (gdval gdval)) 62 63 (if gdval 64 (begin 65 ;; update the state values 54 55 (define (ode:rkg4 oderuntime indep tstart tstop initialtstep) 56 57 (define maxerror (oderuntime 'maxerror)) 58 (define hierror? (oderuntime 'hierror?)) 59 (define lowerror? (oderuntime 'lowerror?)) 60 61 (define printq (oderuntime 'print)) 62 (define evalrhs (oderuntime 'eval)) 63 (define isstate? (oderuntime 'isstate?)) 64 (define isasgn? (oderuntime 'isasgn?)) 65 (define solveenvref (oderuntime 'solveenvref)) 66 (define solveenvset! (oderuntime 'solveenvset!)) 67 (define evalenvref (oderuntime 'evalenvref)) 68 (define evalenvset! (oderuntime 'evalenvset!)) 69 70 (define svecderividx (oderuntime 'svecderividx)) 71 (define svecvalueidx (oderuntime 'svecvalueidx)) 72 (define svecrelerridx (oderuntime 'svecrelerridx)) 73 (define svecabserridx (oderuntime 'svecabserridx)) 74 75 (define (svecderiv svec) (f64vectorref svec svecderividx)) 76 (define (svecvalue svec) (f64vectorref svec svecvalueidx)) 77 (define (svecrelerr svec) (f64vectorref svec svecrelerridx)) 78 (define (svecabserr svec) (f64vectorref svec svecabserridx)) 79 80 (define (svecderivset! svec v) (f64vectorset! svec svecderividx v)) 81 (define (svecvalueset! svec v) (f64vectorset! svec svecvalueidx v)) 82 (define (svecrelerrset! svec v) (f64vectorset! svec svecrelerridx v)) 83 (define (svecabserrset! svec v) (f64vectorset! svec svecabserridx v)) 84 85 (define yorder (oderuntime 'order)) 86 (define (makeyvect) (makef64vector yorder 0.0)) 87 88 ;; maximum and minimum ratios by which to increase to decrease the 89 ;; step size 90 (define M (oderuntime' hmaxfactor)) 91 (define m (oderuntime' hminfactor)) 92 (define S (oderuntime' hscalefactor)) 93 (define relmax (oderuntime' relmax)) 94 (define absmax (oderuntime' absmax)) 95 96 (define hubfactor 0.1) 97 (define hlbfactor 100.0) 98 (define hbias 0.5) 99 100 (define yvect (makeyvect)) 101 (define y1vect (makeyvect)) 102 (define nvect (makevector yorder #f)) 103 104 ;; approximation of (y'/y_max)^4 105 (define lam4 #f) 106 107 ;; upper and lower bounds for timestep 108 (define hub #f) 109 (define hlb #f) 110 111 (define (evalderivs!) 112 ((oderuntime 'evalorderforeach) 113 (lambda (order sym rhs) 114 (let ((v (evalrhs rhs)) 115 (svec (solveenvref sym))) 116 (cond ((isstate? sym) (begin 117 (svecderivset! svec v))) 118 ((isasgn? sym) (begin 119 (evalenvset! sym v) 120 (svecvalueset! svec v)))))))) 121 122 (define (evalderivs t) 123 (let ((y1vect (makeyvect))) 124 (let ((oldt (evalenvref indep))) 125 (evalenvset! indep t) 126 ((oderuntime 'evalorderfold) 127 (lambda (order sym rhs i) 128 (let ((v (evalrhs rhs)) 129 (svec (solveenvref sym))) 130 (cond ((isstate? sym) (f64vectorset! y1vect i v))) 131 (fx+ 1 i))) 0) 132 (evalenvset! indep oldt) 133 y1vect))) 134 135 136 137 ;; 138 ;; Based on CVYddNorm from CVODE 139 ;; 140 ;; This routine computes an estimate of the second derivative of y 141 ;; using a difference quotient, and returns its WRMS norm. 142 ;; 143 (define (yddnorm t hg) 144 ;; y < h*y'(t) + y(t) 145 (let ((y (blas:daxpy hg y1vect 1.0 yvect))) 146 ;; f(tn+hg, y, tempv, f_data); 147 ;; yd < f(t+hg, y(t))  y'(t) 148 (let* ((yd (evalderivs (+ t hg))) 149 (yd1 (blas:daxpy yorder 1.0 yd 1.0 y1vect)) 150 (ydd (blas:dscal yorder (/ 1.0 hg) yd1))) 151 152 ;; Estimate y'' 153 (blas:dnrm2 yorder ydd)))) 154 155 156 ;; 157 ;; Based on CVUpperBoundH0 routine from CVODE 158 ;; 159 ;; This routine sets an upper bound on abs(h0) based on 160 ;; tdist = tstop  tstart and the values of y[i] / y'[i] 161 ;; 162 (define (computeh0upperbound tdist) 163 (define vabsmax? (f64vector? absmax)) 164 165 (let ((yabs (f64vectorabs yvect)) 166 (y1abs (f64vectorabs y1vect))) 167 (let ((yabs+tol (if vabsmax? 168 (blas:daxpy yorder hubfactor yabs 1.0 absmax) 169 (blas:daxpy yorder hubfactor yabs absmax (f64vector yorder 1.0))))) 170 171 (let ((y1/y (f64vectordiv y1abs yabs+tol))) 172 (let ((y1/ymaxnorm (f64vectorref y1/y (blas:damax yorder y1/y))) 173 (hub (* hubfactor tdist))) 174 (if (> (* hub hubinv) 1.0) 175 (values (/ 1.0 y1/ymaxnorm) y1/ymaxnorm) 176 (values hub y1/ymaxnorm))))))) 177 178 179 ;; 180 ;; Based on CVHin routine from CVODE 181 ;; 182 ;; This routine computes a tentative initial step size h0. 183 ;; 184 ;; The algorithm used seeks to find h0 as a solution of 185 ;; 186 ;; (WRMS norm of (h0^2 ydd / 2)) = 1, 187 ;; 188 ;; where ydd = estimated second derivative of y. 189 ;; 190 (define (initialh) 191 (define maxiters 4) 192 193 (let ((tdiff ( tstop tstart))) 194 (let ((sign (negative? tdiff)) 195 (tdist (abs tdiff)) 196 (tround (let ((tround (* uround (max (abs tstart) (abs tstop))))) 197 (if (< tdist (* 2.0 tround)) (error "tstop too close to tstart")) 198 tround))) 199 ;; Set lower and upper bounds on h0, and take geometric mean 200 ;; Exit with this value if the bounds cross each other 201 (let ((hlb (* hlbfactor tround))) 202 (letvalues (((hub y1/ymaxnorm) (computeh0upperbound tdist))) 203 (set! lam4 (expt y1/ymaxnorm 4)) 204 (set! hlb hlb) 205 (set! hub hub) 206 (let ((hg (sqrt (* hlb ulb)))) 207 (if (< hub hlb) 208 (if sign ( hg) hg) 209 ;; Loop up to maxiters times to find h0. 210 ;; Stop if new and previous values differ by a factor < 2. 211 ;; Stop if hnew/hg > 2 after one iteration, as this probably means 212 ;; that the ydd value is bad because of cancellation error. 213 (let ((hnew (let loop ((i 0) (hnew hg) (hprev hg) (hrat 0.0)) 214 (cond ((>= i maxiters) hnew) 215 ((and (> hrat 0.5) (< hrat 2.0)) hnew) 216 ((and (>= i 2) (> hrat 2.0)) hprev) 217 (else 218 (let* ((hgs (* hnew (if sign 1 1))) 219 (yddnrm (yddnorm tstart hgs)) 220 (hnew1 (if (> (* yddnrm hub hub) 2.0) 221 (sqrt (/ 2.0 yddnrm)) (sqrt (* hnew hub))))) 222 (loop (fx+ 1 i) hnew1 hnew (/ hnew1 hnew)))))))) 223 ;; Apply bounds, bias factor, and attach sign. 224 (let ((h0 (* hbias hnew))) 225 (cond ((< h0 hlb) hlb) 226 ((> h0 hub) hub) 227 (else (if sign ( h0) h0)))))))))))) 228 229 230 231 (define (rkg4 startit t tstep gdval h0) 232 233 (let loop ((it startit) (t t) (gdval gdval)) 234 235 (if gdval 236 (begin 237 ;; update the state values 66 238 ((oderuntime 'solveenvstatefold) 67 239 (lambda (sym svec i) … … 71 243 (fx+ i 1)) 0) 72 244 ;; update lambda (eq. 5.16) 73 (set! lam (/ (blas:dnrm2 order y1vect) (blas:dnrm2 order yvect)))245 (set! lam (/ y1norm ynorm)) 74 246 ;; update independent variable 75 247 (solveenvset! indep svecvalueidx t) … … 82 254 83 255 (if (t<tstop tstep t tstop) 84 ( begin256 (let ((h (/ tstep 3))) 85 257 86 258 (if (> (* tstep ( (+ t tstep) tstop)) 0.0) … … 91 263 ((oderuntime 'solveenvstatefold) 92 264 (lambda (sym svec i) 93 (let* ((k1 (* tstep(svecderiv svec)))265 (let* ((k1 (* h (svecderiv svec))) 94 266 (y31 (+ (svecvalue svec) k1))) 95 267 (f64vectorset! k1vect i k1) … … 99 271 100 272 101 (evalenvset! indep (+ t tstep))273 (evalenvset! indep (+ t h)) 102 274 (evalderivs) 103 275 … … 107 279 (lambda (sym svec i) 108 280 (let* ((k1 (f64vectorref k1vect i)) 109 (k2 (* tstep(svecderiv svec)))281 (k2 (* h (svecderiv svec))) 110 282 (y32 (+ (svecvalue svec) k1 k2))) 111 283 (f64vectorset! k2vect i k2) … … 114 286 (cons (fx+ 1 i) (cdr hstlist)))) 0) 115 287 116 (evalenvset! indep (+ t (* Ct32 tstep))) 288 289 (evalenvset! indep (+ t (* h 2.0))) 117 290 (evalderivs) 118 291 … … 123 296 (let* ((k1 (f64vectorref k1vect i)) 124 297 (k2 (f64vectorref k2vect i)) 125 (k3 (* tstep(svecderiv svec)))298 (k3 (* h (svecderiv svec))) 126 299 (y33 (+ (svecvalue svec) k1 k2 k3))) 127 300 (f64vectorset! k3vect i k3) 128 301 (f64vectorset! y33vect i y33) 129 302 (cons (fx+ 1 i) (cdr hstlist)))) 0) 303 304 305 (evalenvset! indep t) 306 (evalderivs) 130 307 131 308 … … 139 316 (cons (fx+ 1 i) (cdr hstlist)))) 0) 140 317 141 (evalenvset! indep (+ t (* Ct21 tstep))) 318 319 (evalenvset! indep (+ t (* h 1.5))) 142 320 (evalderivs) 143 321 … … 147 325 (lambda (sym svec i) 148 326 (let* ((k1 (f64vectorref k1vect i)) 149 (k4 (* tstep(svecderiv svec)))327 (k4 (* h (svecderiv svec))) 150 328 (y22 (+ (svecvalue svec) (* 1.5 k1) (* 1.5 k4)))) 151 329 (f64vectorset! k4vect i k4) … … 153 331 (cons (fx+ 1 i) (cdr hstlist)))) 0) 154 332 333 (evalenvset! indep t) 334 (evalderivs) 155 335 156 336 ;; step y11 … … 187 367 (vectorset! nvect i (f64vector y0 hy1 h2y2 h3y3)))))) 0) 188 368 189 190 ;; Advance state to next step 191 ((oderuntime 'solveenvstatefold) 192 (lambda (sym svec i) 193 (let* ((y33 (f64vectorref y33vect i))) 194 (evalenvset! sym y33) 195 (cons (fx+ 1 i)))) 0) 196 197 198 (evalderivs) 199 369 ;; compute new estimate of lambda 370 371 ;; define relative error as (max ( hub (* h^4 lam4)) ( (* h^4 lam4) hlb)) 200 372 201 373 (let ((err (maxerror))) … … 217 389 (evalenvset! sym (svecvalue svec)))) 218 390 (loop it t #f))) 219 (else (loop (fx+ 1 it) (+ t tstep) #t))))))) 391 (else (loop (fx+ 1 it) (+ t tstep) #t)))))))) 220 392 221 (let* ((ynorm ...) 222 (y4norm (* (expt lam 4) ynorm)) 223 (hm (* S (expt (/ eps (* Cm y4norm)) 0.25)))) 224 (values it t tstep (list hm nvect))))) 393 (values it t tstep (list nvect)))
Note: See TracChangeset
for help on using the changeset viewer.