Index Exceeds Matrix on Runge Kutta Method

Currently trying to create a four order Runge Kutta script to help solve a forced linear oscillator equation x¨(t) + bx˙(t) + x(t) = cos(Ω0t) I keep getting the error Index exceeds matric dimensions. Error in the Function and also an error in the first order of runge kutta line. Any help will be appreciated thanks.

 % Define Parameters 
h=0.1;
Omega=0.5;
Tend=50;
nsteps=Tend/h;
m = 1.0;
b = 0.0;
k = 1.0;
F = 0.0;
T=1;
x¨(t) + bx˙(t) + x(t) = cos(Ω0t)
%Inital conditions
y(1)=1;
y(2)=500;
t(1)=1;
%Define Function Handle
f= @(t,y) [y(2); cos(Omega*T)-b*y(2)-t(1)];
%Function
for i=1:ceil(Tend/h)

    k1=f(t(i),      y(i)         );
    k2=f(t(i)+0.5*h,y(i)+0.5*k1*h);
    k3=f(t(i)+0.5*h,y(i)+0.5*k2*h);
    k4=f(t(i)    +h,y(i)+    k3*h);
    y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);
end
%Plot Results
plot(t,y)
tlabel('t')
ylabel('y')

end

[Eamon]

The Ω in the line

x¨(t) + bx˙(t) + x(t) = cos(Ω0t)

is undefined. You could try replacing Ω with "omega." Also, here's an example of some other 4th order Runge Kutta method code that might be useful.

[Adam Valli]

I have changed that, thank you.

Answer 1

[James Tursa]

I took your code and made some changes to get it to run and produce a plot. Changes are noted in comments:

   % Define Parameters 
  h=0.1;
  Omega=0.5;
  Tend=50;
  nsteps=Tend/h;
  m = 1.0;
  b = 0.0;
  k = 1.0;
  F = 0.0;
  T=1;
  %x¨(t) + bx˙(t) + x(t) = cos(Ω0t)
  %Inital conditions
  y = zeros(2,nsteps+1);  % <-- made this a matrix
  y(1,1)=1;  % <-- 1st column of y matrix is the initial condition
  y(2,1)=500;
  t = linspace(1,Tend,nsteps+1);  % <-- made this a vector
  %Define Function Handle
  f= @(t,y) [y(2); cos(Omega*T)-b*y(2)-y(1)];  % <-- changed last t(1) to y(1) to match your DE
  %Function
  for i=1:nsteps  % <-- Changed indexing limits
      k1=f(t(i),      y(:,i)         );  % <-- changed y(i) to y(:,i) since the "state" is a column vector
      k2=f(t(i)+0.5*h,y(:,i)+0.5*k1*h);
      k3=f(t(i)+0.5*h,y(:,i)+0.5*k2*h);
      k4=f(t(i)    +h,y(:,i)+    k3*h);
      y(:,i+1)=y(:,i)+h/6*(k1+2*k2+2*k3+k4);  % <-- and the "next state" is y(:,i+1)
  end
  %Plot Results
  plot(t,y(1,:))  % <-- plot only the y part, not the y' part
  xlabel('t')  % <-- changed tlabel to xlabel
  ylabel('y')

FYI, the error you were getting was because you were passing in y(i) (a single scalar) to your f function, but the f function assumes that the y passed in is a 2-element vector. Hence the change to make y a 2xN matrix and treat the columns of y as the states at any particular time.

[Adam Valli]

Thank you. I see where I've gone wrong now thanks to your explanation.