Raymond Toy wrote:
> Sam Steingold wrote:
>> how do I solve an 1st order ODE IVP on maxima numerically?
>
> load(dynamics);
>
> There's an rk function for solving odes. I've never used it. "?? rk"
> gives some info. Looking in dynamics.mac gives an example:
>
>
> results: rk(t-x^2,x,1,[t,0,8,0.1])$
> plot2d([discrete, results]));
>
> Don't know if that's right or not.
runge-kutta sounds good but it has been running on my amd64 ws for 20+ minutes
now... (maxima/sbcl):
rk(dy/dx=cos(x*y),y,0,[x,0,5,0.1]);
in clisp
(defun ode1 (f &key (y0 0) (x0 0) (dx 0.1) (x1 1))
"Solve a first order ODE y'=f(x,y) for y(X0)=Y0 upto X1 with step DX.
Return the list of (x y) pairs."
(loop :for x :from x0 :to x1 :by dx
:for yprim = (funcall f x0 y0) :then (funcall f x y)
:for y = y0 :then (+ y (* dx yprim))
:collect (cons x y)))
did the job in a jiffy:
(defparameter *solutions6*
(mapcar (lambda (dx)
(cons (prin1-to-string dx)
(cllib:ode1 (lambda (x y) (cos (* x y))) :dx dx :x1 6)))
'(0.1 0.01 0.001 0.0001 0.00001)))
(cllib:plot-lists-arg *solutions6*
:title "y'=cos(xy)" :ylabel "y" :xlabel "x"
:legend '(:top :right :box))
thanks!