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 first rewrite the ODE as set of first order ODE and only then implement RK similarly to the above. There will be a \(k_i\) parameters for each first order ODE in the set. For this example, since we have two first order ode’s (since it is second order ODE), we can call the parameters for the first 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 first step is to convert this to two first 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.