[ create a new paste ] login | about

Link: http://codepad.org/gSEoxetW    [ raw code | output | fork ]

Scheme, pasted on Apr 22:
; the Runge–Kutta methods are for the approximation of solutions of 
; ordinary differential equations.

(require (lib "6.ss" "srfi"))
(require (lib "40.ss" "srfi"))
(require (lib "43.ss" "srfi"))

(define integrate-system
  (lambda (system-derivative initial-state h)
    (let ((next (runge-kutta-4 system-derivative h)))
      (letrec ((states
                (cons initial-state
                      (delay (map-streams next
                                          states)))))
        states))))

(define runge-kutta-4
  (lambda (f h)
    (let ((*h (scale-vector h))
          (*2 (scale-vector 2))
          (*1/2 (scale-vector (/ 1 2)))
          (*1/6 (scale-vector (/ 1 6))))
      (lambda (y)
        ;; y is a system state
        (let* ((k0 (*h (f y)))
               (k1 (*h (f (add-vectors y (*1/2 k0)))))
               (k2 (*h (f (add-vectors y (*1/2 k1)))))
               (k3 (*h (f (add-vectors y k2)))))
          (add-vectors y
            (*1/6 (add-vectors k0
                               (*2 k1)
                               (*2 k2)
                               k3))))))))

(define elementwise
  (lambda (f)
    (lambda vectors
      (generate-vector
        (vector-length (car vectors))
        (lambda (i)
          (apply f
                 (map (lambda (v) (vector-ref  v i))
                      vectors)))))))

(define generate-vector
  (lambda (size proc)
    (let ((ans (make-vector size)))
      (letrec ((loop
                (lambda (i)
                  (cond ((= i size) ans)
                        (else
                         (vector-set! ans i (proc i))
                         (loop (+ i 1)))))))
        (loop 0)))))

(define add-vectors (elementwise +))

(define scale-vector
  (lambda (s)
    (elementwise (lambda (x) (* x s)))))

(define map-streams
  (lambda (f s)
    (cons (f (head s))
          (delay (map-streams f (tail s))))))

(define head car)
(define tail
  (lambda (stream) (force (cdr stream))))

(define damped-oscillator
  (lambda (R L C)
    (lambda (state)
      (let ((Vc (vector-ref state 0))
            (Il (vector-ref state 1)))
        (vector (- 0 (+ (/ Vc (* R C)) (/ Il C)))
                (/ Vc L))))))

(define the-states
  (integrate-system
     (damped-oscillator 10000 1000 .001)
     '#(1 0)
     .01))

(letrec ( 
   (n 100) 
   (loop (lambda (s) 
      (if (not (= n 0)) 
          (begin (set! n (- n 1)) 
                 (write (head s)) 
                 (newline)
                 (loop (tail s))))))) 
  (loop the-states)) 


Output:
#2(1 0)
#2(0.998950533570875 9.994835082916667e-06)
#2(0.9978022717932012 1.997868135089848e-05)
#2(0.9965554281807733 2.9950551909982803e-05)
#2(0.9952102258871526 3.9909462049570005e-05)
#2(0.9937668976737287 4.985442933866221e-05)
#2(0.9922256858768516 5.978447372177803e-05)
#2(0.9905868423740402 6.969861761453393e-05)
#2(0.9888506285492711 7.959588599888321e-05)
#2(0.987017315257352 8.947530651800312e-05)
#2(0.9850871827873835 9.93359095708213e-05)
#2(0.9830605208253149 0.0001091767284061722)
#2(0.9809376284155985 0.00011899679921657462)
#2(0.9787188139219463 0.00012879516123162125)
#2(0.9764043949871944 0.00013857085681097134)
#2(0.9739946984922812 0.0001483229315369376)
#2(0.9714900605143414 0.0001580504343066583)
#2(0.9688908262839245 0.00016775241742384626)
#2(0.9661973501413404 0.00017742793669010525)
#2(0.9634099954921378 0.00018707605149580582)
#2(0.9605291347617215 0.0001966958249105115)
#2(0.9575551493491138 0.0002062863237729468)
#2(0.9544884295798644 0.0002158466187804987)
#2(0.9513293746581157 0.0002253757845782429)
#2(0.9480783926178291 0.0002348728998474867)
#2(0.9447359002731759 0.0002443370473938197)
#2(0.941302323168102 0.00025376731423466434)
#2(0.9377780955250695 0.0002631627916863184)
#2(0.9341636601929825 0.00027252257545048)
#2(0.9304594685943022 0.0002818457657002486)
#2(0.9266659806713586 0.00029113146716559275)
#2(0.9227836648318641 0.0003003787892182768)
#2(0.9188129978936352 0.00030958684595623924)
#2(0.9147544650285292 0.0003187547562874139)
#2(0.9106085597056015 0.00032788164401298696)
#2(0.9063757836334912 0.0003369666379100813)
#2(0.9020566467020399 0.0003460088718138613)
#2(0.8976516669231515 0.0003550074846990495)
#2(0.8931613703708995 0.0003639616207608487)
#2(0.8885862911208882 0.00037287042949526083)
#2(0.8839269711888746 0.00038173306577879613)
#2(0.8791839604686585 0.0003905486899475647)
#2(0.8743578166692465 0.00039931646787574323)
#2(0.8694491052512983 0.00040803557105341013)
#2(0.8644583993628617 0.00041670517666374116)
#2(0.8593862797744035 0.00042532446765955906)
#2(0.854233334813143 0.00043389263283923)
#2(0.8490001602966964 0.0004424088669218999)
#2(0.8436873594660389 0.0004508723706220639)
#2(0.838295542917791 0.00045928235072346165)
#2(0.8328253285358383 0.0004676380201522927)
#2(0.8272773414222911 0.00047593859804974435)
#2(0.821652213827791 0.00048418330984382606)
#2(0.8159505850811727 0.0004923713873205035)
#2(0.8101731015184898 0.0005005020686941261)
#2(0.8043204164114096 0.0005085745986771421)
#2(0.7983931898949874 0.0005165882285490934)
#2(0.7923920888948264 0.0005245422162248865)
#2(0.7863177870536316 0.0005324358263223308)
#2(0.7801709646571668 0.0005402683302289399)
#2(0.77395230855962 0.0005480390061679897)
#2(0.7676625121083883 0.0005557471392638268)
#2(0.7613022750682882 0.0005633920216064225)
#2(0.7548723035452003 0.0005709729523151652)
#2(0.7483733099091572 0.0005784892376018873)
#2(0.7418060127168813 0.0005859401908331199)
#2(0.7351711366337821 0.0005933251325915696)
#2(0.7284694123554208 0.0006006433907368139)
#2(0.7217015765284511 0.0006078943004652073)
#2(0.7148683716710432 0.0006150772043689949)
#2(0.7079705460928016 0.0006221914524946278)
#2(0.7010088538141838 0.0006292364024002746)
#2(0.6939840544854279 0.0006362114192125252)
#2(0.6868969133049998 0.0006431158756822811)
#2(0.6797482009375667 0.0006499491522398282)
#2(0.6725386934315061 0.0006567106370490864)
#2(0.6652691721359594 0.0006633997260610324)
#2(0.6579404236174388 0.0006700158230662913)
#2(0.6505532395759958 0.0006765583397468905)
#2(0.6431084167609602 0.0006830266957271751)
#2(0.6356067568862592 0.0006894203186238769)
#2(0.6280490665453246 0.0006957386440953359)
#2(0.620436157125597 0.0007019811158898676)
#2(0.6127688447226373 0.0007081471858932745)
#2(0.6050479500538524 0.0007142363141754958)
#2(0.5972742983718456 0.0007202479690363932)
#2(0.5894487193774007 0.0007261816270506678)
#2(0.5815720471321079 0.0007320367731119058)
#2(0.5736451199706413 0.0007378129004757479)
#2(0.5656687804126972 0.0007435095108021801)
#2(0.5576438750746021 0.0007491261141969424)
#2(0.5495712545805995 0.0007546622292520515)
#2(0.5414517734738246 0.0007601173830854356)
#2(0.5332862901269765 0.0007654911113796764)
#2(0.5250756666526971 0.0007707829584198573)
#2(0.5168207688136656 0.0007759924771305131)
#2(0.5085224659324176 0.0007811192291116803)
#2(0.5001816308008994 0.0007861627846740428)
#2(0.4917991395897653 0.0007911227228731737)
#2(0.4833758717574279 0.0007959986315428668)


Create a new paste based on this one


Comments: