How to use ODE 45 to integrate equations of motion?

13 次查看（过去 30 天）
Stephen Mixon
Stephen Mixon 2019-10-28
评论： Cyrus Tirband 2019-10-29
I am modeling a 3d projectile and need help using the ode45 command. I need it to integrate the equations of motion based on my simulation and for it to end when z=0. I need to take this numerical results and sumperimpose them in a plot with my previous data. I need to plot the errors in my data beween my parameter set and this ode 45 set.
This is my code for the simulation
figure('name','Projectile Motion')
xlabel('X (m)')
ylabel('Y (m)')
zlabel('Z (m)')
% This sets the viewing angle
view([1 -1 1])
% This sets the limits on the x- and y-axes
xlim([0,3000])
ylim([0,3000])
zlim([0,1000])
box('on')
grid('on')
hold on
% az is azimuth angle, el is elevation
V0 = 150;
az = 45;
el= 30;
t = [0:0.01:25]';
vx0 = V0*cosd(az)*cosd(el);
vy0 = V0*sind(az)*cosd(el);
vz0 = V0*sind(el);
x = vx0.*t;
y = vy0.*t;
z = vz0*t-0.5*9.81.*t.^2;
[~,ind]=min(abs(z));
plot3(x,y,z)
%Time of flight
T =(2*vz0)/9.81
  7 个评论
显示 5更早的评论
Stephen Mixon
Stephen Mixon 2019-10-28
Write MATLAB code that uses the ode45 package to numerically integrate the equations of
motion and simulate the flight of the projectile. The integration should stop when the projectile
reaches . The simulation and subsequent post-processing should give you all the
z=0
parameters that you obtained analytically in Part 1a.
Compare the numerical results to the analytical results for each parameter computed in Part 1a
by superimposing plots. Time should be on the x-axis and the parameter on the y-axis. Use
dashed lines for analytical results and a continuous line of a different color for the numerical
results.
darova
darova 2019-10-28
What should be the result of integration? What is the parameter? t?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Cyrus Tirband
Cyrus Tirband 2019-10-28
First, you need to write down your equations of motion. In your simple case, you have:
with initial conditions
Of course, these integrals are fairly easy to solve and you don't really need to verify these with a numerical solver, but for your assignment you are asked to compare them.
Knowing the equations, the ode function becomes easy to set up. Each of them are independent, and almost the same, so
function dzdt = odefunc(t,z,g)
dzdt = [y(2), -g];
end
similar for x, y (or set g to zero)
then call the solver
[t,z] = ode45(@(t,z), odefunc(t,z,g), [0 tend],[0; vz0]);
The solutions can be seen to overlap perfectly.
ode_example.png
  2 个评论
Stephen Mixon
Stephen Mixon 2019-10-28
I used the code you gave me and inseted inital conditions to 0 and got this error
Function definitions in a script must appear at the end of the file.
Move all statements after the "odefunc" function definition to before the first local
function definition.
here is the updated code
figure('name','Projectile Motion')
xlabel('X (m)')
ylabel('Y (m)')
zlabel('Z (m)')
% This sets the viewing angle
view([1 -1 1])
% This sets the limits on the x- and y-axes
xlim([0,3000])
ylim([0,3000])
zlim([0,1000])
box('on')
grid('on')
hold on
% az is azimuth angle, el is elevation
V0 = 150;
az = 45;
el= 30;
t = [0:0.01:25]';
vx0 = V0*cosd(az)*cosd(el);
vy0 = V0*sind(az)*cosd(el);
vz0 = V0*sind(el);
x = vx0.*t;
y = vy0.*t;
z = vz0*t-0.5*9.81.*t.^2;
[~,ind]=min(abs(z));
plot3(x,y,z)
%Time of flight
T =(2*vz0)/9.81;
x_a=0;
y_a=0;
z_a=0;
x_v=vx0;
x_v=vy0;
z_v=vz0;
x0=0;
y0=0;
z0=0;
function dzdt = odefunc(t,z,g)
dzdt = [y(2), -g];
end
[t,z] = ode45(@(t,z), odefunc(t,z,g), [0 tend],[0; vz0]);
Cyrus Tirband
Cyrus Tirband 2019-10-29
The error message is fairly unambiguous. You cannot define functions in the middle of the script, you have to put them at the end.
If you're using MATLAB 2016a or earlier, you'll have to define them in their own file, or make the main script a function.
Additionally, you'll have to define a few parameters before my snippets will generate the solution.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by