matlab编程求加速度积分

Question

d2x/dt2=0.435(U-dx/dt)[(U-dx/dt)^2+(dz/dt)^2]^(1/2)加速度方程

d2z/dt2=1.1856(-dz/dt)[(U-dx/dt)^2+(dz/dt)^2]^(1/2)-9.8加速度方程
U=20(z/10)^0.12
根据方程积分求运动轨迹，以及x、y方向的速度
t=0时x=0，z=500，v=0，时间间隔t=0.1s，求到z=0停止
谢谢，不胜感激

Answer 1

这个问题我之前回答过（编号2117613008218288267），当时因为题主给的数据有问题没能彻底解决。现在给的条件基本上一个没问题了，给你新版的程序：

function main
X0 = [0;0;500;0];
[t,X] = ode45(@ode,0:0.1:1000,X0,odeset('Events',@events)); 
ystr = {'x', 'dx/dt', 'z', 'dz/dt'}; 
for i=1:4 
    subplot(2,2,i) 
    plot(t(1:end-1),X(1:end-1,i)); 
    xlabel('{\itt} / s') 
    ylabel(ystr{i}) 
end 
function dX = ode(t, X)
% X(1)=x, X(2)=dx/dt, X(3)=z, X(4)=dz/dt
% x = X(1);
dx = X(2);
z = X(3);
dz = X(4);
U = 20*(z/10)^0.12;
dX = [ ... 
    dx; ...
    0.435*(U-dx)*((U-dx)^2+dz^2)^(1/2); ... 
    dz; 
    1.1856*(-dz)*((U-dx)^2+dz^2)^(1/2) - 9.8 ... 
    ];
function [value,isterminal,direction] = events(t,X)
value = X(3);
isterminal = 1;
direction = 0;

和之前的差别，除了模型按照现在新的模型之外，最主要是增加了事件函数，以实现z到达0就能够结束。另外，微分方程的写法之前用匿名函数，现在改用具名函数，更清晰一些。仿真的最后一个点可能出现复数，所以绘图时删除最后一个点。