What is wrong with this code? Need help.

Main code:
clear all
close all
clc
M = [1 0 0
0 1 0
0 0 1];
x0=[-0.5; 0.2; 0.5];
options = odeset('Mass',M,'RelTol',1e-12,'AbsTol',[1e-14 1e-14 1e-14], 'Vectorized','on');
global t x y z dt alpha
dt=0.01;
for alpha=0.8800:0.0005:1.6000
alpha
clear n
clear m
[t,x]=ode15s(@equations,0:dt:500,x0);
n=length(x(:,1));
m=floor(n/2);
y=diff(x(m,n,1))/dt; %% Error Message at this line: Index exceeds matrix dimensions
z=diff(y)/dt;
k=1;
clear aa
for i=m:n
t0=t(i); %Comp. local max. pts for m<t<n
option=optimset('display','off');
zer=fsolve(@differ,t0,option);
if interp1(t(m:n-2),z,zer)>0
aa(k)=interp1(t(m:n),x(m:n,1),zer);
k=k+1;
end
end
kmax=k-1;
h=plot(alpha.*ones(1,kmax),aa,'r.');
hold on
set(h,'MarkerSize',0.1);
end
Function equations.m:
function xdot=equations(t,x)
global alpha
a=0.0005; b=0.01; eps=0.01; beta=-1; R=0.3;
xdot(1) = (-x(2) + alpha*x(1)^2 + beta*x(1)^3)/eps;
xdot(2) = x(1) - x(3) - R*x(2);
xdot(3) = a - b*x(2);
xdot = xdot';
end
Function differ.m:
function f=differ(a)
global t x dt
s=diff(x(m:n,1))/dt;
f=interp1(t(m:n-1),s,aa);
end
The author of this program plot the attached picture. But when I tried cannot get that plot. Please help!!!!

回答(1 个)

y=diff(x(m,n,1))/dt; %% Error Message at this line: Index exceeds matrix dimensions
Yes, see the reason:
>> m
m =
25000
>> n
n =
50001
>> whos x
Name Size Bytes Class Attributes
x 50001x3 1200024 double global
The size of x is 5000x1 and you are trying to acess of x as follows (In first iterations)
x(25000,50001,1)
%..^.........^ no x values
Thats why it shows the error index exceeds matrix dimensions. That means you trying to acess the indices more than it have.

3 个评论

Respected Sir,
Can you fix it and plot the graph as attached in the pdf.
Thanks
Considering your function differ() uses
s=diff(x(m:n,1))/dt;
then I speculate that your earlier line
y=diff(x(m,n,1))/dt;
should be
y=diff(x(m:n,1))/dt;
Respected Sir,
I did tried for that also. It gives error. Can you help and plot. Plot is given in pdf file.
Thanks

请先登录,再进行评论。

类别

帮助中心File Exchange 中查找有关 Mathematics 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by