Index Exceeds Matrix on Runge Kutta Method

4 次查看（过去 30 天）
Adam Valli
Adam Valli 2016-10-13
评论： Adam Valli 2016-10-20
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
  2 个评论
Eamon
Eamon 2016-10-13
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

James Tursa
James Tursa 2016-10-13
编辑：James Tursa 2016-10-13
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.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Runge Kutta Methods 的更多信息

标签

尚未输入任何标签。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by