solveodeusingmatlab.txt

《大学数学实验》中一道解微分方程的作业题

function dx=huojian(t,x)
% x(1)是高度, x(2)是速度
f1=x(2);
f2=(32000-(1400-18*t)*9.8-0.4*x(2)^2)/(1400-18*t);
dx=[f1;f2];


function dx=huojian1(t,x)
f1=x(2);
f2=-(320*9.8+0.4*x(2)^2)/320;
dx=[f1;f2];


ts=0:0.1:60;
x0=[0,0];
[t,x]=ode23('huojian',ts,x0); % 计算高度函数和速度函数
h=x(601,1)                    % 引擎关闭瞬间火箭的高度
v1=x(601,:);                 
va=huojian(60,v1);
v=va(1)                        % 引擎关闭瞬间火箭的速度
a=va(2)                       % 引擎关闭瞬间火箭的加速度                      

ts1=61:0.1:100;
x01=[0,v];
[t1,x1]=ode23('huojian1',ts1,x01);     % 计算引擎关闭后的高度函数和速度函数
m=min(find(x1(:,2)<=0));               % 计算达到最高点时的时间
h1=x1(m,1);                          
h1=h+h1                                % 火箭达到最高点时的高度
v2=x1(m,:);
va1=huojian1(m,v2);
v1=va1(1)                             % 火箭达到最高点时的速度, 准确值为0        
a1=va1(2)                             % 火箭达到最高点时的加速度, 准确值为-9.8


ts2=[ts ts1(1,1:m)];
hh1=x1(1:m)+h;

h2=[x(:,1)' hh1];
subplot(2,2,1)
plot(ts2,h2)
axis([20 80 0 15000])                % 画出高度随时间变化的图像 

for k=1:length(ts)
    va2(:,k)=huojian(t(k),x(k,:));    % 计算引擎关闭前各时刻的速度和加速度
end
for k=1:m
    va3(:,k)=huojian1(t1(k),x1(k,:));  % 计算引擎关闭后各时刻的速度和加速度
end
va4=[va2 va3];

subplot(2,2,2)
plot(ts2,va4(1,:))                       % 画出速度随时间变化的图像 
axis([0 80 0 300])

subplot(2,2,3)
plot(ts2,va4(2,:))
axis([0 80 -100 20])                      % 画出加速度随时间变化的图像