Changeset 40140 in project
 Timestamp:
 05/25/21 16:11:01 (3 weeks ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

wiki/eggref/5/csm
r40139 r40140 346 346 (export integratesystem head tail) 347 347 (include "integrateimplementation.scm")) 348 % cat all.options 348 % cat integrateimplementation.scm 349 ; The procedure integratesystem integrates the system 350 ; yâ²k = fk(y1, y2, . . . , yn), k = 1, . . . , n 351 ; of differential equations with the method of RungeKutta. 352 ; The parameter systemderivative is a function that 353 ; takes a system state (a vector of values for the state vari 354 ; ables y1, . . . , yn) and produces a system derivative (the val 355 ; ues yâ²1, . . . , yâ²n). The parameter initialstate provides 356 ; an initial system state, and h is an initial guess for the 357 ; length of the integration step. 358 ; The value returned by integratesystem is an infi 359 ; nite stream of system states. 360 361 (define (integratesystem systemderivative initialstate h) 362 (let ((next (rungekutta4 systemderivative h))) 363 (letrec ((states 364 (cons initialstate 365 (delay (mapstreams next states))))) 366 states))) 367 368 (define (rungekutta4 f h) 369 (let ((*h (scalevector h)) 370 (*2 (scalevector 2)) 371 (*1/2 (scalevector (/ 1 2))) 372 (*1/6 (scalevector (/ 1 6)))) 373 (lambda (y) 374 ;; y is a system state 375 (let* ((k0 (*h (f y))) 376 (k1 (*h (f (addvectors y (*1/2 k0))))) 377 (k2 (*h (f (addvectors y (*1/2 k1))))) 378 (k3 (*h (f (addvectors y k2))))) 379 (addvectors y 380 (*1/6 (addvectors k0 381 (*2 k1) 382 (*2 k2) 383 k3))))))) 384 385 (define (elementwise f) 386 (lambda vectors 387 (generatevector 388 (vectorlength (car vectors)) 389 (lambda (i) 390 (apply f 391 (map (lambda (v) (vectorref v i)) 392 vectors)))))) 393 394 (define (generatevector size proc) 395 (let ((ans (makevector size))) 396 (letrec ((loop 397 (lambda (i) 398 (cond ((= i size) ans) 399 (else 400 (vectorset! ans i (proc i)) 401 (loop (+ i 1))))))) 402 (loop 0)))) 403 404 (define addvectors (elementwise +)) 405 406 (define (scalevector s) 407 (elementwise (lambda (x) (* x s)))) 408 409 (define (mapstreams f s) 410 (cons (f (head s)) 411 (delay (mapstreams f (tail s))))) 412 413 (define head car) 414 415 (define (tail stream) 416 (force (cdr stream))) 417 ; % cat all.options 349 418 program example 350 419 </enscript>
Note: See TracChangeset
for help on using the changeset viewer.