### 4.21 How to implement Runge-Kutta to solve diﬀerential equations?

4.21.1 First order ODE
4.21.2 second order ODE

#### 4.21.1 First order ODE

Solve $$\frac {dx}{dt} = x(t) \sin (t)$$ with initial conditions $$x(0)=\frac {1}{10}$$. The exact solution is $$\frac {1}{10} e^{1-\cos (t)}$$

Find $$x(t)$$ for 5 seconds, using $$\Delta _t=0.001$$

Matlab

function [x,t] = nma_RK4_2(  )
%solve dx/dt= x*sin(t); with x(0)=1/10 which has solution
%x(t)=1/10*exp(1-cos(t));

delT = 0.001; %time grid
t    = linspace(0,5,round(5/delT)); %time
N    = length(t);
x    = zeros(N,1);
x(1) = 0.1;  %initial conditions
for i = 1:(N-1)
k1     = delT * f( t(i)          , x(i) );
k2     = delT * f( t(i)+1/2*delT , x(i)+1/2*k1 );
k3     = delT * f( t(i)+1/2*delT , x(i)+1/2*k2 );
k4     = delT * f( t(i)+delT     , x(i)+k3 );
x(i+1) = x(i)+1/6*(k1+2*k2+2*k3+k4);
end

function r=f(t,x) %change RHS as needed
r=x*sin(t);
end
end


Now call the above function and plot the solution

[x,t]=nma_RK4_2();
plot(t,x)



#### 4.21.2 second order ODE

To solve second (and higher order ODE’s), we have to ﬁrst rewrite the ODE as set of ﬁrst order ODE and only then implement RK similarly to the above. There will be a $$k_i$$ parameters for each ﬁrst order ODE in the set. For this example, since we have two ﬁrst order ode’s (since it is second order ODE), we can call the parameters for the ﬁrst oder ODE $$k_i$$ and for the second order ODE $$z_i$$. For higher order, we can use a matrix to keep track of the RK parameters.

Solve $$x''(t)=x(t) + x'(t)^2 - 5 t x'(t)$$ with initial conditions $$x(0)=2,x'(0)=0.01$$. The right hand side is $$f(x,x'(t),t)=x(t)+ x'(t)^2-5 t x'(t)$$. This is non-linear second order ode.

The ﬁrst step is to convert this to two ﬁrst order odes. Let $$x=x(t)$$ and $$v(t)=x'(t)$$, hence $$x'(t)=v(t)$$ and $$v'(t)=x''(t)=f(x,v,t)$$.

Matlab

function [x,v,t] = nma_RK4_3(  )
%solve x''(t)=x - 5 t x'(t) + (x'(t))^2 with x(0)=2,x'(0)=0.01,

delT = 0.001; %time grid
t    = linspace(0,5,round(5/delT)); %time
N    = length(t);
x    = zeros(N,1);
v    = zeros(N,1);
x(1) = 2;  %initial conditions
v(1) = 0.01;

for i = 1:(N-1)
k1     = delT * v(i);
z1     = delT * f( t(i)          , x(i)        , v(i) );

k2     = delT * (v(i) + 1/2* z1);
z2     = delT * f( t(i)+1/2*delT , x(i)+1/2*k1 , v(i)+1/2*z1);

k3     = delT * (v(i) + 1/2* z2);
z3     = delT * f( t(i)+1/2*delT , x(i)+1/2*k2 , v(i)+1/2*z2);

k4     = delT * (v(i) + z3);
z4     = delT * f( t(i)+delT     , x(i)+k3      , v(i)+z3);

x(i+1) = x(i)+1/6*(k1+2*k2+2*k3+k4);
v(i+1) = v(i)+1/6*(z1+2*z2+2*z3+z4);
end

function r=f(t,x,v)
r=x - 5*t*v + v^2;
end
end


The above function is called the solution is plotted

[x,v,t]=nma_RK4_3();
plot(t,x);
plot(t,v);



Here is plot result

A better way to do the above, which also generalized to any number of system of equations is to use vectored approach. Writing $\left [\dot {x}\right ] = \left [f(t,x,x',\dots )\right ]$ Where now $$[x]$$ is the derivative of state vector and $$\left [f(t,x,x',\dots )\right ]$$ is the right hand side, in matrix form. The above second order is now solved again using this method. This shows it is simpler approach than the above method.

function [x,t] =nma_RK4_4(  )
%solve x''(t)=x - 5 t x'(t) + (x'(t))^2 with x(0)=2,x'(0)=0.01,

delT = 0.001; %time grid
t    = linspace(0,5,round(5/delT)); %time
N    = length(t);
x     = zeros(2,N);
x(1,1) = 2;  %initial conditions
x(2,1) = 0.01;

for i = 1:(N-1)
k1     = delT * f( t(i)          , x(:,i));
k2     = delT * f( t(i)+1/2*delT , x(:,i)+1/2*k1);
k3     = delT * f( t(i)+1/2*delT , x(:,i)+1/2*k2);
k4     = delT * f( t(i)+delT     , x(:,i)+k3);

x(:,i+1) = x(:,i)+1/6*(k1+2*k2+2*k3+k4);
end

function r=f(t,x)
r=[x(2)
x(1)-5*t*x(2)+x(2)^2];
end
end


The above is called as follows

[x,t]=nma_RK4_4();
plot(t,x(1,:)) %plot the position x(t) solution
plot(t,x(2,:)) %plot the velocity x'(t) solution



Same plots result as given above.