Index in position 1 is invalid. Array indices must be positive integers or logical values. My Events

I am currently working on an inverted pendulum problem. I have used the My event funtion to stop the penduum action at 45degree (please see code and function below). The highlighted line is causing problems. It is my understanding this is the strating conditions of the pendulum, but if i set to [0,0] it comes up with the above error. The current conditions are not correct but i am unable to change them and run the programme.

function zdot=crash(t,z,g,r,m,k,a)
zdot=zeros(2,1);
zdot(1)=z(2);
zdot(2)=(r^2*z(2)^2-r*a*cos(z(1))-r*g*sin(z(1)))/(k^2-2*r^2);
end

function [value,isterminal,direction] = dashboard(t,z)
  minz = (pi/4);   %for example
  value = z(:,1) < minz; 
  isterminal = 1; %stops the integration
  direction = 0; % The zero can be approached from either direction, -1 decreasing function, +1 increasing fucniton
end

clear;clc
%crash conditions
g=9.81; r=0.4; k=0.5; m=40; a=98.1; 
 
tspan1=(0:0.002:5);
z0=[pi/2 0];

 
z0=[0, (a/r)^0.5];  <<<<<<<<<<<<<<<<<<<<----------------------------------------THIS LINE
Opt    = odeset('Events', @dashboard);
[t,z]=ode45(@(t,z)crash(t,z,g,r,k,a,m),tspan1,z0,Opt);
 
ind = interp1(z(:,1),1:length(z(:,1)),pi/4,'nearest');
 
V_g=z(ind,2)*r
 
V_h=z(ind,2)*.75

Answer 1

[Ameer Hamza]

The ode solver does not cause the issue. When z0 = [0 0], the pendulum never reaches the angle of pi/4, and then you run the line

ind = interp1(z(:,1),1:length(z(:,1)),pi/4,'nearest');

and that outputs NaN because pi/4 is not in the range of elements in z(:,1). Then you try to index using NaN value

z(ind,2)*r

and MATLAB's gives an error. You need to consider how to specify the maximum angle (instead of pi/4) for a given initial condition.

[Matt Egan]

Thanks, that was something i was missing but not the main problem. If i set z0=[0,0] it get the error message. For some reason I have to assign a value (such as (a/r)^0.5). Why cant i set my initial conditions to 0, 0?

[Matt Egan]

Or maybe the major problem is the below changes the orignal function entirely, bringing completly differant results. I though this would just enable me to stop the pendulum at 45 degrees;

Opt    = odeset('Events', @dashboard);
[t,z]=ode45(@(t,z)crash(t,z,g,r,k,a,m),tspan1,z0,Opt);

[Ameer Hamza]

Matt, As i already mentioned, this is not an issue with ode45 solver. If you look at the error message, it will point it to this line

V_g=z(ind,2)*r

The dashboard function is correct. This behaviour is the characteristic of your ODE. For example, Do you expect the pendulum to reach pi/4 angle if initial state is [0 0]?

[Matt Egan]

I understand your first point and have adjjusted my calcs and code accordingly. As such the system does reach pi/4 with initial conditions [0,0]. The major problem is that by calling the even function, my original 'crash' function becomes changed. Please try running the below, without the even function and view the 'z' matrix in workspace. Then add the even function and the answers come up completely differant.

clear;clc
 
%crash conditions
g=9.81; r=0.4; k=0.5; m=40; a=98.1; 
 
tspan1=(0:0.002:1);

z0=[0,0];

[t,z]=ode45(@(t,z)crash(t,z,g,r,m,k,a),tspan1,z0);

clear;clc
 
%crash conditions
g=9.81; r=0.4; k=0.5; m=40; a=98.1; 
 
tspan1=(0:0.002:1);

z0=[0,0];

Opt    = odeset('Events', @dashboard);
[t,z]=ode45(@(t,z)crash(t,z,g,r,k,a,m),tspan1,z0,Opt);

Thanks for your help mate.

[Ameer Hamza]

Matt, when I run your first code the ode45 fails at t = 0.848, probably because of a singularity

Warning: Failure at t=8.489838e-01.  Unable to meet integration tolerances
without reducing the step size below the smallest value allowed (1.776357e-15)
at time t. 
> In ode45 (line 360)
  In test_fun (line 8) 

Also, according to the documentation, the event function 'dashboard' should not effect the crash function. It just tell the ode45 when to stop the integration.

Also, note that the two code you have are not exactly same

[t,z]=ode45(@(t,z)crash(t,z,g,r,m,k,a),tspan1,z0);
[t,z]=ode45(@(t,z)crash(t,z,g,r,k,a,m),tspan1,z0,Opt);
                               %^ different inputs

[Ameer Hamza]

It will be useful if you can provide the mathematical form of your equations. It will be easier to suggest a solution.

[Matt Egan]

Attached is the calculations, its been done ussing vector kinematic equations simulating a car crash.

When i put the first one in i get no error codes and the oedulum reaches pi/4 in a time that seems reasonable. The second equation shows an angualar velocity (z(2)) far to tiny to be realistic meaning it never reaches pi/4. The only thing that has changed when i do that is adding the event.

As for rearanging the second ode to match the first it simply wont run. I am so confused

[Matt Egan]

Why does rearanging the inputs have such a great differance?

[Ameer Hamza]

The inputs define the parameters in your ODE, so rearranging one can have a drastic effect. For example, the fifth input to crash function is mass (m), which has a value of 40. But if you pass k as fifth input, the ODE will be solved if the mass is 0.5 and k is 40. It can completely change the solution of the ODE. Try to run both codes with the correct order if input values.

[Matt Egan]

I can see why the first script doenst work for you, it comes up with what i see as the correct value. I have now rearranged the script to the below:

clear;clc
 
%crash conditions
g=9.81; r=0.4;  m=40; k=0.5;  a=98.1; 
 
tspan1=(0:0.002:1);

z0=[0,0];

Opt    = odeset('Events', @dashboard);
[t,z]=ode45(@(t,z)crash(t,z,g,r,m,k,a),tspan1,z0,Opt);

The variables are now all in the correct order (sames as the crash function) however i get this error message and no values for 'z' in the workspace:

Index exceeds the number of array elements (1).

Error in odezero (line 142)
  if any(isterminal(indzc))

Thanks

[Ameer Hamza]

Change the line in your dashboard function from

value = z(:,1) < minz;

to

value = z(1) < minz;

[Matt Egan]

Ameer my man thats done the trick! Really can't thnk you enough for your help.

[Ameer Hamza]

Good to hear that it worked. I am glad to be of help.