Changeset 40140 in project


Ignore:
Timestamp:
05/25/21 16:11:01 (2 months ago)
Author:
felix winkelmann
Message:

csm: missing example code

File:
1 edited

Legend:

Unmodified
Added
Removed
  • wiki/eggref/5/csm

    r40139 r40140  
    346346  (export integrate-system head tail)
    347347  (include "integrate-implementation.scm"))
    348 % cat all.options
     348% cat integrate-implementation.scm
     349; The procedure integrate-system integrates the system
     350; y′k = fk(y1, y2, . . . , yn), k = 1, . . . , n
     351; of differential equations with the method of Runge-Kutta.
     352; The parameter system-derivative 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 initial-state 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 integrate-system is an infi-
     359; nite stream of system states.
     360
     361(define (integrate-system system-derivative initial-state h)
     362  (let ((next (runge-kutta-4 system-derivative h)))
     363    (letrec ((states
     364               (cons initial-state
     365                     (delay (map-streams next states)))))
     366      states)))
     367
     368(define (runge-kutta-4 f h)
     369  (let ((*h (scale-vector h))
     370        (*2 (scale-vector 2))
     371        (*1/2 (scale-vector (/ 1 2)))
     372        (*1/6 (scale-vector (/ 1 6))))
     373    (lambda (y)
     374      ;; y is a system state
     375      (let* ((k0 (*h (f y)))
     376             (k1 (*h (f (add-vectors y (*1/2 k0)))))
     377             (k2 (*h (f (add-vectors y (*1/2 k1)))))
     378             (k3 (*h (f (add-vectors y k2)))))
     379        (add-vectors y
     380          (*1/6 (add-vectors k0
     381                             (*2 k1)
     382                             (*2 k2)
     383                             k3)))))))
     384
     385(define (elementwise f)
     386  (lambda vectors
     387    (generate-vector
     388      (vector-length (car vectors))
     389      (lambda (i)
     390        (apply f
     391               (map (lambda (v) (vector-ref v i))
     392                    vectors))))))
     393
     394(define (generate-vector size proc)
     395  (let ((ans (make-vector size)))
     396    (letrec ((loop
     397               (lambda (i)
     398                 (cond ((= i size) ans)
     399                       (else
     400                         (vector-set! ans i (proc i))
     401                         (loop (+ i 1)))))))
     402      (loop 0))))
     403
     404(define add-vectors (elementwise +))
     405
     406(define (scale-vector s)
     407  (elementwise (lambda (x) (* x s))))
     408
     409(define (map-streams f s)
     410  (cons (f (head s))
     411        (delay (map-streams f (tail s)))))
     412
     413(define head car)
     414
     415(define (tail stream)
     416  (force (cdr stream)))
     417; % cat all.options
    349418-program example
    350419</enscript>
Note: See TracChangeset for help on using the changeset viewer.