MATLAB Answers

Error in ode45 while doing numerical integration

조회 수: 2(최근 30일)
Elie Hatem
Elie Hatem 2021년 6월 1일
댓글: Elie Hatem 2021년 6월 2일
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.

채택된 답변

Elie Hatem
Elie Hatem 2021년 6월 2일
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

추가 답변(1개)

James Tursa
James Tursa 2021년 6월 1일
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
Elie Hatem
Elie Hatem 2021년 6월 2일
Thank you for your comment. I managed to make it work thanks to you. I answered my question again with the solution.

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by