Hello,
For some differential equations I tried, Jaime's rk routine was
somewhat slow.
So I wrote a simple implementation of the runge-kutta method for
autonomous systems using the 'reverse list trick' which Stavros
suggested before for the routine to produce successive positive terms
in a list.
I am including the code for Runge-Kutta and examples in the file RK_2.mac
below.
I would appreciate any suggestions for how to speed up the
computations.
Note that there are two routines for Runge-Kutta.
The first routine implements one step of integration of the ODE. The
second routine produces a list of length n so that if the step_size is
h, then the total time of integration is nh.
Thanks again for any suggestions
-sen
/********* begin file RK_2.mac ***********************/
/* 4th order Runge-Kutta for autonomous systems: x is an n-dim vector,
and VF(x) is an n-dim vector field */
rk_ode(VF,x,h):=block([k1,k2,k3,k4],
k1: h*VF(x), k2: h*VF(x+k1/2), k3: h*VF(x+k2/2), k4: h*VF(x+k3),
x: x + k1/6 + k2/3 + k3/3 + k4/6
)$
rk_usg: print("rk_ode(VF,x,h) gives one iteration of
the Euler method for solving autonomous ODE's. The vector field is
VF,
which has to be already defined as a vector function of the vector x.
")$
ode_examples_usg: print("
Example:
1. The Lorenz Differential Equations:
par: [sigma,rho,beta];
Lorenz(x):= [par[1]*(x[2]-x[1]), x[1]*(par[2] -x[3]),x[1]*x[2] -
par[3]*x[3]];
To get one iteration using the rk_ode method type
rk_ode(Lorenz,x,h)$
");
Lorenz_usg: print("The Lorenz differential equations have the form
x' = sigma*(y-x)
y' = x*(rho -z)
z' = x*y = beta*z
We use these as par: [sigma, rho, beta]
")$
par: [10.0, 28.0, 2.667];
Lorenz(x):= [par[1]*(x[2]-x[1]), x[1]*(par[2] -x[3]),x[1]*x[2] -
par[3]*x[3]];
rk_list(VF,init_point,step_size,NumPoints):=
block([list_1:[], list_2:[], y:init_point],
for i:1 thru NumPoints do
(y: rk_ode(VF,y,step_size), list_1: cons(y, list_1)),
list_2: reverse(list_1)
)$
rk_list_usg: print(" rk_list(VF,init_point,step_size,NumPoints)
gives a list of NumPoints iterates of VF using rk_ode for each
iterate. ");
print(" To use draw3d to plot the output
type
r_1: rk_list(Lorenz,[.1,0,0],.01,500)$
draw3d(points(r_1))$
");
/******* end file RK_2.mac *************/
--
---------------------------------------------------------------------------
| Sheldon E. Newhouse | e-mail: sen1 at math.msu.edu |
| Mathematics Department | |
| Michigan State University | telephone: 517-355-9684 |
| E. Lansing, MI 48824-1027 USA | FAX: 517-432-1562 |
---------------------------------------------------------------------------