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

[Mirko Vukovic]

Liam,

That is good news.  You caught me as I was debugging as well.

I was originally doing the problem on clisp+cygwin.  I am going to continue
on sbcll+linux.  In addition, I rewrote the problem to use the vanderpol
equations that are already present in the ode-example.lisp file.

speculation: On clisp+cygwin, it looks as if marray values are not being
passed in and out of apply-step.  For example, if I pre-load dydt-out with
some arbitrary values before calling apply-step, it will contain those same
values after the call.  I would have expected it to contain the values at
the new time step.

On clisp+cygwin I did not install the fsbv library: I understood that it was
not essential for gsll

Again, I will redo this on sbcl, and compare traces.

On Thu, Dec 3, 2009 at 10:15 PM, Liam Healy <lhealy at common-lisp.net> wrote:

> 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))
> >
> > |#
> >
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.common-lisp.net/pipermail/gsll-devel/attachments/20091203/af2c2e2b/attachment.html>