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+?(x),\text{\hspace{1em}(6.26)}$$
式中：$y,?\in R^m;A$為$m$階方陣。若矩陣$A$的特征值${\lambda }_i(i=1,2,?,m)$滿足關系
$$\{\begin{array}{l}\mathrm{Re}{\lambda }_i<0,i=1,2,?,m,\\ \underset{1?i?m}{\max }|\mathrm{Re}{\lambda }_i|?\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]