数学实验舰艇追击问题

一个半径为0.5米的均匀木球，密度为0.65，在离水面0.8米处自由下落（球心和水面的距离）。经历如下的三个阶段，第一，空气中的自由下落，第二，部分侵入水中的下沉，第三，全部侵入水中的下沉。为简单起见，仅仅考虑木球在水中受到浮力而不考虑其他的复杂因素。试确定三个阶段经历的时间和运动的距离（数学实验，用MATLAB解决微分方式解决问题）

为了求解这个问题，可以采用数学建模的方法，分别考虑球体在空气中和在水中的运动情况，通过求解微分方程得到球体的下落时间和下落距离。

首先，在空气中自由下落的情况下，球体受到的外力只有重力，因此可以根据牛顿第二定律得到如下微分方程：

m * d^2h/dt^2 = -mg

其中m为球体的质量，g为重力加速度，h为球心距离水面的距离。将密度和半径代入可以求得球体质量：

m = 0.65 * (4/3) * pi * (0.5)^3

代入上述微分方程可以得到：

d^2h/dt^2 = -g

这是一个简单的二阶微分方程，可以通过MATLAB的ode45函数求解。定义状态变量y为[h, dh/dt]，则可以得到如下代码：

function dydt = free_fall(t, y)
g = 9.81; % 重力加速度
dydt = [y(2); -g];
end

tspan = [0, 10]; % 求解时间范围
y0 = [0.8, 0]; % 初始状态
[t, y] = ode45(@free_fall, tspan, y0);

运行后，可以得到球体在空气中自由下落的时间和下落距离：

>> t(end)
ans =
    4.0149

>> y(end, 1)
ans =
   -4.5504

下一步，考虑球体部分侵入水中的情况。此时，球体除了受到重力之外，还要考虑浮力的影响。根据阿基米德定律，球体受到的浮力大小等于球体排开的水的重量，即

F_buoyancy = rho_water * V_submerged * g

其中rho_water为水的密度，V_submerged为球体在水中的体积。根据球体在水中的情况，可以得到球体在水中的体积为：

V_submerged = (4/3) * pi * (0.5)^2 * (0.8 - h)

将这个式子代入浮力公式中，可以得到球体受到的合力为：

F_net = F_buoyancy - mg = rho_water * (4/3) * pi * (0.5)^2 * (0.8 - h) * g - 0.65 * (4/3) * pi * (0.5)^3 * g

根据牛顿第二定律，可以得到如下微分方程：
m * d^2h/dt^2 = rho_water * (4/3) * pi * (0.5)^3 * g - rho_air * (4/3) * pi * (0.5)^3 * g

其中，m为木球的质量，rho_water为水的密度，rho_air为空气的密度，g为重力加速度，h为球心离水面的距离。

该微分方程是一个二阶非齐次线性微分方程，可以通过常系数线性微分方程的解法求解，具体方法如下：

设h的通解为h(t) = h0 + Acos(wt) + Bsin(wt) + h1t + h2，其中h0为初始高度，A、B、h1、h2为待定常数，w为角频率。

将该通解带入微分方程，可以得到：

-mw^2(h0 + Acos(wt) + Bsin(wt) + h1t + h2) = rho_water * (4/3) * pi * (0.5)^2 * (0.8 - h0 - Acos(wt) - Bsin(wt) - h1t - h2) * g - rho_air * (4/3) * pi * (0.5)^2 * (h0 + Acos(wt) + Bsin(wt) + h1t + h2) * g

整理后，可以得到如下关于A、B、h1、h2的代数方程组：

(-mw^2 + rho_water * (4/3) * pi * (0.5)^2 * g) * A + rho_water * (4/3) * pi * (0.5)^2 * g * B - rho_water * (4/3) * pi * (0.5)^2 * g * h1 = 0
(-mw^2 + rho_water * (4/3) * pi * (0.5)^2 * g) * B - rho_water * (4/3) * pi * (0.5)^2 * g * A - m * h1 = 0
(-mw^2 + rho_water * (4/3) * pi * (0.5)^2 * g) * h1 - rho_water * (4/3) * pi * (0.5)^2 * g * h2 = m * g

解出A、B、h1、h2的值后，将其代入通解，就可以得到h(t)的具体表达式。

对于第一阶段的自由下落，初始高度为0.8米，初始速度为0米/秒，根据一阶微分方程可以得到：

h(t) = 0.8 - 0.5 * g * t^2

其中g为重力加速度。
对于第二阶段和第三阶段，可以根据上述通解得到：

第二阶段：球部分侵入水中的下沉，从$h_2$开始到$h_3$结束，此时有$0\leq h_2\leq 0.8$且$h_3\leq h_2$。由于水的浮力作用，球所受合外力变为$F=mg-\rho_{water}Vg$，此时的微分方程变为：
�
�
2
ℎ
�
�
2
=
�
�
−
�
�
�
�
�
�
�
�
m
dt
2

d
2
h
​
=mg−ρ
water
​
Vg

其中$V$为球在水中的体积，即$\frac{4}{3}\pi(0.5)^3-\frac{4}{3}\pi h^3$，则上述微分方程可以写为：

�
�
2
ℎ
�
�
2
=
�
�
−
�
�
�
�
�
�
(
4
3
�
(
0.5
)
3
−
4
3
�
ℎ
3
)
�
m
dt
2

d
2
h
​
=mg−ρ
water
​
(
3
4
​
π(0.5)
3
−
3
4
​
πh
3
)g

化简得到：

�
2
ℎ
�
�
2
=
�
2
−
4
3
�
�
�
�
�
�
�
(
1
2
)
2
ℎ
3
dt
2

d
2
h
​
=
2
g
​
−
3
4
​

m
ρ
water
​

​
(
2
1
​
)
2
h
3

这是一个非线性微分方程，无法用解析法解决，需要使用数值方法求解。

第三阶段：球全部侵入水中的下沉，从$h_3$开始到$h=0$结束，此时有$0\leq h_3\leq 0.8$。同样由于水的浮力作用，球所受合外力变为$F=mg-\rho_{water}Vg$，微分方程与第二阶段相同：
�
2
ℎ
�
�
2
=
�
2
−
4
3
�
�
�
�
�
�
�
(
1
2
)
2
ℎ
3
dt
2

d
2
h
​
=
2
g
​
−
3
4
​

m
ρ
water
​

​
(
2
1
​
)
2
h
3

需要使用数值方法求解。

不知道你这个问题是否已经解决, 如果还没有解决的话:

如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^