rk (runge-kutta) of dynamics package really slow?

The rk package could be made faster in several ways.
Using "ev" is probably a bad idea; subst could be used and should be faster.
It is also possible to compile rk, but unless there are some declarations,
it may not help much.
rk could also be written in lisp directly.

RJF 

> -----Original Message-----
> From: maxima-bounces at math.utexas.edu 
> [mailto:maxima-bounces at math.utexas.edu] On Behalf Of Joachim De Beule
> Sent: Friday, May 09, 2008 4:55 AM
> To: maxima at math.utexas.edu
> Subject: rk (runge-kutta) of dynamics package really slow?
> 
> 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:
> 
> ==============================================
> load("dynamics")
> 
> 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
>