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

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

[Walter Roberson]

What problem are you encountering? What is your code, and what error message are you encountering?

[Tariq Malik]

This is the code that I am using.

intmin=0;
intmax=1;
numnodes=5; 
f=@(t,x) t+3*(sin(t));
inival=1;
h=(intmax-intmin)/(numnodes-1);
t=zeros(1,numnodes);
x=zeros(1,numnodes);
X=zeros(1,numnodes);
t(1)=intmin;
x(1)=inival;
X(1)=inival;
for i=2:numnodes; 
    t(i)=t(i-1)+h; 
    k1=f(t(i-1),x(i-1)); 
    k2=f(t(i-1)+h/2,x(i-1)+(h/2)*k1);
    x(i)= x(i-1)+h*k2;
    X(i)=X(i-1)+h*f(t(i-1),X(i-1));
end
figure
plot(t,x,'.-',t,X,'-*')
hold on 
syms t x(t) sym(f) 
eqn=diff(x,t)==f(t,x(t));
cond1=x(intmin)==inival;
odesol=dsolve(eqn,cond1);
odesolfun=@(t) eval(odesol);
tt=linspace(intmin,intmax,100*(ceil(intmax-intmin)));
xx=odesolfun(tt);
plot(tt,xx,'r')
hold off

[Tariq Malik]

Am I making the Problem in regarding to typing the Equation?

Answer 1

[Abraham Boayue]

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;

[Tariq Malik]

intmin=0;

intmax=1;

inival1=0;

inival2=0;

numnodes = 10;

t(1) = intmin;

x1(1)= inival1;

x2(1)= inival2;

t =zeros(1,numnodes);

x1=zeros(1,numnodes);

x2=zeros(1,numnodes);

h=(intmax-intmin)/(numnodes-1);

f1 = @(t,x1,x2) x2;

f2 = @(t,x1,x2) 6*x2-13*x1+t+3*sin(t);

numnodes1=300; 

a=zeros(1,numnodes1);

b1=zeros(1,numnodes1);

b2=zeros(1,numnodes1);

a(1)=intmin;

b1(1)=inival1;

b2(1)=inival2;

g=(intmax-intmin)/(numnodes1-1);

F1 = @(a,b1,b2) b2;

F2 = @(a,b1,b2) 6*b2-13*b1+a+3*sin(a);

for i = 2:numnodes 

t(i) = t(i-1)+h; 

k1 = f2(t(i-1),x1(i-1),x2(i-1)); 

k2 = f2(t(i-1)+h/2,x1(i-1)+(h/2)*k1,x2(i-1)+(h/2)*k1); 

k3 = f2(t(i-1)+h/2, x1(i-1)+(h/2)*k2,x2(i-1)+(h/2)*k2); 

x1(i) = x1(i-1)+(h/6)*(k1+4*k2+k3); 

x2(i) = x2(i-1)+(h/6)*(k1+4*k2+k3); 

end

for i = 2:numnodes1

a(i) = a(i-1)+g; 

k1 = F2(a(i-1),b1(i-1),b2(i-1));

   k2 = F2(a(i-1)+g/2,b1(i-1)+(g/2)*k1,b2(i-1)+(g/2)*k1);

   k3 = F2(a(i-1)+g/2, b1(i-1)+(g/2)*k2,b2(i-1)+(g/2)*k2); 

b1(i) = b1(i-1)+(g/6)*(k1+4*k2+k3);

   b2(i) = b2(i-1)+(g/6)*(k1+4*k2+k3); 

end 

figure

plot(t,x1,'.-',t,x2,'-*')

hold on

syms t x(t)

sym(x)

eqn=diff(x,t,2) == diff(x,t)-13*x+t+3*sin(t);

cond1==1;

cond2==0;

odesol=dsolve(eqn,cond1,cond2);

odesolfun=@(t) eval(odesol);

tt=linspace(intmin,intmax,100*(ceil(intmax-intmin)));

xx=odesolfun(tt);

legend("Runge-Kutta 2", "0.1", "Exact"); 

plot(tt,xx,'r')

hold off