Hi,
I am trying to solve a second order ODE by calling a stepper for a number of fixed-size steps. I am getting an error about conversion to foreign type c-pointer.
I modeled my call to `apply-step' based on `apply-evolution'. Looking at the source files, that is likely not to work because `apply-evolution' uses c-pointers, and `apply-step' makes direct variable references. But I don't know how to deal with those - my C and CFFI knowledge are non-existent.
Anyway, here is a sample code that I hope someone can improve:
;;;; 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.1d0) (stepper-type +step-rk2+)) (let ((stepper (make-ode-stepper stepper-type 2 #'sin-ode #'sin-ode-jacobian)) (dep (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 dep 0) 0d0 (maref dep 1) 1d0) (dotimes (i 10) (incf time delta-t) (apply-step stepper time dep delta-t yerr dydt-in dydt-out))))
|#
BTW, once this is working, I will submit it as a patch to the ode-examples.lisp file
Thanks,
Mirko
OK fixed, thanks for the report. You might want to return a value from your calculation, viz.
(let ((time 0d0) (delta-t 0.1d0) (stepper-type +step-rk2+)) (let ((stepper (make-ode-stepper stepper-type 2 #'sin-ode #'sin-ode-jacobian)) (dep (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 dep 0) 0d0 (maref dep 1) 1d0) (dotimes (i 10) (incf time delta-t) (apply-step stepper time dep delta-t yerr dydt-in dydt-out)) dep))
Liam
On Mon, Nov 30, 2009 at 8:30 AM, Mirko Vukovic mirko.vukovic@gmail.com wrote:
Hi,
I am trying to solve a second order ODE by calling a stepper for a number of fixed-size steps. I am getting an error about conversion to foreign type c-pointer.
I modeled my call to `apply-step' based on `apply-evolution'. Looking at the source files, that is likely not to work because `apply-evolution' uses c-pointers, and `apply-step' makes direct variable references. But I don't know how to deal with those - my C and CFFI knowledge are non-existent.
Anyway, here is a sample code that I hope someone can improve:
;;;; 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.1d0) (stepper-type +step-rk2+)) (let ((stepper (make-ode-stepper stepper-type 2 #'sin-ode #'sin-ode-jacobian)) (dep (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 dep 0) 0d0 (maref dep 1) 1d0) (dotimes (i 10) (incf time delta-t) (apply-step stepper time dep delta-t yerr dydt-in dydt-out))))
|#
BTW, once this is working, I will submit it as a patch to the ode-examples.lisp file
Thanks,
Mirko
Gsll-devel mailing list Gsll-devel@common-lisp.net http://common-lisp.net/cgi-bin/mailman/listinfo/gsll-devel
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))
|#
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@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))
|#
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@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@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
(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))
|#
-
Liam Healy -
Mirko Vukovic