书上课后例题,导弹飞行问题的a,c两题。
clc;clear;close all;
% 定义初始条件
c = 0.002; % 拉力系数
g = 32.2; % 重力加速度
v0 = 600; % 初始速度
a = pi/4; % 初始角度
v0x = v0*cos(a); % 水平分量
v0y = v0*sin(a); % 垂直分量
y0 = 0; % 初始高度
t0 = 0; % 初始时间
yinit = [v0x;v0y;y0;t0]; % 初始条件向量
% 定义求解区间
x0 = 0; % 初始位置
xfinal = 1000; % 最终位置
% 调用ode45函数求解
[x,y] = ode45(@missile,[x0,xfinal],yinit);
% 第一题:画出导弹飞行轨迹
ind = find(y(:,3)>0); % 找到y>0的索引
plot(x(ind),y(ind,3),'b-') % 画出x-y曲线
xlabel('水平距离 x (ft)') % x轴标签
ylabel('垂直高度 y (ft)') % y轴标签
title('导弹飞行轨迹') % 标题
grid on % 网格线
% 第三题:求y=0时的x值及飞行时间
x0 = interp1(y(:,3),x,0); % 用插值法求y=0时的x值
t0 = interp1(y(:,3),y(:,4),0); % 用插值法求y=0时的t值
fprintf('导弹在 y=0 时的水平距离为 %.2f ft,飞行时间为 %.2f s\n',x0,t0)
% 定义导弹飞行轨迹的四个一阶方程
function dydx = missile(x,y)
% 参数设置(所有单位均为英尺)
c = 0.002; % 拉力系数
g = 32.2; % 重力加速度
v0 = 600; % 初始速度
a = pi/4; % 初始角度
% 条件检查
if abs(y(1)) < v0*10^(-6) % 如果水平速度太小,终止程序
error('水平速度太小,无法继续计算')
end
% 四个一阶方程
dydx = zeros(4,1); % 初始化输出向量
dydx(1) = -c*y(1); % dvx/dx=-cv
dydx(2) = -(g+c*y(2)*v0/y(1)); % dvy/dx=-(g+cv*vy)/vx
dydx(3) = y(2)/y(1); % dy/dx=vy/vx
dydx(4) = 1/y(1); % dt/dx=1/vx
end