Liam,<br><br>That is good news. You caught me as I was debugging as well. <br><br>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.<br>
<br>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.<br>
<br>On clisp+cygwin I did not install the fsbv library: I understood that it was not essential for gsll<br><br>Again, I will redo this on sbcl, and compare traces.<br><br><div class="gmail_quote">On Thu, Dec 3, 2009 at 10:15 PM, Liam Healy <span dir="ltr"><<a href="mailto:lhealy@common-lisp.net">lhealy@common-lisp.net</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Mirko,<br>
<br>
Here is what I get when I trace sine-ode<br>
0: (SIN-ODE 0.0d0 0.0d0 1.0d0)<br>
0: SIN-ODE returned 1.0d0 -0.0d0<br>
Y-in: 0.0100, 0.0d0, 1.0d0<br>
dydt-in: 0.0100, 1.0d0, -0.0d0<br>
0: (SIN-ODE 0.015d0 0.005d0 1.0d0)<br>
0: SIN-ODE returned 1.0d0 -0.005d0<br>
0: (SIN-ODE 0.02d0 0.01d0 0.9999d0)<br>
0: SIN-ODE returned 0.9999d0 -0.01d0<br>
0: (SIN-ODE 0.02d0 0.009999833333333333d0 0.99995d0)<br>
0: SIN-ODE returned 0.99995d0 -0.009999833333333333d0<br>
Y-out: 0.0100, 0.009999833333333333d0, 0.99995d0<br>
dydt-out: 0.0100, 0.99995d0, -0.009999833333333333d0<br>
<br>
Y-in: 0.0200, 0.009999833333333333d0, 0.99995d0<br>
dydt-in: 0.0200, 0.99995d0, -0.009999833333333333d0<br>
0: (SIN-ODE 0.025d0 0.014999583333333334d0 0.9999000008333333d0)<br>
0: SIN-ODE returned 0.9999000008333333d0 -0.014999583333333334d0<br>
0: (SIN-ODE 0.03d0 0.01999833335d0 0.9997500066666667d0)<br>
0: SIN-ODE returned 0.9997500066666667d0 -0.01999833335d0<br>
0: (SIN-ODE 0.03d0 0.019998666683333333d0 0.9998000058333055d0)<br>
0: SIN-ODE returned 0.9998000058333055d0 -0.019998666683333333d0<br>
Y-out: 0.0200, 0.019998666683333333d0, 0.9998000058333055d0<br>
dydt-out: 0.0200, 0.9998000058333055d0, -0.019998666683333333d0<br>
<br>
#<VECTOR-DOUBLE-FLOAT #(0.019998666683333333d0 0.9998000058333055d0)><br>
<br>
which superficially looks right to me.<br>
<div><div></div><div class="h5"><br>
<br>
On Tue, Dec 1, 2009 at 10:23 PM, Mirko Vukovic <<a href="mailto:mirko.vukovic@gmail.com">mirko.vukovic@gmail.com</a>> wrote:<br>
> Hi Liam (& others)<br>
><br>
> I am getting funny behavior with the ODE stepper:<br>
><br>
> (As refresher, I am solving a 2nd order system that should reproduce a<br>
> sinusoid)<br>
><br>
> A trace of the dydt function (sin-ode below) shows that it is always called<br>
> with arguments (time 0 0), even though the initial conditions are (time 0 1)<br>
> (see the let form).<br>
><br>
> It looks as if the initial conditions (stored in vector y below) are not<br>
> passed as arguments to the dydt function (sin-ode). Or am I missing<br>
> something?<br>
><br>
> I tried both the rk2 and rk4 steppers and different step sizes.<br>
><br>
> Thanks,<br>
><br>
> Mirko<br>
><br>
><br>
> Here is the source:<br>
><br>
> ;;;; Example of a simple time-stepper<br>
> (defun sin-ode (time y z)<br>
> "Define ODE for a sinusoid<br>
><br>
> y''=-y<br>
><br>
> or as a system:<br>
><br>
> y'=z<br>
> z'=-y<br>
><br>
> with initial conditions:<br>
><br>
> y0=0<br>
> z0=1 "<br>
> (declare (ignore time))<br>
> (values z (- y)))<br>
><br>
> (defun sin-ode-jacobian (time y z)<br>
> (declare (ignore time y z))<br>
> (values 0d0 0d0<br>
> 0d0 1d0<br>
> -1d0 0d0))<br>
> #|<br>
> (let ((time 0d0)<br>
> (delta-t 0.01d0)<br>
> (stepper-type +step-rk2+))<br>
> (let ((stepper (make-ode-stepper stepper-type 2 #'sin-ode<br>
> #'sin-ode-jacobian))<br>
> (y (make-marray 'double-float :dimensions 2))<br>
> (dydt-in (make-marray 'double-float :dimensions 2))<br>
> (dydt-out (make-marray 'double-float :dimensions 2))<br>
> (yerr (make-marray 'double-float :dimensions 2)))<br>
> (setf (maref y 0) 0d0 ;; y(0)=0<br>
> (maref y 1) 1d0);; y'(0)=1<br>
> (multiple-value-bind (yp zp) (sin-ode time (maref y 0) (maref y 1))<br>
> (setf (maref dydt-out 0) yp<br>
> (maref dydt-out 1) zp))<br>
> (dotimes (i 2)<br>
> (incf time delta-t)<br>
> (setf (maref dydt-in 0) (maref dydt-out 0)<br>
> (maref dydt-in 1) (maref dydt-out 1))<br>
> (format t "Y-in: ~12,4f, ~a, ~a~%" time (maref y 0) (maref y 1))<br>
> (format t "dydt-in:~12,4f, ~a, ~a~%" time (maref dydt-in 0) (maref<br>
> dydt-in 1))<br>
> (apply-step stepper time y delta-t yerr dydt-in dydt-out)<br>
> (format t "Y-out: ~12,4f, ~a, ~a~%" time (maref y 0) (maref y 1))<br>
> (format t "dydt-out:~12,4f, ~a, ~a~%~%" time (maref dydt-out 0) (maref<br>
> dydt-out 1)))<br>
> y))<br>
><br>
> |#<br>
><br>
</div></div></blockquote></div><br>