MATLAB: Problem in simulation ODE45 Index exceeds array bounds.

differential equationsMATLABode45simulation

Hi. I am new in Matlab and I am trying to solve the following differential equation.
function [xdot] = HW5P3(t,x)% x is theta
%Function equations
% x(1)=theta,x(2)=thetadot,x(3)=h,x(4)=lambda
v=12.*sin(2.*pi.*(1./2).*t); %sine source
L=1-(2./pi).*abs(x(1));
hdot=-((1./pi).*(sign(x(1))./(L.^2)).*(x(4)).^2)-10.*x(3)-x(1);
thetadot=1000.*x(3);
lambdadot=v-(100./L).*x(4);
xdot=[hdot;thetadot;lambdadot];
end
[t,x]=ode45('HW5P3',[0 5],[1,0,0]);

Best Answer

 I added more elements to the code and also new conditions but this one worked for me.
    % Function
function [xdot y] = HW5P3(t,x)% x is theta
%Function equations
global a b B J K vo R Lo thetaout
theta=x(1);
h=x(2);
lambda=x(3);
v=(t<0.5).*(0)+((t>=0.5)&(t<4.5)).*(vo*sin(2*pi*(1/0.1)*t))+(t>=4.5).*(0); %Voltage source starts at t=0.5sec
L=(theta<-thetaout).*(Lo)+((theta>=-thetaout)&(theta<thetaout)).*(a+b.*abs(theta))+(theta>=thetaout).*(Lo); %Inductance change with theta
Ldot=b*sign(theta);%Derivative of L
lambdadot=v-R*(lambda/L);
hdot=((Ldot/(2*(L)^2))*lambda^2)-((B/J)*h)-K*theta;
thetadot=h/J;
xdot=[thetadot;hdot;lambdadot];
y=[L;v;Ldot];
end
      % code
global a b B J K vo R Lo thetaout
a=1;
b=-2/pi;
B=0.01;%Damping factor for bearings
J=0.001;%Inertia of the shaft
K=1;%Rotational Spring of the shaft
Lo=0.1;%minimum inductance when the iron of the blade is outside the gap
vo=120;%source voltage pk-pk
R=10; %The resistance value was to high and the source wasn't giving enough current. I changed it to 10ohms
thetaout=1.4137;%limit of the blade to work by the effect of L(theta) out of this angle L=Lo
lambdai=0;%Initial condition
thetai=1.4137;%Initial COndition

hi=0;%Initial COndition
[t,x]=ode45('HW5P3',[0 5],[thetai,hi,lambdai]);  %calling of function
L=length(t);
N = length(t);
for i=1:N %getting values of output from the function for every time point
[xd,y] = HW5P3(t(i),x(i,:));
Ln(i) = y(1);%L(theta)
vs(i) = y(2);%Voltage source
Lndot(i) = y(3);%Derivative of L(theta)
end
theta=x(:,1);
h=x(:,2);
lambda=x(:,3);
thetadot=h/J;
is=lambda/L;
Trotspring=K*theta;
Tbearing=B*thetadot;
figure(1);
subplot(3,1,1)
plot(t,theta), xlabel('time[sec]'),ylabel('Angular displacement [rad]')
subplot(3,1,2)
plot(t,thetadot), xlabel('time[sec]'),ylabel('Angular velocity [rad/sec]')
subplot(3,1,3)
plot(t,Ln,t,theta), xlabel('time[sec]'),legend('Ln [H]','Angular displacement [rad]')
figure(2);
subplot(3,1,1)
plot(theta,Ln), xlabel('Angular displacement [rad]'),ylabel('L [H]')
subplot(3,1,2)
plot(t,theta,t,is*1000), xlabel('time[sec]'),legend('Angular displacement [rad]','Current [mA]')
subplot(3,1,3)
plot(t,Trotspring,t,Tbearing), xlabel('time[sec]'),legend('Spring Torque [N]','Bearings Tor
Related Question