or is there an easier way to do it!?
Hey, I am currently trying to program a system of ODE's in a 4th order Runge-Kutta Method in maple. Maple doesn't seem to like my code! Whats wrong
restart;
with(plots):
f:=(x,y,t)->x(t)+2*y-t;
x0:=-1.0;
g:=(x,y,t)->2*y(t)-x-5*t;
y0:=2.0;
t0:=0;
xf:=0.5;
n:=2;
The following procedure estimates the solution of ordinary differential equations at a point xf.
n = number of steps
x0 = boundary condition for x
y0 = boundary condition for y
xf = value at which solution is desired
f = differential equation in the form
d
--- x(t) = f(x, y,t)
dt
g = differential equation in the form
d
--- y(t) = g(x, y,t)
dt
RK4th:=proc(n,x0,y0,xf,f,t0,g)
local X,Y,t,h,i,k1,k2,k3,k4,g1,g2,g3,g4:
h:=(xf-x0)/n:
X:=y0:
Y:=x0:
t:=t0:
for i from 0 by 1 to n-1 do
k1:=f(X,y0+i*h,t);
g1:=g(x0+i*h,Y,t);
k2:=f(X+0.5*k1*h,y0+i*h+0.5*g1*h,t+0.5*h);
g2:=g(x0+i*h+0.5*g1*h,y0+0.5*k1*h,t+0.5*h);
k3:=f(X+0.5*k2*h,y0+i*h+0.5*g2*h,t+0.5*h);
g3:=g(x0+i*h+0.5*g2*h,Y+0.5*k2*h,t+0.5*h);
k4:=f(X+k3*h,y0+i*h+g3*h,t+h);
g4:=g(x0+i*h+g3*h,Y+k3*h,t+h);
Y:=Y+1/6*(g1+2*g2+2*g3+g4)*h;
X:=X+1/6*(k1+2*k2+2*k3+k4)*h;
end do:
return (Y):
return (X):
end proc:
ODE1:=diff(y(t),t)=g(x,y,t);
ODE2:=diff(x(t),t)=f(x,y,t);
soln1:=(dsolve({ODE1,y(t0)=y0}));
assign(soln1):
y(t);
soln2:=(dsolve({ODE2,x(t0)=x0}));
assign(soln2):
x(t);
EV1:=evalf(subs(t=xf,y(t)));
EV2:=evalf(subs(t=xf,x(t)));
X[0]:=x0;
Y[0]:=y0;
h:=(xf-x0)/n;
for i from 1 by 1 to n do
X[i]:=x0+i*h:
Y[i]:=y0+i+h:
t[i]:=t0+i+h:
k1:=f(X[i-1],Y[i-1],t[i-1]);
g1:=g(X[i-1],Y[i-1],t[i-1]);
k2:=f(X[i-1]+0.5*k1*h,Y[i-1]+0.5*g1*h,t[i-1]+0.5*h);
g2:=g(X[i-1]+0.5*g1*h,Y[i-1]+0.5*k1*h,t[i-1]+0.5*h);
k3:=f(X[i-1]+0.5*k2*h,Y[i-1]+0.5*g2*h,t[i-1]+0.5*h);
g3:=g(X[i-1]+0.5*g2*h,Y[i-1]+0.5*k2*h,t[i-1]+0.5*h);
k4:=f(X[i-1]+k3*h,Y[i-1]+g3*h,t[i-1]+h);
g4:=g(X[i-1]+g3*h,Y[i-1]+k3*h,t[i-1]+h);
Y[i]:=Y[i-1]+1/6*(g1+2*g2+2*g3+g4)*h;
X[i]:=X[i-1]+1/6*(k1+2*k2+2*k3+k4)*h;
end do:
data:=[seq([x[i],Y[i]],i=0..n)];
Not sure if this helps, but have you tried to use Sage?
desolve_system_rk4 uses a wrapper around maxima, and in the documentation, the example they give to solve a Lotka Volterra systems seems syntactically relatively easy.