MATLAB Answers

Error in ode45 while doing numerical integration

3 views (last 30 days)
Elie Hatem
Elie Hatem on 1 Jun 2021
Commented: Elie Hatem on 2 Jun 2021
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.

Accepted Answer

Elie Hatem
Elie Hatem on 2 Jun 2021
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

More Answers (1)

James Tursa
James Tursa on 1 Jun 2021
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);
  4 Comments
Elie Hatem
Elie Hatem on 2 Jun 2021
Thank you for your comment. I managed to make it work thanks to you. I answered my question again with the solution.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!