MATLAB: Runge-Kutta 3

errorsMATLABrunge kutta

So I'm doing a project for school, where i have to solve the x''-x'-x=sin(t) equation with Runge-Kutta.

intervalmin=0;

intervalmax=1;

numnodes=5;

f=@(t,x) [x(2); x(2)+x(1)+sin(t)];

inival=0;

h=0.1;

t=zeros(1,numnodes);

t(1)=intervalmin;

x=zeros(1,numnodes);

X=zeros(1,numnodes);

x(1)=inival;

X(1)=inival;

for i = 2:numnodes

k1 = f(t(i-1),x(i-1));

k2 = f(t(i-1)+h/2,x(i-1)+(h/2)*k1);

k3 = f(t(i-1)+h,x(i-1)-h*k1+2*h*k2);

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

X(i)=(X(i-1)+h/6)*((f(t(i-1),X(i-1)))+4*(f(t(i-1)+h/2,X(i-1)+(h/2)*k1))+(f(t(i-1)+h,X(i-1)-h*k1+2*h*k2)));

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

end

figure

plot(t,x,'.-b',t,X,'-.r');

The code is this so far, but it has multiple errors like:

Index exceeds the number of array elements (1).

Error in Untitled5>@(t,x)[t(2);t(2)+t(1)+sin(x)]

Error in Untitled5 (line 14) k1=f(t(i-1),x(i-1));

I'm pretty new to matlab, and just can't seem to find a way to make it work.

Best Answer

Here is the mistake

SOlution

Related Solutions

MATLAB: Second Order ODE with Runge Kutta 3 “K’s problem”

You are using the same k's for x1 and x2, which is incorrect. That is, you are using f2( ) to calculate the k's, but these should only apply to the 2nd derivative, not to the 1st derivative. Your 1st derivative function f1( ) isn't even used, but it should be.

I would suggest you set up your state as a 2-element vector instead of separate variables x1 and x2. This will help you to write vectorized code for your RK4 scheme, and will also match what you would need to do when you move to the MATLAB function ode45( ). So, only have one derivative function, f( ), and have it work with a 2-element state vector. E.g.,

f = @(t,x) [ x(2); (-cos(t))./((1+sin(t)).^2) ];

Then use this to calculate your k's and to update your state at each step. Instead of individual scalars at each step, you will have 2-element vectors at each step. Maybe store them by columns. E.g.,

This

x1=zeros(1,numnodes);
x2=zeros(1,numnodes);

becomes this

x=zeros(2,numnodes); % 2-element state vectors stored by columns

And this

x1(1)=inival1;
x2(1)=inival2;

becomes this

x(:,1)=[inival1;inival2];

And this

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

becomes this

k1 = f(t(i-1),x(:,i-1));

And this

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

becomes this

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

etc.

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

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;

Best Answer

Related Solutions

MATLAB: Second Order ODE with Runge Kutta 3 “K’s problem”

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

Related Question