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

[Mirko Vukovic]

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/20091201/8612f0b0/attachment.html>