MATLAB: Equation: X”-6x’+13x = t+3sin(t) Initial Value: x(0)=1 t є [0,1] Method: Runge-Kutta II Step Sizes: h=0.1 , h=0.03

matlark2runge kutta 2

I want to solve it by the Matlab only. but facing the Problem . Can someone Please help me out?

Best Answer

The first thing you need to do is to write the ode as two first order equations and use the code below. You will be required to supply two initial conditions for the 1s order equations. Use the one that you are given plus another of your choice.
function [t,x,y,N] = Runge2_2eqs(f1,f2,to,tfinal,xo,yo,h)
     % This function implements the Rk2 method. 
     t    =  to;
     N = ceil((tfinal-to)/h);
     x = zeros(1,N); 
     y = zeros(1,N) ;
     x(1) = xo;
     y(1) = yo;
       for i = 1:N
          t(i+1) = t(i)+h;
          Sx1 = f1(t(i),x(i),y(i));
          Sy1 = f2(t(i),x(i),y(i));
          Sx2 = f1(t(i)+h, x(i)+Sx1*h, y(i)+Sy1*h);
          Sy2 = f2(t(i)+h, x(i)+Sx1*h, y(i)+Sy1*h);
         x(i+1) = x(i) + h/2*(Sx1+Sx2);
         y(i+1) = y(i) + h/2*(Sy1+Sy2);
       end
end
This is the mfile.
xo = 1;
yo = 0;
h = [.1 0.03];
to = 0;
tfinal = 20;
M = ceil((tfinal-to)/h(2));
dx1 =  @(t,x1,x2) x2;
dx2 =  @(t,x1,x2) 6*x2 -13*x1 + t + 3*sin(t);
% When you reduce the equation to two first order, x will be the solution
% of the ode, i.e x'' and y is its derivative, x'.
for i = 1: length(h)
if (i== 1)   % This for the case when h = 0.1
 [t,x,y,N] = Runge2_2eqs(dx1,dx2,to,tfinal,xo,yo,h(i));
   y1  = x;
   y2  = y;
else       % and for the case when h = 0.03
   [t,x,y,N] = Runge2_2eqs(dx1,dx2,to,tfinal,xo,yo,h(1));
    x3 = x;   
    x4 = y;
end
end
t1 = t(1):(t(end)-t(1))/(M-1):t(end);
figure(1);
subplot(121)
plot(t1,y1, '-o')
hold on
plot(t1,y2,'-o')
legend('Dfx1','Dfx2')
title('Solution to two systems of ODEs using RK2, h= 0.1')
xlabel('x')
ylabel('y')
xlim([to tfinal])
grid
subplot(122)
plot(t,x3,'linewidth',2,'color','b')
hold on
plot(t,x4,'linewidth',2,'color','r')
legend('Dfx1','Dfx2')
title('Solution to two systems of ODEs using RK2, h = 0.03')
xlabel('x')
ylabel('y')
xlim([to tfinal])
grid
% Using ode 45 just to prove that the solution with RK2 is correct.
F = @(t,y) [ y(2); (6*y(2) -13*y(1) + t + 3*sin(t)) ];
t0 = 0;
tf = 20;
delta = (tf-t0)/(201-1); 
tspan = t0:delta:tf;
ic = [1 0];
[t,y] = ode45(F, tspan, ic);
figure
plot(t,y(:,1),'-o')
hold on
plot(t,y(:,2),'-o')
a = title('Using ode45');
legend('x','x_{prime}');
set(a,'fontsize',14);
a = ylabel('y');
set(a,'Fontsize',14);
a = xlabel('t [0 20]');
set(a,'Fontsize',14);
xlim([t0 tf])
grid
grid minor;
Related Question