Function Issue with ODE45

Hello,

I have this function which I am trying to pass to ODE45. I do not understand why my X array is all zeros. Does anyone know why? Thank you.

clear
clc
 
tspan = [0,6];
options = odeset('RelTol', 1e-7, 'AbsTol', 1e-7);
 
%I.C's
x0= [0;0;0];
rdot0= [0;0;0];
psi0= [0;0;0];
phi0= [0;0;0];
theta0=[0;0;0];
p0= [0;0;0];
q0= [0;0;0];
r0= [0;0;0];
 
[t,X] = ode45(@xdot,tspan,[x0,rdot0,psi0,theta0,phi0,p0,q0,r0],options);
plot(X)
 
 
 
 
function x = xdot(t,x)
 
g = 9.8;
m = 0.450;
l = 0.225;
k = 2.98e-6;
b = 1.14e-6;
 
omegaih=4*m*g;
omega1h=m*g;
omega2h=m*g;
omega3h=m*g;
omega4h=m*g;
 
if t <= 1
    omega1 = omegaih + 70*sin((2*pi*t)/4);
    omega2 = omegaih + 70*sin((2*pi*t)/4);
    omega3 = omegaih + 70*sin((2*pi*t)/4);
    omega4 = omegaih + 70*sin((2*pi*t)/4);
    
elseif t <= 2
    omega1 = omegaih - 77*sin((2*pi*t)/4);
    omega2 = omegaih - 77*sin((2*pi*t)/4);
    omega3 = omegaih - 77*sin((2*pi*t)/4);
    omega4 = omegaih - 77*sin((2*pi*t)/4);
 
elseif t <= 3
    omega1=m*g;
    omega2 = sqrt(omega2h^2 - 70^2*sin((2*pi*(t-2))/4));
    omega3=m*g;
    omega4 = sqrt(omega4h^2 + 70^2*sin((2*pi*(t-2))/4));
 
elseif t <= 4
    omega1=m*g;
    omega2 = sqrt(omega2h^2 + 70^2*sin((2*pi*(t-2))/4));
    omega3=m*g;
    omega4 = sqrt(omega4h^2 - 70^2*sin((2*pi*(t-2))/4));
 
elseif t <= 5
    omega1 = sqrt(omega1h^2 - 70^2*sin((2*pi*(t-4))/4));
    omega2=m*g;
    omega3 = sqrt(omega3h^2 + 70^2*sin((2*pi*(t-4))/4));
    omega4=m*g;
    
elseif t <= 6
    omega1 = sqrt(omega1h^2 + 70^2*sin((2*pi*(t-4))/4));
    omega2=m*g;
    omega3 = sqrt(omega3h^2 - 70^2*sin((2*pi*(t-4))/4));
    omega4=m*g;
 
end
 
xdot=zeros(12,1);
 
T = 4*k*(omega1^2 +omega2^2 + omega3^2 +omega4^2);
 
 
%State Space
xdot(1) = x(4);
xdot(2) = x(5);
xdot(3) = x(6);
xdot(4) = (T*(sin(x(7))*sin(x(9)) + cos(x(7))*cos(x(9))*sin(x(8))))/m;
xdot(5) = -(T*(cos(x(9))*sin(x(7)) - cos(x(7))*sin(x(8))*sin(x(9))))/m;
xdot(6) =  (T*cos(x(7))*cos(x(8)))/m;
 
 
xdot(7) = x(10) + (x(12)*cos(x(7))*sin(x(8)))/(cos(x(8))*cos(x(7))^2 + cos(x(8))*sin(x(7))^2) + (x(11)*sin(x(7))*sin(x(8)))/(cos(x(8))*cos(x(7))^2 + cos(x(8))*sin(x(7))^2);
xdot(8) = (x(11)*cos(x(7))/(cos(x(7)))^2 + sin(x(7))^2) - (x(12)*sin(x(7)))/(cos(x(7))^2 + sin(x(7))^2);
xdot(9) = (x(12)*cos(x(7)))/(cos(x(8))*cos(x(7))^2 + cos(x(8))*sin(x(7))^2) + (x(11)*sin(x(7)))/(cos(x(8))*cos(x(7))^2 + cos(x(8))*sin(x(7))^2);
 
xdot(10) =(3166346726764097*omega4^2)/4722366482869645213696 - (79*x(8)*x(9))/20000 - (3166346726764097*omega2^2)/4722366482869645213696;
xdot(11) =(3166346726764097*omega3^2)/4722366482869645213696 + (79*x(7)*x(9))/20000 - (3166346726764097*omega1^2)/4722366482869645213696;
xdot(12) =(1345874447617849*omega1^2)/1180591620717411303424 + (1345874447617849*omega3^2)/1180591620717411303424 - (1345874447617849*omega2^2)/1180591620717411303424 - (1345874447617849*omega4^2)/1180591620717411303424;
 
 
end
 

Answer 1

[Walter Roberson]

x0= [0;0;0];
rdot0= [0;0;0];
psi0= [0;0;0];
phi0= [0;0;0];
theta0=[0;0;0];
p0= [0;0;0];
q0= [0;0;0];
r0= [0;0;0];

That is 8 variables with 3 values each

[t,X] = ode45(@xdot,tspan,[x0,rdot0,psi0,theta0,phi0,p0,q0,r0],options);

The [x0,rdot0,psi0,theta0,phi0,p0,q0,r0] part is going to construct a 3 x 8 array of initial conditions -- a total of 24 initial conditions

xdot=zeros(12,1);
...
xdot(12) =(1345874447617849*omega1^2)/1180591620717411303424 + (1345874447617849*omega3^2)/1180591620717411303424 - (1345874447617849*omega2^2)/1180591620717411303424 - (1345874447617849*omega4^2)/1180591620717411303424;

You are constructing 12 derivatives, not 24. The number of derivatives constructed must match the number of initial conditions.

if t <= 1
    omega1 = omegaih + 70*sin((2*pi*t)/4);

The mathematics of ode45() and all the other ode*() routines requires that the ode function you pass in must be twice differentiable -- once to construct a Jacobi, and a second time to construct the Hessian. It is not obvious that you have taken care for those implied derivatives to be continuous at each of the time breakpoints in your function.

If your ode function has continuous first derivative (at each of the breakpoints) but not continuous second derivative, then most of the time ode45 will not notice, and will instead produce wrong answers without noticing that the answers are wrong.

If your ode function does not have continuous first derivative (at each of the breakpoitns), then whether ode45 notices or not will depend upon how large the discontinuity in the derivatives is. Sometimes it will be able to treat it as approximately a steep slope and continue, and other times it will emit an error message and stop. In the cases where it manages to continue... the answers will likely be wrong. In such situations, you are lucky if it stops with an error message, because then you will not get fooled into thinking the results were meaningful when they are not.

You should be breaking your system up and calling ode45() once for each segment; you can use the same ode function, but the first tspan should be [0 1], the second [1 2], and so on to [5 6]. You would copy the boundary conditions out of the results of the previous one to become the boundary conditions for the next one.

[Steven Lord]

Once you've modified your function so it returns the same number of elements in the output argument of your xdot function as you have elements in your initial condition vector, call it with the starting time (t = 0) and the initial condition vector. Because your initial condition vector is all 0's I suspect your function will return all 0's.

If you start at the starting line (initial condition all 0's) and leave your engine off (your function returns all 0's), it's going to be a long time before you finish the race.

[John Biscanti]

Thank you for the response @Steven Lord. I am sorry, I do not understand what you mean by call it with the starting time. Do you mean in tspan of ode45?

[John Biscanti]

Also, thank you for the response and help @Walter Roberson. I am implementing your suggestions now. They have helped me make more sense of this now.

[Walter Roberson]

Suppose you call your function passing in time 0 and boundary conditions all 0. Is there any force, or does it sit there not moving? If it sits there then your function will produce all 0 that will become input for the next time: if you have a system sitting still at the start of time 1 millisecond, then do your equations have a time-dependent force that would move the object? or is it going to continue sitting still?