Numerical Methods of Solving 1st order System of ODEs:  Runga-Kutta Method

The initial value problem

Consider the system of DEs 

x’=2x+y
y’=x+2y
x(0)=3, y(0)=-2

This is easy to find the exact solution, so we do not need to use numerical techniques.

Here’s the eigensystem method solution

Clear[MatA]
MatA={{2,1},{1,2}};
Eigensystem[MatA]

Clear[c1,c2]
c1{1,1}Exp[3t] +c2{-1,1}Exp[1t]

{-c2 (E^t)+c1 (E^(3 t)),c2 (E^t)+c1 (E^(3 t))}/.{t->0}

Solve[{c1-c2,c1+c2}=={3,-2}]

{-c2 (E^t)+c1 (E^(3 t)),c2 (E^t)+c1 (E^(3 t))}/.{c1->(1/2),c2->-(5/2)}

This last function is the solution to the IVP.

Here’s the automatic solution.

Clear[t,x, y]
DSolve[{x'[t]==2x[t]+y[t],y'[t] == x[t]+2y[t]},{x[t], y[t]},t]

Expand[{{x[t]->(1/2) (E^t) (1+(E^(2 t))) C[1]+(1/2) (E^t) (-1+(E^(2 t))) C[2],y[t]->(1/2) (E^t) (-1+(E^(2 t))) C[1]+(1/2) (E^t) (1+(E^(2 t))) C[2]}}]

The unique solution to the initial value problem:

DSolve[{x'[t]==2x[t]+y[t],y'[t] == x[t]+2y[t],
x[0]==3,y[0]==-2
},{x[t], y[t]},t]

Expand[{(1/2) (E^t) (5+(E^(2 t))),(1/2) (E^t) (-5+(E^(2 t)))}]

Clear[xsoln,ysoln]
xsoln[t_]:=((5 (E^t))/2)+((E^(3 t))/2)
ysoln[t_]:=((-5 (E^t))/2)+((E^(3 t))/2)

Clear[p1x]
p1x=Plot[xsoln[t],{t,0.0,5.0}]

Clear[p1y]
p1y=Plot[ysoln[t],{t,0.0,5.0}]

ParametricPlot[{xsoln[t],ysoln[t]},{t,0.0,5.0}]

iterations with  step size h=.05

Clear[n,h,
t,x,y,
fx,fy,
m1x,m2x,m3x,m4x,
m1y,m2y,m3y,m4y,
numIterations]

t[0]=0.0;
x[0]=3.0;
y[0]=-2.0;

h=.1;
(*  x'=2x+y
    y'=x+2y  *)

fx[x_,y_,t_]:=2x + y
fy[x_,y_,t_]:=x + 2y

numIterations = 10;

For[n=0,n≤numIterations,

m1x = fx[ x[n], y[n],t[n] ];
m2x  = fx[ x[n] + h/2 , y[n] + m1x*h/2,t[n]];
m3x  = fx[ x[n] + h/2 , y[n] + m2x*h/2,t[n]];
m4x = fx[ x[n] + h , y[n] + m3x*h,t[n]];

m1y = fy[ x[n], y[n],t[n]];
m2y  = fy[ x[n] + h/2 , y[n] + m1y*h/2,t[n]];
m3y  = fy[ x[n] + h/2 , y[n] + m2y*h/2,t[n]];
m4y  = fy[ x[n] + h , y[n] + m3y*h,t[n]];

n=n+1;
t[n]=t[n-1] + h;
x[n] = x[n-1] + (m1x+2*m2x+2*m3x+m4x)*h/6;
y[n] = y[n-1] + (m1y+2*m2y+2*m3y+m4y)*h/6     ]

Clear[pointlistT]
pointlistT=Table[{t[i],x[i],y[i]},{i,0,n}]
Clear[pointlist]
pointlist=Table[{x[i],y[i]},{i,0,n}]

Clear[p2]
p2=ListPlot[pointlist,Joined->True,PlotStyle->{Dashing[{0.15,0.05}],Hue[0.05]}]

ParametricPlot[{xsoln[t],ysoln[t]},{t,0.0,1.1},AxesOrigin->{0,0}]

Compare with the known solution

Because this is an instructional example, we calculate the actual value and relative error.

TableForm[Table[{y[i],soln[x[i]], Abs[(soln[x[i]] - y[i])/soln[x[i]] ] },{i,1,n}]]

Clear[pointlist]
pointlist=Table[{x[i],y[i]},{i,0,n}]

Clear[p2]
p2=ListPlot[pointlist,Joined->True,PlotStyle->{Dashing[{0.15,0.05}],Hue[0.05]}];
Show[{p1,p2}]

Compare with direction field

Even when the DE is hard or impossible to solve, we can compare it to the direction field.

Clear[p3]
p3=VectorPlot[{fx[x,y,0],fy[x,y,0]},{x,-0.5,5.5},{y,-0.5,4.5}];
Show[{p1,p2,p3},Axes->True]

Show[{p1,p2,p3},Axes->True,PlotRange->{{-0.5,5.5},{-0.5,4.5}},AxesOrigin->{0,0}]