Error in ode45 while doing numerical integration

Hello,

I am trying to implement ode45 to integrate my numerical function by following the ODE with Time-Dependent Terms in the ode45 documentation.

I have a function with the form:

Where and are 1xn vectors that change at each time step and is a constant 1xn vector with each element being 9.81 .

As suggested in the documentation, I created the vector for in the main script as follows:

gt = linspace(0,0.1,25);
g = -tan(phi).*(zdd + gravity);

I am not sure about gt. I chose it arbitrarily.

Then, I defined my ode function as follows:

function dydt = myode(t,gt,g)
g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t
dydt = g; % Evaluate ODE at time t

Then, in the main script again, I tried to integrate using ode45 as follows:

tspan = [0 1.99]; % time span
ic = 0; % initial condition
opts = odeset('RelTol',1e-2,'AbsTol',1e-4); % ode options
[t,y] = ode45(@(t) myode(t,gt,g), tspan, ic, opts);

However, if I try to run the script, I get an error saying that there are too many input arguments. Can someone tell me what I am doing wrong?

I understand that I can use cumtrapz() in my case to solve to do the numerical integration. However, my professor asked me to use ode45() since it offers a very accurate integration.

Moreover, I would like to point out that my main goal is to obtain . So, how can I use ode45() twice in this case?

If you are interested to find out how phi, zdd and gravity are defined, please refer to the code below:

I am solving an optimization problem using fmincon() to find the 9 parameters

z_hover1, z_start, z_end, z_hover2, phi_start, phi_end, t1, t2 ,t3

where t1,t2 and t3 are the times required to perform the trajectory in each of the three different phases of the flipping trajectory.

After finding the 9 parameters above, I can use them to create a flipping trajectory for a drone using a separate script. It should also be noted that I used the following function sevaral times in the script to come in order to find the trajectory of the position, velocity, acceleration, jerk and snap:

function traj = trajectory(start,goal,time)
 
    global step;
    [A,b] = build_matrix(start,goal,time);
    coeff_p = fliplr((A\b)');    
 
    t = step:step:time;
    % t = 0:step:time;
    
    coeff_v = polyder(coeff_p); % coefficients of the velocity polynomial
    coeff_a = polyder(coeff_v); % coefficients of the acceleration polynomial
    coeff_j = polyder(coeff_a); % coefficients of the jerk polynomial (derivative of the acceleration)
    coeff_s = polyder(coeff_j); % coefficients of the snap polynomial (derivative of the jerk)
    
    % Resulting trajectory
    traj = [polyval(coeff_p,t);
            polyval(coeff_v,t);
            polyval(coeff_a,t);
            polyval(coeff_j,t);
            polyval(coeff_s,t)];
    
end

So, the script which is then used to find the trajetory along y, z and phi is:

%% Reaching phase

z1_start = [z_hover1 0 0 0 0]; % initial condition of z, z_dot, z_ddot, z_dddot, z_ddddot in the first phase
z1_end = [z_start ((z_start-z_end)/t2 + g*t2/2) -g 0 0]; % final condition of z, z_dot, z_ddot, z_dddot, z_ddddot in the first phase
z1 = trajectory(z1_start,z1_end,t1);

phi1_start = [0 0 0 0 0]; % initial condition of phi, phi_dot, phi_ddot, phi_dddot, phi_ddddot in the first phase
phi1_end = [phi_start (phi_end-phi_start)/t2 0 0 0]; % Final condition of phi, phi_dot, phi_ddot, phi_dddot, phi_ddddot in the first phase
phi1 = trajectory(phi1_start,phi1_end,t1);

%% Flip phase

coeff_z2 = [-g/2 ((z_end-z_start)/t2+g*t2/2) z_start];
coeff_zd2 = polyder(coeff_z2);
coeff_zdd2 = polyder(coeff_zd2);
z2 = polyval(coeff_z2,step:step:t2); % trajectory of z in the second phase
zd2 = polyval(coeff_zd2,step:step:t2); % trajectory of zd in the second phase
zdd2 = polyval(coeff_zdd2,step:step:t2); % trajectory of zdd in the second phase

coeff_phi2 = [(phi_end-phi_start)/t2 phi_start];
coeff_phid2 = polyder(coeff_phi2);
coeff_phidd2 = polyder(coeff_phid2);
phi2 = polyval(coeff_phi2,step:step:t2); % trajectory of phi in the second phase
phid2 = polyval(coeff_phid2,step:step:t2); % trajectory of phid in the second phase
phidd2 = polyval(coeff_phidd2,step:step:t2); %  trajectory of phidd in the second phase


%% Recovery phase

z3_start = [z_end ((z_start-z_end)/t2-g*t2/2) -g 0 0]; % Initial condition of z, z_dot, z_ddot, z_dddot, z_ddddot in the third phase
z3_end = [z_hover2 0 0 0 0];  % Final condition of z, z_dot, z_ddot, z_dddot, z_ddddot in the third phase
z3 = trajectory(z3_start,z3_end,t3);

phi3_start = [phi_end (phi_end-phi_start)/t2 0 0 0]; % initial condition of phi, phi_dot, phi_ddot, phi_dddot, phi_ddddot in the third phase
phi3_end = [2*flips*pi 0 0 0 0]; % final condition of phi, phi_dot, phi_ddot, phi_dddot, phi_ddddot in the third phase
phi3 = trajectory(phi3_start,phi3_end,t3);

%% Global trajectory (concatenating all the three phases to obtain the global trajectory)

z = [z1(1,:) z2 z3(1,:)];
zd = [z1(2,:) zd2 z3(2,:)];
zdd = [z1(3,:) zdd2 z3(3,:)];

phi = [phi1(1,:) phi2 phi3(1,:)];
phid = [phi1(2,:) phid2 phi3(2,:)];
phidd = [phi1(3,:) phidd2 phi3(3,:)];

gravity = g*ones(size(z));

%% ydd must be integrated twice to find the trajectory along y
ydd = -tan(phi).*(zdd + gravity);

As you can see, all that is left is to integrate ydd twice in order to the trajectory along y. And, this is where I am facing my problem.

Thank you for taking the time to read my question.

Answer 1

[Elie Hatem]

I managed to make it work using by first checking the length of my phi, zdd and gravity vectors which were equal to 199.

Then, I after the first integration, I checked the length of T_yd which was 41. So I used this value to define ht for the second ode function which is just the output of the first ode45.

%% integrate ydd to obtain yd
total_time = t1+t2+t3
gt = linspace(0,total_time,199);
g = -tan(phi).*(zdd + gravity);


tspan = [0 total_time]; % time span
ic = 0; % initial condition
opts = odeset('RelTol',1e-2,'AbsTol',1e-4); % ode options
[T_yd,yd] = ode45(@(T_yd,yd) myode_ydd(T_yd,gt,g), tspan, ic, opts);


% plot(t,yd)

%% integrate yd to obtain y

ht = linspace(0,total_time,41);
h = yd;

tspan = [0 total_time]; % time span
ic = 0; % initial condition
opts = odeset('RelTol',1e-2,'AbsTol',1e-4); % ode options
[T_y,y] = ode45(@(T_y,y) myode_yd(T_y,ht,h), tspan, ic, opts);

figure
plot(T_y,y)
title('Postion y(t) with ode45'), grid

Answer 2

[James Tursa]

The function handle you pass to ode45( ) needs to have a (t,y) signature. E.g.,

[t,y] = ode45(@(t,y) myode(t,gt,g), tspan, ic, opts);

[Elie Hatem]

Thank you for your response. I tried to replace my line of code with the one above. However, I got the following erros:

Error using interp1>reshapeAndSortXandV (line 424)

X and V must be of the same length.

Error in interp1 (line 93)

[X,V,orig_size_v] = reshapeAndSortXandV(varargin{1},varargin{2});

Error in myode (line 2)

g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t

Error in tutorial>@(t,y)myode(t,gt,g) (line 9)

[t,y] = ode45(@(t,y) myode(t,gt,g), tspan, ic, opts);

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in tutorial (line 9)

[t,y] = ode45(@(t,y) myode(t,gt,g), tspan, ic, opts);

[James Tursa]

You don't show us how phi and zdd and gravity are generated, so we have no way of checking your code to see if they are compatible with the size of gt.

[Elie Hatem]

I updated the post if you want to have a look at the remainder of the code. Thank you.

[Elie Hatem]

Thank you for your comment. I managed to make it work thanks to you. I answered my question again with the solution.