Numerical solution of ordinary differential equtions
Cauchy's initial value problem
Mirko Navara (navara@math.feld.cvut.cz) & Ale\305\241 N\304\233me\304\215ek (nemecek@math.feld.cvut.cz)
Program inicialization
res := proc (xk, k, metoda1, metoda2) local j,metoda3,metoda4; global x,y,h;
if nargs>5 then metoda3:=args[5]; metoda4:=args[6] fi;
h:=(xk-x0)/k;
x:=array(0..k); for j from 0 to k do x[j]:=x0+j*h od;
y:=array(0..k); y[0]:=y0;
userinfo(4, res, `x =`, x0, `y =`, y0);
for j from 1 to k
do y[j]:=metoda1(x, y, j-1, h, metoda2, metoda3, metoda4);
userinfo(4, res, `x =`, x[j], `y =`, y[j])
od;
y[k]
end: ##
Runge-Kutta methods (single-step methods):
Common way call:
RK := proc (x,y,j,h,metoda2) metoda2(x[j],y[j],h) end: ##
Procedures for separete single-step methods:
## Euler's Method
Euler := proc (x,y,h) local k1;
k1:=f(x,y); userinfo(8, res, `k1 =`, k1);
y+h*k1 end: ## Euler's Method
## Euler's Method - 1st modification
EulerMod1 := proc (x,y,h) local k1,k2,y2;
k1:=f(x,y); y2:=y+h/2*k1; userinfo(16, res, `y2 =`, y2);
k2:=f(x+h/2,y2); userinfo(8, res, `k1 =`, k1, `k2 =`, k2);
y+h*k2 end: ## Euler's Method - 1st modification
## Euler's Method - 2nd modification (Heune's method)
EulerMod2 := proc (x,y,h) local k1,k2,kk,y2;
k1:=f(x,y); y2:=y+h*k1; userinfo(16, res, `y2 =`, y2);
k2:=f(x+h,y2); userinfo(8, res, `k1 =`, k1, `k2 =`, k2);
kk:=(k1+k2)/2; userinfo(12, res, `kk =`, kk);
y+h*kk end: ## Euler's Method - 2nd modification (Heune's method)
## Runge-Kutta method 3rd order
RungeKutta3 := proc (x,y,h) local k1,k2,k3,kk,y2,y3;
k1:=f(x,y); y2:=y+h/2*k1;
k2:=f(x+h/2, y2); y3:=y+h*(2*k2-k1);
k3:=f(x+h, y3); userinfo(16, res, `y2 =`, y2, `y3 =`, y3);
userinfo(8, res, `k1 =`, k1, `k2 =`, k2, `k3 =`, k3);
kk:=(k1+4*k2+k3)/6; userinfo(12, res, `kk =`, kk);
y+h*kk end: ## Runge-Kutta method 3rd order
## Runge-Kutta method 4th order
RungeKutta4 := proc (x,y,h) local k1,k2,k3,k4,kk,y2,y3,y4;
k1:=f(x,y); y2:=y+h/2*k1;
k2:=f(x+h/2, y2); y3:=y+h/2*k2;
k3:=f(x+h/2, y3); y4:=y+h*k3;
k4:=f(x+h, y4);
userinfo(16, res, `y2 =`, y2, `y3 =`, y3, `y4 =`, y4);
userinfo(8, res, `k1 =`, k1, `k2 =`, k2, `k3 =`, k3, `k4 =`, k4);
kk:=(k1+2*k2+2*k3+k4)/6; userinfo(12, res, `kk =`, kk);
y+h*kk end: ## Runge-Kutta method 4th order
Explicit methods (metoda2 is a start single-step method):
## Adams-Bashforth method 2nd order
AB2 := proc (x,y,j,h,metoda2)
if j<1 then metoda2(x[j],y[j],h)
else y[j]+h/2*(3*f(x[j],y[j])-f(x[j-1],y[j-1])) fi
end: ## Adams-Bashforth method 2nd order
## Adams-Bashforth method 3rd order
AB3 := proc (x,y,j,h,metoda2)
if j<2 then metoda2(x[j],y[j],h)
else y[j]+h/12*(23*f(x[j],y[j])-16*f(x[j-1],y[j-1])+5*f(x[j-2],y[j-2])) fi
end: ## Adams-Bashforth method 3rd order
## Adams-Bashforth method 4th order
AB4 := proc (x,y,j,h,metoda2)
if j<3 then metoda2(x[j],y[j],h)
else y[j]+h/24*(55*f(x[j],y[j])-59*f(x[j-1],y[j-1])+37*f(x[j-2],y[j-2])
-9*f(x[j-3],y[j-3])) fi
end: ## Adams-Bashforth method 4th order
Predictor-corrector methods (PECE):
metoda2 ... starting single-step method
metoda3, metoda4 ... describe the predictor as in procedure 'res'
## Adams-Moulton method 2nd order
AM2 := proc (x,y,j,h,metoda2,metoda3,metoda4) local ypred;
ypred:=metoda3(x, y, j, h, metoda4);
userinfo(6, res, `x =`, x[j+1], `ypred =`, ypred);
y[j]+h/2*(f(x[j+1],ypred)+f(x[j],y[j]))
end: ## Adams-Moulton method 2nd order
## Adams-Moulton method 3rd order
AM3 := proc (x,y,j,h,metoda2,metoda3,metoda4) local ypred;
if j<1 then metoda2(x[j],y[j],h)
else ypred:=metoda3(x, y, j, h, metoda4);
userinfo(6, res, `x =`, x[j+1], `ypred =`, ypred);
y[j]+h/12*(5*f(x[j+1],ypred)+8*f(x[j],y[j])-f(x[j-1],y[j-1]))
fi
end: ## Adams-Moulton method 3rd order
## Adams-Moulton method 4th order
AM4 := proc (x,y,j,h,metoda2,metoda3,metoda4) local ypred;
if j<2 then metoda2(x[j],y[j],h)
else ypred:=metoda3(x, y, j, h, metoda4);
userinfo(6, res, `x =`, x[j+1], `ypred =`, ypred);
y[j] + h/24*(9*f(x[j+1],ypred) + 19*f(x[j],y[j]) - 5*f(x[j-1],y[j-1]) + f(x[j-2],y[j-2]))
fi
end: ## Adams-Moulton method 4th order
Standard procedures of Maple:
## symbolic solution
MSymb := proc (xk, k) local g,gg,j; global x,y,h;
gg := dsolve({diff(Y(X),X)=f(X,Y(X)), Y(x0)=y0}, Y(X), explicit);
userinfo(16, res, `exact solution:`, gg);
g := unapply(op(2, gg), X);
if nargs>1 then
h:=(xk-x0)/k;
x:=array(0..k); for j from 0 to k do x[j]:=x0+j*h od;
y:=array(0..k);
for j from 0 to k
do y[j]:=g(x[j]);
userinfo(4, res, `x =`, evalf(x[j]), `y =`, evalf(y[j]))
od;
fi;
g(xk);
end: ## symbolic solution
## numerical solution
MNum := proc (xk, k)
## 3rd argument is one method from: rkf45, dverk78, classical
local g,gg,j,metoda; global x,y,h;
if nargs<3 then metoda:=rkf45 else metoda:=args[3] fi;
gg := dsolve({diff(Y(X),X)=f(X,Y(X)), Y(x0)=y0}, Y(X),
type=numeric, method=metoda);
g := unapply('op(2, op(2, gg(t)))', t);
if nargs>1 then
h:=(xk-x0)/k;
x:=array(0..k); for j from 0 to k do x[j]:=x0+j*h od;
y:=array(0..k);
for j from 0 to k
do y[j]:=g(x[j]);
userinfo(4, res, `x =`, x[j], `y =`, y[j])
od;
fi;
g(xk);
end: ## numerical solution
Display points:
Use additional plot options (style = ..., color = ...).
plotxy := proc(x,y,k) local j;
plot([seq([x[j],y[j]], j=0..k)], seq(args[j], j=4..nargs))
end:
Plot options for all following graphical commands.
plots[setoptions](thickness=2);
Graphical representation
pole:=DEtools[DEplot]({diff(Y(X),X)=f(X,Y(X))}, {Y(X)}, X=x0..xk, [[Y(x0)=y0]], Y=-2..15):
If you leave out 4th parameter [Y(x0) = y0] the automatic solution will be not displayed. Use additional options DEplot.
presres:=DEtools[DEplot]({diff(Y(X),X)=f(X,Y(X))}, {Y(X)}, X=x0..xk, [[Y(x0)=y0]], Y=-2..15, arrows='none'):
pole;
Calculation and examples
Demonstration examples
k:=5;
Single-step method (Euler)
res(xk, k, RK, Euler);
Explicit method (Adams-Bashforth method 2nd order, start - Euler's Method - 1st modification)
res(xk, k, AB2, EulerMod1);
Predictor-corrector method (corrector - Adams-Moulton method 4th order, start - Runge-Kutta method 3rd order; predictor - Heune's method)
res(xk, k, AM4, RungeKutta3, RK, EulerMod2);
Symbolic solution
MSymb(xk);
MSymb(xk, k); # with the possibility of additional informations
Numerical solution
MNum(xk);
MNum(xk, k); # with ability additionaly informations
MNum(xk, k, dverk78); # set other methods, see Help
Graphical demonstration of the last solution
plots[setoptions](thickness=1, font=[TIMES, ROMAN, 5]);
Tile vertically
reseni:=plotxy(x, y, k, style=POINT, symbol=POINT, color=BLUE, title="PECE", titlefont=[TIMES, ROMAN, 5]):
p1 := plots[display]([presres, reseni]):
p2 := plots[display]([pole, reseni]):
plots[display](array(1..2, [p1, p2]));
Go to demonstration examples.