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)

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.