研究受空气阻尼的抛体运动。其中:抛体质量为m.初速度为v0，所受空气阻力R的大小与速率v的n次方成正比，，其中b是阻尼系数。程序中取m=1，空气阻尼为[0. 0.2，0.2]，即分别为无阻尼，阻力正比于v，和阻力正比于v^2。请写出用 ode45 指令求解的该问题的 matiab 程序。 要求:1.写出问题的坐标系，质点的微分方程和改写后的一阶微分方程组; 2.面图显示最后的 y-x，y-t，dx/dt-t，和 dx/dt-t 的图形

1. 坐标系

我们可以使用笛卡尔坐标系，将抛体的位置表示为(x,y)，速度表示为(dx/dt,dy/dt)。

2. 微分方程

考虑受空气阻尼的抛体运动，我们可以将运动方程表示为一个微分方程组：

m*d^2x/dt^2 = -b*(dx/dt)^n
m*d^2y/dt^2 = -b*(dy/dt)^n - mg

其中，m是抛体的质量，b是空气阻力系数，g是重力加速度，n表示空气阻力与速度的n次方成正比。

我们将该微分方程组改写成一阶微分方程组的形式，即：

dx/dt = v_x
dv_x/dt = -b/m*(v_x)^n

dy/dt = v_y
dv_y/dt = -b/m*(v_y)^n - g

其中，v_x和v_y分别为x和y方向上的速度。

3. Matlab程序

基于以上微分方程组，我们可以编写如下的Matlab程序来求解受空气阻尼的抛体运动：

function [t,x,y,vx,vy] = projectile_motion(m,b,g,v0,theta,n)
% 定义运动方程
f = @(t,x) [x(2); -b/m*x(2)^n; x(4); -b/m*x(4)^n - g];

% 定义初值条件
x0 = [0; v0*cos(theta); 0; v0*sin(theta)];

% 求解微分方程组
[t,x] = ode45(f, [0, 10], x0);

% 计算位置和速度
y = x(:,3);
x = x(:,1);
vy = x(:,4);
vx = x(:,2);

% 绘制图形
subplot(2,2,1);
plot(x,y);
xlabel('x (m)');
ylabel('y (m)');
title('Projectile motion with air resistance');

subplot(2,2,2);
plot(t,y);
xlabel('t (s)');
ylabel('y (m)');
title('y-t graph');

subplot(2,2,3);
plot(t,vx);
xlabel('t (s)');
ylabel('dx/dt (m/s)');
title('dx/dt-t graph');

subplot(2,2,4);
plot(t,vy);
xlabel('t (s)');
ylabel('dy/dt (m/s)');
title('dy/dt-t graph');
end

其中，m、b、g分别为抛体的质量、空气阻力系数和重力加速度；v0和theta分别为抛体的初速度和发射角度；n表示空气阻力与速度的n次方成正比。

我们将运动方程写成一个匿名函数f，并定义初值条件x0，然后使用ode45函数求解微分方程组。最后，我们计算抛体的位置和速度随时间的变化情况，并绘制出四个图形，分别为最终的y-x图、y-t图、dx/dt-t图和dy/dt-t图。

例如，我们可以在Matlab命令行中调用函数，例如：

[t,x,y,vx,vy] = projectile_motion(1, 0.1, 9.8, 20, pi/4, 2);

这里我们假设抛体的质量为1kg，空气阻力系数为0.1，重力加速度为9.8m/s^2，初速度为20m/s，发射角度为45度，空气阻力与速度的平方成正比。运行后，Matlab会返回时间和位置随时间的变化情况，并绘制出四个图形。