# [Maxima] rk (runge-kutta) of dynamics package really slow?

Sheldon Newhouse sen1 at math.msu.edu
Fri May 9 13:09:49 CDT 2008

```Joachim De Beule wrote:
> I tried to numerically solve a system of ODE's with maxima's rk function but
> it takes ages to terminate. When I solve the same set of ODE's in Octave, it
> almost immediately finishes. I prefer to use maxima however, so if anyone has
> a suggestion on how to increase performance I would be very grateful!
>
> This is what I did in maxima:
>
> ==============================================
>
> tau_m: 1.0
> s_m: 1.0
>
>
> eqm: s_m + f_1*mf_1 + f_2*mf_2 + f_3*mf_3 - m*(tau_m + (f_1 + f_2 + f_3)^2 +
> (f_1 - mf_1) + (f_2 - mf_2) + (f_3 - mf_3))\$
> eqf_1: -m*f_1 + (m + f_1)*mf_1\$
> eqf_2: -m*f_2 + (m + f_2)*mf_2\$
> eqf_3: -m*f_3 + (m + f_3)*mf_3\$
> eqmf_1: m*f_1 - (m + f_1)*mf_1 - mf_1*(f_2 + f_3) + f_1*( mf_2 + mf_3)\$
> eqmf_2: m*f_2 - (m + f_2)*mf_2 - mf_2*(f_1 + f_3) + f_2*( mf_1 + mf_3)\$
> eqmf_3: m*f_3 - (m + f_3)*mf_3 - mf_3*(f_2 + f_1) + f_3*( mf_2 + mf_1)\$
>
> msol: rk([eqm,eqf_1,eqf_2,eqf_3,eqmf_1,eqmf_2,eqmf_3],
> [m,f_1,f_2,f_3,mf_1,mf_2,mf_3],[0.7,0.5,0.3,0.2,0.15,0.15,0.2],
> [t,0,50,0.05])\$
> ==============================================
>
> and in octave:
>
> ==============================================
> global tau_m = 1.0;
> global s_m = 1.0;
>
> function xdot = f (x, t)
>
> global s_m;
> global tau_m;
>
> m = x(1);
> f_1 = x(2);
> f_2 = x(3);
> f_3 = x(4);
> mf_1 = x(5);
> mf_2 = x(6);
> mf_3 = x(7);
>
>
> xdot(1) = s_m + f_1*mf_1 + f_2*mf_2 + f_3*mf_3 - m*(tau_m + (f_1 + f_2 +
> f_3)^2 + (f_1 - mf_1) + (f_2 - mf_2) + (f_3 - mf_3));
> xdot(2) = -m*f_1 + (m + f_1)*mf_1;
> xdot(3) = -m*f_2 + (m + f_2)*mf_2;
> xdot(4) = -m*f_3 + (m + f_3)*mf_3;
> xdot(5) = m*f_1 - (m + f_1)*mf_1 - mf_1*(f_2 + f_3) + f_1*( mf_2 + mf_3);
> xdot(6) = m*f_2 - (m + f_2)*mf_2 - mf_2*(f_1 + f_3) + f_2*( mf_1 + mf_3);
> xdot(7) = m*f_3 - (m + f_3)*mf_3 - mf_3*(f_2 + f_1) + f_3*( mf_2 + mf_1);
>
> endfunction
>
> x0 = [0.7; 0.5; 0.3; 0.2; 0.15; 0.15; 0.2];
>
> t = linspace (0, 50, 1000)';
>
> x = lsode("f", x0, t);
> =============================================
>
> thanks, Joachim.
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>
>
I ran into this problem too, so I wrote my own  RK 4th order routine
(for autonomous equations.  It works much faster.  I am attaching the
code. Try it and let me know if it works.

-sen
-------------- next part --------------
A non-text attachment was scrubbed...
Name: RK_4.mac
Type: image/x-pict
Size: 3181 bytes
Desc: not available
Url : http://www.math.utexas.edu/pipermail/maxima/attachments/20080509/c62632c4/attachment.bin
```