Question

一个简单的matlab作业，请大家帮下忙。

求代码，最好带解释的，谢谢！
物体由高度h处以速度v0水平抛出，如图11-21所示（pp.257）。空气阻力F大小可视为与速度v大小的n次方成正比，即F = kmvn ，其中m为物体的质量，k为常系数。阻力F方向与速度v方向相反。求物体的运动方程和轨迹。（h，v0，k，m的值可自定）
1）分别画出n=1,2,3时的物体运动轨迹图，并给出落地所需时间。
2）画出n=1时物体飞行的动画。有动画，有真相。
3）作业考核采用面试方式，通过即可得分。独立完成。

Answer 1

clear all;

close all;

h = 10;

v0 = 50;

k = 0.02;

m = 5;

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);

[t1,y1] = ode45(@MFy1,[0,1.440],[h,0,0,v0]);

% y第一列是y向位移，第二列是y向速度，第三列是x向位移，第四列是x向速度

[t2,y2] = ode45(@MFy2,[0,1.477],[h,0,0,v0]);

[t3,y3] = ode45(@MFy3,[0,1.765],[h,0,0,v0]);

plot(y1(:,3),y1(:,1),'r',y2(:,3),y2(:,1),'g',y3(:,3),y3(:,1),'b')

grid on;

legend('n=1','n=2','n=3');

xlabel('x向位移');

ylabel('y向位移');

给出函数MFy1，类似的MFy2，MFy3自己写一下。

function dy = MFy1(t,y)

dy = zeros(4,1);

dy(1) = y(2);%y向速度

dy(2) = -9.8 + abs(y(2))^1*0.02;%y向加速度

dy(3) = y(4);%x向速度

dy(4) = -abs(y(4))^1*0.02;%x向加速度

动画嘛稍后给出；

这里值得一提的是落地时间是估算出来的，如何设置计算终止条件还有待完善。这种解微分方程问题如果放到simlink里会容易很多，但是由于长时间未用simlink我一时也不能弄出来，你可以自己试试。

动画的函数，k是演示的次数

function f = anim_down(k)

h = 10;

v0 = 50;

ht = h;

tend = 1.4;

while ht >= 1e-5

    [t1,y1] = ode45(@MFy1,[0,tend],[h,0,0,v0]);

    ht = y1(end,1);

    tend = tend + 0.005;

end

t2 = t1([2:end])-t1([1:end-1]);

t2 = [0;t2];

plot(y1(:,3),y1(:,1),'Color',[1,1,1]);

grid on;

xlabel('x向位移');

ylabel('y向位移');

hold on;

h = line('Color',[1 0 0],'marker','.','markersize',24,'erasemode','xor');

n = length(t1);i = 1;j = 1;

while 1

    set(h,'xdata',y1(i,3),'ydata',y1(i,1));

    if i > 1

        plot(y1(i-1:i,3),y1(i-1:i,1),'Color',[0,0,1]);

    end

    drawnow;

    pause(t2(i));

    i = i+1;

    if i > n

        i = 1;

        j = j+1;

        if j>k

            break;

        end

    end

end