[Gsll-devel] ODE apply-step: argument passing problem

[Liam Healy]

Mirko,

Here is what I get when I trace sine-ode
  0: (SIN-ODE 0.0d0 0.0d0 1.0d0)
  0: SIN-ODE returned 1.0d0 -0.0d0
Y-in:       0.0100, 0.0d0, 1.0d0
dydt-in:      0.0100, 1.0d0, -0.0d0
  0: (SIN-ODE 0.015d0 0.005d0 1.0d0)
  0: SIN-ODE returned 1.0d0 -0.005d0
  0: (SIN-ODE 0.02d0 0.01d0 0.9999d0)
  0: SIN-ODE returned 0.9999d0 -0.01d0
  0: (SIN-ODE 0.02d0 0.009999833333333333d0 0.99995d0)
  0: SIN-ODE returned 0.99995d0 -0.009999833333333333d0
Y-out:       0.0100, 0.009999833333333333d0, 0.99995d0
dydt-out:      0.0100, 0.99995d0, -0.009999833333333333d0

Y-in:       0.0200, 0.009999833333333333d0, 0.99995d0
dydt-in:      0.0200, 0.99995d0, -0.009999833333333333d0
  0: (SIN-ODE 0.025d0 0.014999583333333334d0 0.9999000008333333d0)
  0: SIN-ODE returned 0.9999000008333333d0 -0.014999583333333334d0
  0: (SIN-ODE 0.03d0 0.01999833335d0 0.9997500066666667d0)
  0: SIN-ODE returned 0.9997500066666667d0 -0.01999833335d0
  0: (SIN-ODE 0.03d0 0.019998666683333333d0 0.9998000058333055d0)
  0: SIN-ODE returned 0.9998000058333055d0 -0.019998666683333333d0
Y-out:       0.0200, 0.019998666683333333d0, 0.9998000058333055d0
dydt-out:      0.0200, 0.9998000058333055d0, -0.019998666683333333d0

#<VECTOR-DOUBLE-FLOAT #(0.019998666683333333d0 0.9998000058333055d0)>

which superficially looks right to me.


On Tue, Dec 1, 2009 at 10:23 PM, Mirko Vukovic <mirko.vukovic at gmail.com> wrote:
> Hi Liam (& others)
>
> I am getting funny behavior with the ODE stepper:
>
> (As refresher, I am solving a 2nd order system that should reproduce a
> sinusoid)
>
> A trace of the dydt function (sin-ode below) shows that it is always called
> with arguments (time 0 0), even though the initial conditions are (time 0 1)
> (see the let form).
>
> It looks as if the initial conditions (stored in vector y below) are not
> passed as arguments to the dydt function (sin-ode).  Or am I missing
> something?
>
> I tried both the rk2 and rk4 steppers and different step sizes.
>
> Thanks,
>
> Mirko
>
>
> Here is the source:
>
> ;;;; Example of a simple time-stepper
> (defun sin-ode (time y z)
>   "Define ODE for a sinusoid
>
> y''=-y
>
> or as a system:
>
> y'=z
> z'=-y
>
> with initial conditions:
>
> y0=0
> z0=1 "
>   (declare (ignore time))
>   (values z (- y)))
>
> (defun sin-ode-jacobian (time y z)
>   (declare (ignore time y z))
>   (values 0d0  0d0
>       0d0  1d0
>       -1d0 0d0))
> #|
> (let  ((time 0d0)
>        (delta-t 0.01d0)
>        (stepper-type +step-rk2+))
>   (let ((stepper (make-ode-stepper stepper-type 2 #'sin-ode
> #'sin-ode-jacobian))
>     (y (make-marray 'double-float :dimensions 2))
>     (dydt-in (make-marray 'double-float :dimensions 2))
>     (dydt-out (make-marray 'double-float :dimensions 2))
>     (yerr (make-marray 'double-float :dimensions 2)))
>     (setf (maref y 0) 0d0 ;; y(0)=0
>       (maref y 1) 1d0);; y'(0)=1
>     (multiple-value-bind (yp zp) (sin-ode time (maref y 0) (maref y 1))
>       (setf (maref dydt-out 0) yp
>         (maref dydt-out 1) zp))
>     (dotimes (i 2)
>       (incf time delta-t)
>       (setf (maref dydt-in 0) (maref dydt-out 0)
>         (maref dydt-in 1) (maref dydt-out 1))
>       (format t "Y-in: ~12,4f, ~a, ~a~%" time (maref y 0) (maref y 1))
>       (format t "dydt-in:~12,4f, ~a, ~a~%" time (maref dydt-in 0) (maref
> dydt-in 1))
>       (apply-step stepper time y delta-t yerr dydt-in dydt-out)
>       (format t "Y-out: ~12,4f, ~a, ~a~%" time (maref y 0) (maref y 1))
>       (format t "dydt-out:~12,4f, ~a, ~a~%~%" time (maref dydt-out 0) (maref
> dydt-out 1)))
>     y))
>
> |#
>