数学建模算法-day[14]

6.2 传染病预测问题

问题提出
世界上存在很多传染病，如何根据其传播机理预测疾病得传染范围及染病人数等，对传染病的控制意义十分重大。

1.指数传播模型
基本假设

(1) 所研究的区域是一封闭区域，在一个时期内人口总量相对稳定，不考虑人口的迁移（迁入或迁出）。
(2) $t$时刻染病人数$N(t)$是随时间连续变化的、可微的函数。
(3) 每个病人在单位时间内的有效接触（足以使人致病）或传染的人数为$\lambda$（$\lambda >0$为常数）。

模型建立与求解

记$N(t)$为$t$时刻染病人数，则$t+\Delta t$时刻的染病人数为$N(t+\Delta t)$。从$t\to t+\Delta t$时间内，净增加的染病人数为$N(t+\Delta t)-N(t)$，根据基本假设(3),有
$$N(t+\Delta t)-N(t)=\lambda N(t)\Delta t$$
若记$t=0$时刻染病人数为$N_0$，则由基本假设(2)，在上式两端同时除以$\Delta t$,并令$\Delta t\to 0$,得传染病染病人数的微分方程预测模型：
$$\{\begin{array}{rl}& \frac{dN(t)}{dt}=\lambda N(t),t>0,\\ & N(0)=N_0.\end{array}\text{\hspace{1em}(6.13)}$$
利用mathematica就可以求出该模型的解析解为
$$N(t)=N_0{e}^{\lambda t}.$$

DSolve[{n'[t] == \[Lambda]*n[t], n[0] == n0}, n[t], t]

结果分析与评价

模型结果显示传染病的传播是按指数函数。。。

这一节主播学过了，不想过多介绍，有空了再说

6.3 常微分方程的求解

这一节在韩中庚老师的书里是介绍了matlab的，不过本人matlab不是非常适合这种符号计算的，于是本人在这里同时也介绍mathematica的用法，例子还是书上的。

6.3.1 常微分方程的符号解

Matlab符号运算工具箱提供了功能强大的求常微分方程符号解函数dsolve，Methematica提供了功能强大的求常微分方程符号解函数DSolve

例 6.3 求解微分方程$y^{\prime }=-2y+2x^2+2x,y(0)=1$。

解 求得的符号解为$y=e^{-2x}+x^2$。

mathematica代码

DSolve[{y'[x] == -2*y[x] + 2*x^2 + 2*x, y[0] == 1}, y[x], x]

matlab代码

clc;clear
syms y(x)
y=dsolve(diff(y)==-2*y+2*x^2+2*x,y(0)==1)

例 6.4 求解二阶线性微分方程$y^{\prime \prime }-2y^{\prime }+y=e^x,y(0)=1,y^{\prime }(0)=-1$。

解 求得二阶线性微分方程的解为$y=e^x+\frac{x^2{e}^x}{2}-2x{e}^x$。

mathematica代码

DSolve[{y''[x] - 2 y'[x] + y[x] == Exp[x], y[0] == 1, y'[0] == -1}, y[x], x]

matlab代码

clc;clear
syms y(x)
dy=diff(y);
y=dsolve(diff(y,2)-2*diff(y)+y==exp(x),y(0)==1,dy(0)==-1)

例 6.5 已知输入信号为$u(t)=e^{-t}\cos t$，试求下面微分方程的解。
$$y^{(4)}(t)+10y^{(3)}(t)+35y^{\prime \prime }(t)+50y^{\prime }(t)+24y(t)=u^{\prime \prime }(t),y(0)=0,y^{\prime }(0)=-1,y^{\prime \prime }(0)=1,y^{\prime \prime \prime }(0)=1$$

mathematica代码

DSolve[{D[y[t], {t, 4}] + 10 y'''[t] + 35 y''[t] + 50 y'[t] + 24 y[t] == u''[t], y[0] == 0, y'[0] == -1, y''[0] == 1, y'''[0] == 1}, y[t], t]

matlab代码

clc;clear
syms y(t) u(t)
dy=diff(y,1);
dy2=diff(y,2);
dy3=diff(y,3);
u(t)=exp(-t)*cos(t);
y=dsolve(diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y,1)+24*y==diff(u,2),y(0)==0,dy(0)==-1,dy2(0)==1,dy3(0)==1)

下面给出求常微分方程组符号解的例子。

例6.6 试求解下列柯西问题：
$$\{\begin{array}{rl}& \frac{dx}{dt}=Ax,\\ & x(0)=[1,1,1{]}^T\end{array}$$
的解，其中$x(t)=[x_1(t),x_2(t),x_3(t){]}^T,\mathbf{A}=\left[\begin{array}{ccc}3& -1& 1\\ 2& 0& -1\\ 1& -1& 2\end{array}\right]$。

解 求得的解为
$$\mathbf{x}(t)=\left[\begin{array}{c}\frac{4}{3}{\mathrm{e}}^{3t}-\frac{1}{2}{\mathrm{e}}^{2t}+\frac{1}{6}\\ \frac{2}{3}{\mathrm{e}}^{3t}-\frac{1}{2}{\mathrm{e}}^{2t}+\frac{5}{6}\\ \frac{2}{3}{\mathrm{e}}^{3t}+\frac{1}{3}\end{array}\right]$$

本人不会用mathematica求这类问题了，555.没学过

matlab代码

clc;clear
syms x(t) [3,1]  %定义符号向量,x(t)后有空格
A=[3 -1 1;2 0 -1;1 -1 2];
[s1,s2,s3]=dsolve(diff(x)==A*x,x(0)==ones(3,1))

例6.7 试解初值问题
$$\left[\begin{array}{c}{x}_1^{\prime }(t)\\ x_2^{\prime }(t)\\ x_3^{\prime }(t)\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 2& 1& -2\\ 3& 2& 1\end{array}\right]\left[\begin{array}{c}{x}_1(t)\\ x_2(t)\\ x_3(t)\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ {\mathrm{e}}^t\cos 2t\end{array}\right],\left[\begin{array}{c}{x}_1(0)\\ x_2(0)\\ x_3(0)\end{array}\right]=\left[\begin{array}{c}0\\ 1\\ 1\end{array}\right]$$
解 所得的解为

$$\mathbf{x}(t)=\left[\begin{array}{c}{x}_1(t)\\ x_2(t)\\ x_3(t)\end{array}\right]=\left[\begin{array}{c}0\\ {\mathrm{e}}^t\left[\cos (2t)-\sin (2t)-\frac{t\sin (2t)}{2}\right]\\ {\mathrm{e}}^t\left[\cos (2t)+\frac{5\sin (2t)}{4}+\frac{t\cos (2t)}{2}\right]\end{array}\right]$$

matlab代码

clc;clear
syms x(t) [3,1]
A=[1 0 0;2 1 -2;3 2 1];
[s1,s2,s3]=dsolve(diff(x)==A*x+[0,0,exp(t)*cos(2*t)]',x(0)==[0,1,1]')

例6.8 试求常微分方程组
$$\{\begin{array}{rl}& f^{\prime \prime }+g=3\\ & g^{\prime }+f^{\prime }=1\end{array}$$
在初边值条件$f^{\prime }(1)=0,f(0)=0,g(0)=0$时的解。

解 求得的解为
$$\begin{gathered} f(x)=x-\frac{\mathrm{e}-3}{{\mathrm{e}}^2+1}{\mathrm{e}}^x+\frac{\mathrm{e}(3\mathrm{e}+1)}{{\mathrm{e}}^2+1}{\mathrm{e}}^{-x}-3, \\ g(x)=\frac{\mathrm{e}-3}{{\mathrm{e}}^2+1}{\mathrm{e}}^x-\frac{3{\mathrm{e}}^2+\mathrm{e}}{{\mathrm{e}}^2+1}{\mathrm{e}}^{-x}+3. \end{gathered}$$

matlab代码

clc;clear
syms f(x) g(x)
df=diff(f,1);
[f,g]=dsolve(diff(f,2)+g==3,diff(g)+diff(f)==1,df(1)==0,f(0)==0,g(0)==0)

mathematica代码

DSolve[{f''[x] + g[x] == 3, g'[x] + f'[x] == 1, f'[1] == 0, f[0] == 0,g[0] == 0}, {f[x], g[x]}, x]

6.3.2 初值问题的Matlab数值解

matlab的工具箱提供了几个解常微分方程数值解的函数，如ode45,ode23,ode113,其中：ode45采用四五阶龙格-库塔方法(以下简称RK方法),是解非刚性常微分方程的首选方法；ode23采用二三阶RK方法；ode113采用的是多步法，效率一般比ode45高。

在化学工程及自动控制等领域中，所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。具体地说，对一阶线性常微分方程组
$$\frac{dy}{dx}=Ay+\varphi (x),\text{\hspace{1em}(6.26)}$$
式中：$y,\varphi \in R^m;A$为$m$阶方阵。若矩阵$A$的特征值${\lambda }_i(i=1,2,\cdots ,m)$满足关系
$$\{\begin{array}{l}\mathrm{Re}{\lambda }_i<0,i=1,2,\cdots ,m,\\ \underset{1⩽i⩽m}{\max }|\mathrm{Re}{\lambda }_i|\gg \underset{1⩽i⩽m}{\min }|\mathrm{Re}{\lambda }_i|\end{array}$$
则称方程组(6.26)为刚性方程组或stiff方程组，称数
$$s=\underset{1⩽i⩽m}{\max }|\mathrm{Re}{\lambda }_i|/\underset{1⩽i⩽m}{\min }|\mathrm{Re}{\lambda }_i|$$
为刚性比。理论上的分析表明，求解刚性问题所用的数值方法最好是对步长$h$不做任何限制。

Matlab的工具箱提供连几个解刚性常微分方程的功能函数，如ode15s,ode23s,ode23t,ode23tb。

下节将简单介绍ode45求数值解的用法。今天就到这里吧。

参考文献
[1] 司守奎,孙玺青 数学建模算法与应用第3版[M]