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))

|#