如果这篇文章对你有帮助，欢迎点赞与收藏~
基本概念
微分方程建模是数学建模中一种极其重要的方法，它在解决众多实际问题时发挥着关键作用。这些实际问题的数学表述通常会导致求解特定的微分方程。将各种实际问题转换为微分方程的定解问题主要包括以下几个步骤：
-
确定研究对象：首先，需要明确研究的量，这包括自变量、未知函数和必要的参数等，并选择合适的坐标系。
-
找出基本规律：接下来，需要识别这些量所遵循的基本规律，这些规律可能源自物理、几何、化学或生物学等领域。
-
列出方程和定解条件：利用这些规律来构建方程和相应的定解条件。
在列出方程时，常见的方法有以下几种：
-
直接列方程法：在数学、力学、物理、化学等领域，许多自然现象的规律已被广泛认识，并可直接用微分方程描述。如牛顿第二定律、放射性物质的衰变规律等。常利用这些已知规律直接建立微分方程。
-
微元分析法和积分法：对于一些自然现象，其规律需通过变量的微元间的关系来表达。在这种情况下，不能直接列出自变量与未知函数及其变化率之间的关系，而是通过微元分析法，利用已知规律建立变量间的微元关系，再通过取极限得到微分方程，或通过在任意区域上进行积分来建立微分方程。
-
模拟近似法：在生物学、经济学等领域，许多现象的规律不明确且相当复杂。在这些情况下，需要根据实际资料或实验数据提出各种假设，然后在这些假设下，利用合适的数学方法列出微分方程。
在实际的微分方程建模过程中，这些方法往往是综合应用的。无论采用哪种方法，通常需要根据实际情况做出一定的假设和简化，并将模型的理论或计算结果与实际情况对比验证，以此来修改模型，使其更准确地描述实际问题，从而达到预测和预报的目的。
习题6.4
1. 题目要求
2.解题过程
解：
(1)
对函数值 s ( t ) s(t) s(t) 的影响因素，根据题目分析，可以得到如下两点：
- 下降速度与 s ( t ) s(t) s(t) 的值的大小成正比，设比例系数为 λ \lambda λ ，所以式子中有一个部分为： − λ s ( t ) -\lambda s(t) −λs(t) 。
- 增长速度与 a ( t ) a(t) a(t) 的值的大小成正比，设比例系数为 μ \mu μ ，所以式子中有一个部分为： + μ a ( t ) +\mu a(t) +μa(t) 。
但是，题目中说到，广告只能影响这种商品在市场上尚未饱和的部分(设饱和量为M)，所以有两个选择：分段函数 或者 渐进函数。
我选择使用渐进函数，渐近线为M，增长部分改为： μ a ( t ) ( 1 − s ( t ) M ) \mu a(t)\left(1-\frac{s(t)}{M}\right) μa(t)(1−Ms(t)) 。
这样，当 s=M 或者 a=0 时，销量都不会上涨，符合题意。
将所有分析合并，建立如下模型：
d
s
d
t
=
μ
a
(
t
)
(
1
−
s
(
t
)
M
)
−
λ
s
(
t
)
\frac{\mathrm{d} s}{\mathrm{d} t}=\mu a(t)\left(1-\frac{s(t)}{M}\right)-\lambda s(t)
dtds=μa(t)(1−Ms(t))−λs(t)
(2)
广告宣传只进行有限时间
τ
\tau
τ ，且广告费为常数a，所以建立下列函数
a
(
t
)
=
{
a
/
c
,
0
<
t
⩽
τ
,
0
,
t
⩾
τ
.
a(t)= \begin{cases}a/c, & 0<t \leqslant \tau, \\ 0, & t \geqslant \tau . \end{cases}
a(t)={a/c,0,0<t⩽τ,t⩾τ.
带入(1)中公式，并移项，有：
d
s
d
t
+
(
λ
+
μ
a
M
τ
)
s
=
μ
a
τ
,
0
<
t
<
τ
,
{\frac{\mathrm{d}s}{\mathrm{d}t}}+\left(\lambda\,+{\frac{\mu a}{M\tau}}\right)s\ ={\frac{\mu a}{\tau}},0\ <t<\tau,
dtds+(λ+Mτμa)s =τμa,0 <t<τ,
(3)
先讨论 0 < t ⩽ τ 0<t \leqslant \tau 0<t⩽τ :
为了运算方便，进行换元：
b
=
λ
+
μ
a
M
τ
,
c
=
μ
a
τ
b = {\lambda\,+\,{\frac{\mu a}{M\tau}} , \,\,\, c={\frac{\mu a}{\tau}}}\\
b=λ+Mτμa,c=τμa
则上述式子化简为：
d
s
d
t
+
b
s
=
c
,
{\frac{\mathrm{d}s}{\mathrm{d}t}\,+\,b s\,=\,c\,,}
dtds+bs=c,
(4)
这是一个一阶非齐次线性微分方程.
假设初始值为
s
0
s_0
s0 ,可以得出结果为:
s
(
t
)
=
c
b
(
1
−
e
−
b
t
)
+
s
0
e
−
b
t
,
0
<
t
⩽
τ
s(t)= \frac{c}{b}\left(1-\mathrm{e}^{-b t}\right)+s_{0} \mathrm{e}^{-b t}, 0<t \leqslant \tau
s(t)=bc(1−e−bt)+s0e−bt,0<t⩽τ
(5)
再讨论 t ⩾ τ t \geqslant \tau t⩾τ :
式子为：
d
s
d
t
=
−
λ
s
,
{\frac{\mathrm{d}s}{\mathrm{d}t}}=-\,\lambda s,
dtds=−λs,
可以得出结果为：
s
(
t
)
=
s
(
τ
)
e
λ
(
τ
−
t
)
,
t
⩾
τ
.
s(t)= s(\tau) \mathrm{e}^{\lambda(\tau-t)}, t \geqslant \tau .
s(t)=s(τ)eλ(τ−t),t⩾τ.
(6)
合并所有结果，得到：
s
(
t
)
=
{
c
b
(
1
−
e
−
b
t
)
+
s
0
e
−
b
t
,
0
<
t
⩽
τ
,
s
(
τ
)
e
λ
(
τ
−
t
)
,
t
⩾
τ
.
s(t)= \left\{\begin{array}{ll} \frac{c}{b}\left(1-\mathrm{e}^{-b t}\right)+s_{0} \mathrm{e}^{-b t}, & 0<t \leqslant \tau, \\ s(\tau) \mathrm{e}^{\lambda(\tau-t)}, & t \geqslant \tau . \end{array}\right.
s(t)={bc(1−e−bt)+s0e−bt,s(τ)eλ(τ−t),0<t⩽τ,t⩾τ.
3.结果
s
(
t
)
s(t)
s(t) 的变化符合下列函数：
s
(
t
)
=
{
c
b
(
1
−
e
−
b
t
)
+
s
0
e
−
b
t
,
0
<
t
⩽
τ
,
s
(
τ
)
e
λ
(
τ
−
t
)
,
t
⩾
τ
.
s(t)= \left\{\begin{array}{ll} \frac{c}{b}\left(1-\mathrm{e}^{-b t}\right)+s_{0} \mathrm{e}^{-b t}, & 0<t \leqslant \tau, \\ s(\tau) \mathrm{e}^{\lambda(\tau-t)}, & t \geqslant \tau . \end{array}\right.
s(t)={bc(1−e−bt)+s0e−bt,s(τ)eλ(τ−t),0<t⩽τ,t⩾τ.
其中，
b
=
λ
+
μ
a
M
τ
,
c
=
μ
a
τ
b = {\lambda\,+\,{\frac{\mu a}{M\tau}} , \,\,\, c={\frac{\mu a}{\tau}}}\\
b=λ+Mτμa,c=τμa
习题6.5(1)
1. 题目要求
2.解题过程
解：
Matlab的工具箱提供了几个解常微分方程的功能函数，如ode45 、ode23 、ode113 。其中，ode45 采用四五阶龙格库塔方法（以下简称RK方法），是解非刚性常微分方程的首选方法； ode23 采用二三阶龙格库塔方法。
高阶常微分方程，必须做变量替换， 化成一阶微分方程组，才能使用Matlab求数值解。
设
y
1
=
y
,
y
2
=
y
′
y_1=y,y_2=y^{\prime}
y1=y,y2=y′ ,则原式化为：
{
y
1
′
=
y
2
,
y
1
(
π
2
)
=
2
,
y
2
′
=
(
n
2
x
2
−
1
)
y
1
−
y
2
x
,
y
2
(
π
2
)
=
−
2
π
.
\left\{\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}\right.
{y1′=y2,y2′=(x2n2−1)y1−xy2,y1(2π)=2,y2(2π)=−π2.
接下来就可以用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
′
y_1=y,y_2=y^{\prime}
y1=y,y2=y′ ,则原式化为：
{
y
1
′
=
y
2
,
y
1
(
0
)
=
1
,
y
2
′
=
−
y
1
c
o
s
x
,
y
2
(
0
)
=
0.
\left\{\begin{array}{ll} y_{1}^{\prime}= \,y_{2}\,,&{{y_{1}(0)~=~1\,,}} \\ y_{2}^{\prime}= -\,y_{1}\mathrm{cos}x\,,&y_{2}(0)~=\,0. \end{array}\right.
{y1′=y2,y2′=−y1cosx,y1(0) = 1,y2(0) =0.
接下来就可以用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方向的速度：
d
x
d
t
=
−
v
2
cos
θ
,
d
y
d
t
=
v
1
−
v
2
sin
θ
.
{\frac{\mathrm{d}x}{\mathrm{d}t}}=-\,v_{2}\cos\theta\,,{\frac{\mathrm{d}y}{\mathrm{d}t}}=v_{1}-\,v_{2}\sin\theta.
dtdx=−v2cosθ,dtdy=v1−v2sinθ.
根据几何关系，我们可以求出：
cos
θ
=
x
x
2
+
y
2
,
sin
θ
=
y
x
2
+
y
2
,
\cos\theta\,=\,{\frac{x}{\sqrt{x^{2}\,+\,y^{2}}}},\sin\theta\,=\,{\frac{y}{\sqrt{x^{2}\,+\,y^{2}}}},
cosθ=x2+y2x,sinθ=x2+y2y,
联立（14）（15）式，得到：
d
x
d
t
=
−
v
2
x
x
2
+
y
2
,
d
y
d
t
=
v
1
−
v
2
y
x
2
+
y
2
,
{\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}}}},
dtdx=−v2x2+y2x,dtdy=v1−v2x2+y2y,
将（16）的两个式子相除，得到：
d
y
d
x
=
−
k
x
2
+
y
2
x
+
y
x
,
0
<
x
<
d
\frac{\mathrm{d}y}{\mathrm{d}x}=-\,k\,{\frac{{\sqrt{x^{2}+y^{2}}}}{x}}+{\frac{y}{x}},0<x<d
dxdy=−kxx2+y2+xy,0<x<d
这是一个一阶齐次微分方程，于是得到初值问题：
{
d
y
d
x
=
−
k
x
2
+
y
2
x
+
y
x
,
0
<
x
<
d
y
(
d
)
=
0
\left\{\begin{array}{ll} {\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}\right.
{dxdy=−kxx2+y2+xy,0<x<dy(d)=0
为了求解这个微分方程，我试图使用Matlab的dsolve函数解决，但是matlab无法解决这个微分方程。
使用Wolfram Mathematica软件的Dsolve函数，成功解出了微分方程的解为：
y
=
d
2
[
(
x
d
)
1
−
k
−
(
x
d
)
1
+
k
]
,
0
<
x
<
d
.
y\,=\,\frac{d}{2}\,\Bigl[\left(\frac{x}{d}\right)^{1-k}\,-\,\left(\frac{x}{d}\right)^{1+k}\,\Bigr],0<x<d.
y=2d[(dx)1−k−(dx)1+k],0<x<d.
3.第(2)问
根据第一问，可知参数方程为：
{
d
x
d
t
=
−
2
x
x
2
+
y
2
,
x
(
0
)
=
d
d
y
d
t
=
1
−
2
y
x
2
+
y
2
,
y
(
0
)
=
0
\left\{\begin{array}{ll} {\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}\right.
⎩
⎨
⎧dtdx=−x2+y22x,x(0)=ddtdy=1−x2+y22y,y(0)=0
求解的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秒，数值解与解析解的对照图（航行曲线）如下图：
文章来源：https://www.toymoban.com/news/detail-841134.html
如果这篇文章对你有帮助，欢迎点赞与收藏~文章来源地址https://www.toymoban.com/news/detail-841134.html
到了这里，关于【数学建模】《实战数学建模：例题与讲解》第五讲-微分方程建模（含Matlab代码）的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
