【数学建模】《实战数学建模：例题与讲解》第五讲-微分方程建模（含Matlab代码）

如果这篇文章对你有帮助，欢迎点赞与收藏~

基本概念

微分方程建模是数学建模中一种极其重要的方法，它在解决众多实际问题时发挥着关键作用。这些实际问题的数学表述通常会导致求解特定的微分方程。将各种实际问题转换为微分方程的定解问题主要包括以下几个步骤：

1. 确定研究对象：首先，需要明确研究的量，这包括自变量、未知函数和必要的参数等，并选择合适的坐标系。

2. 找出基本规律：接下来，需要识别这些量所遵循的基本规律，这些规律可能源自物理、几何、化学或生物学等领域。

3. 列出方程和定解条件：利用这些规律来构建方程和相应的定解条件。

在列出方程时，常见的方法有以下几种：

1. 直接列方程法：在数学、力学、物理、化学等领域，许多自然现象的规律已被广泛认识，并可直接用微分方程描述。如牛顿第二定律、放射性物质的衰变规律等。常利用这些已知规律直接建立微分方程。

2. 微元分析法和积分法：对于一些自然现象，其规律需通过变量的微元间的关系来表达。在这种情况下，不能直接列出自变量与未知函数及其变化率之间的关系，而是通过微元分析法，利用已知规律建立变量间的微元关系，再通过取极限得到微分方程，或通过在任意区域上进行积分来建立微分方程。

3. 模拟近似法：在生物学、经济学等领域，许多现象的规律不明确且相当复杂。在这些情况下，需要根据实际资料或实验数据提出各种假设，然后在这些假设下，利用合适的数学方法列出微分方程。

在实际的微分方程建模过程中，这些方法往往是综合应用的。无论采用哪种方法，通常需要根据实际情况做出一定的假设和简化，并将模型的理论或计算结果与实际情况对比验证，以此来修改模型，使其更准确地描述实际问题，从而达到预测和预报的目的。

习题6.4

1. 题目要求

2.解题过程

解：

(1)

对函数值 $s(t)$ 的影响因素，根据题目分析，可以得到如下两点：

• 下降速度与 $s(t)$ 的值的大小成正比，设比例系数为 $\lambda$ ，所以式子中有一个部分为：$-\lambda s(t)$ 。
• 增长速度与 $a(t)$ 的值的大小成正比，设比例系数为 $\mu$ ，所以式子中有一个部分为：$+\mu a(t)$ 。

但是，题目中说到，广告只能影响这种商品在市场上尚未饱和的部分(设饱和量为M)，所以有两个选择：分段函数 或者 渐进函数。

我选择使用渐进函数，渐近线为M，增长部分改为：$\mu a(t)\left(1-\frac{s(t)}{M}\right)$ 。

这样，当 s=M 或者 a=0 时，销量都不会上涨，符合题意。

将所有分析合并，建立如下模型：
$$\frac{\mathrm{d}s}{\mathrm{d}t}=\mu a(t)\left(1-\frac{s(t)}{M}\right)-\lambda s(t)$$

(2)

广告宣传只进行有限时间 $\tau$ ，且广告费为常数a，所以建立下列函数
$$a(t)=\{\begin{array}{ll}a/c,& 0<t⩽\tau ,\\ 0,& t⩾\tau .\end{array}$$
带入(1)中公式，并移项，有：
$$\frac{\mathrm{d}s}{\mathrm{d}t}+\left(\lambda +\frac{\mu a}{M\tau }\right)s=\frac{\mu a}{\tau },0<t<\tau ,$$

(3)

先讨论 $0<t⩽\tau$ :

为了运算方便，进行换元：
$$b=\lambda +\frac{\mu a}{M\tau },c=\frac{\mu a}{\tau }$$
则上述式子化简为：
$$\frac{\mathrm{d}s}{\mathrm{d}t}+bs=c,$$

(4)

这是一个一阶非齐次线性微分方程.

假设初始值为 $s_0$ ,可以得出结果为:
$$s(t)=\frac{c}{b}\left(1-{\mathrm{e}}^{-bt}\right)+s_0{\mathrm{e}}^{-bt},0<t⩽\tau$$

(5)

再讨论 $t⩾\tau$ :

式子为：
$$\frac{\mathrm{d}s}{\mathrm{d}t}=-\lambda s,$$
可以得出结果为：
$$s(t)=s(\tau ){\mathrm{e}}^{\lambda (\tau -t)},t⩾\tau .$$
(6)

合并所有结果，得到：
$$s(t)=\{\begin{array}{ll}\frac{c}{b}\left(1-{\mathrm{e}}^{-bt}\right)+s_0{\mathrm{e}}^{-bt},& 0<t⩽\tau ,\\ s(\tau ){\mathrm{e}}^{\lambda (\tau -t)},& t⩾\tau .\end{array}$$

3.结果

$s(t)$ 的变化符合下列函数：
$$s(t)=\{\begin{array}{ll}\frac{c}{b}\left(1-{\mathrm{e}}^{-bt}\right)+s_0{\mathrm{e}}^{-bt},& 0<t⩽\tau ,\\ s(\tau ){\mathrm{e}}^{\lambda (\tau -t)},& t⩾\tau .\end{array}$$
其中，
$$b=\lambda +\frac{\mu a}{M\tau },c=\frac{\mu a}{\tau }$$

习题6.5(1)

1. 题目要求

2.解题过程

解：

Matlab的工具箱提供了几个解常微分方程的功能函数，如ode45 、ode23 、ode113 。其中，ode45 采用四五阶龙格库塔方法（以下简称RK方法），是解非刚性常微分方程的首选方法； ode23 采用二三阶龙格库塔方法。

高阶常微分方程，必须做变量替换， 化成一阶微分方程组，才能使用Matlab求数值解。

设 $y_1=y,y_2=y^{\prime }$ ,则原式化为：
$$\{\begin{array}{ll}{y}_1^{\prime }=y_2,& y_1\left(\frac{\pi }{2}\right)=2,\\ y_2^{\prime }=\left(\frac{n^2}{x^2}-1\right)y_1-\frac{y_2}{x},& y_2\left(\frac{\pi }{2}\right)=-\frac{2}{\pi }.\end{array}$$

接下来就可以用ode45来解决相关问题。

3.程序

求解的MATLAB程序如下：

clc, clear

syms y(x) % 符号变量
dy = diff(y); % 一阶导
d2y = diff(dy); % 二阶导

% 为了分析比较，首先求解符号解
y = dsolve(x^2*d2y+x*dy+(x^2 - 1 / 4)*y, ...
    y(pi/2) == 2, dy(pi/2) == -2/pi); % 直接所有内容
y = simplify(y); % 化简结果
symdisp(y); % 直观输出函数（自建函数）
fplot(y); % 画图
hold on

%龙格库塔方法求解数值解
dy_2 = @(x, y) [y(2); (1 / 4 / x^x - 1) * y(1) - y(2) / x]; % 方程组的右端 注意：用分号隔开！！！

% 求解的区间是 [pi / 2, 15]
% 注意：2和-2/pi是区间左端点的时候的函数值，即x=pi/2时
[x, y] = ode45(dy_2, [pi / 2, 15], [2, -2 / pi]);
plot(x, y(:,1), '*'); % 画图 注意画的是y1

legend('符号解','数值解');

4.结果

微分方程求解结果如下：

通过函数曲线可以看出，解析解和数值解吻合。

习题6.5(2)

1. 题目要求

2.解题过程

解：

设 $y_1=y,y_2=y^{\prime }$ ,则原式化为：
$$\{\begin{array}{ll}{y}_1^{\prime }=y_2,& y_1(0)=1,\\ y_2^{\prime }=-y_1\cos x,& y_2(0)=0.\end{array}$$
接下来就可以用ode45来解决相关问题。

3.程序

求解的MATLAB程序如下：

clc, clear

% 使用幂级数解作为对照
power_series = @(x) 1 - 1 / gamma(3) * x.^2 + 2 / gamma(5) * x.^4 - 9 / gamma(7) * x.^6 + 55 / gamma(9) * x.^8;
x = 0:0.1:4;
y = power_series(x);
plot(x, y);
hold on

%龙格库塔方法求解数值解
dy = @(x, y) [y(2); -y(1) * cos(x)]; % 方程组的右端 注意：用分号隔开！！！

[x, y] = ode45(dy, [0, 4], [1, 0]);
plot(x, y(:, 1), '*'); % 画图 注意画的是y1

legend('幂级数解', '数值解');

4.结果

通过函数曲线可以看出，当x较小的时候，级数解的近似值和数值解吻合较好。

习题6.6

1. 题目要求

2.第(1)问

解：

(1)

如图建立坐标系。

我们可以看出，小船的和速度是水速和静水速度的向量和，其中，水速平行与河岸（y轴），静水速度始终朝向B点。

根据正交分解，我们可以求出x与y方向的速度：
$$\frac{\mathrm{d}x}{\mathrm{d}t}=-v_2\cos \theta ,\frac{\mathrm{d}y}{\mathrm{d}t}=v_1-v_2\sin \theta .$$

根据几何关系，我们可以求出：
$$\cos \theta =\frac{x}{\sqrt{x^2+y^2}},\sin \theta =\frac{y}{\sqrt{x^2+y^2}},$$
联立（14）（15）式，得到：
$$\frac{\mathrm{d}x}{\mathrm{d}t}=-v_2\frac{x}{\sqrt{x^2+y^2}},\frac{\mathrm{d}y}{\mathrm{d}t}=v_1-v_2\frac{y}{\sqrt{x^2+y^2}},$$
将（16）的两个式子相除，得到：
$$\frac{\mathrm{d}y}{\mathrm{d}x}=-k\frac{\sqrt{x^2+y^2}}{x}+\frac{y}{x},0<x<d$$
这是一个一阶齐次微分方程，于是得到初值问题：
$$\{\begin{array}{l}\frac{\mathrm{d}y}{\mathrm{d}x}=-k\frac{\sqrt{x^2+y^2}}{x}+\frac{y}{x},0<x<d\\ y(d)=0\end{array}$$
为了求解这个微分方程，我试图使用Matlab的dsolve函数解决，但是matlab无法解决这个微分方程。

使用Wolfram Mathematica软件的Dsolve函数，成功解出了微分方程的解为：
$$y=\frac{d}{2}[{\left(\frac{x}{d}\right)}^{1-k}-{\left(\frac{x}{d}\right)}^{1+k}],0<x<d.$$

3.第(2)问

根据第一问，可知参数方程为：
$$\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=-\frac{2x}{\sqrt{x^2+y^2}},x(0)=d\\ \frac{\mathrm{d}y}{\mathrm{d}t}=1-\frac{2y}{\sqrt{x^2+y^2}},y(0)=0\end{array}$$
求解的MATLAB程序如下：

clc, clear

d = 100;
v1 = 1;
v2 = 2;
k = 0.5;

% 第一问的解析解
y = @(x)d / 2 * ((x / d).^(1 - k) - (x / d).^(1 + k));
fplot(y, [0, 100]);
hold on

%龙格库塔方法求解数值解
dy = @(t, xy) [-2 * xy(1) / sqrt(xy(1)^2+xy(2)^2); 1 - 2 * xy(2) / sqrt(xy(1)^2+xy(2)^2)]; % 方程组的右端 注意：用分号隔开！！！

[t, xy] = ode45(dy, [0, 60], [100, 0]);
plot(xy(:, 1), xy(:, 2), '*'); % 画图 注意画的是y1

legend('符号解', '数值解');

4.结果

渡河所需时间需要不断地尝试，修改处如下图所示：

测试过程如下图（测试数据为40，很明显未到岸）：

最终，试验得到渡河所用时间为66.6秒，数值解与解析解的对照图（航行曲线）如下图：

如果这篇文章对你有帮助，欢迎点赞与收藏~