matlab和Python常微分方程的非标准有限差分法

标准有限差分方案

遇到标准有限差分法的主要问题是数值不稳定性。如果存在有限差分方程的解与微分方程的任何可能解都不对应，则称微分方程的离散模型具有数值不稳定性。 Mickens指出了一些数值不稳定原因。这些原因是：

1. 用不同阶数的离散导数表示微分方程中的导数。 这种情况的一个例子是在表示逻辑模型时

   $$\frac{\mathbf{d}\mathbf{P}}{\mathbf{d}\mathbf{t}}=\mathbf{P}\left(1-\mathbf{P}\right),\mathbf{P}\left(0\right)={\mathbf{P}}_0$$

   使用中心差分方案：

   $$\frac{{\mathbf{P}}^{\mathbf{j}+1}-{\mathbf{P}}^{\mathbf{j}-1}}{2\mathbf{h}}={\mathbf{P}}^{\mathbf{j}}\left(1-{\mathbf{P}}^{\mathbf{j}}\right)\Rightarrow {\mathbf{P}}^{\mathbf{j}+1}={\mathbf{P}}^{\mathbf{j}-1}+2{\mathbf{h}\mathbf{P}}^{\mathbf{j}}\left(1-{\mathbf{P}}^{\mathbf{j}}\right)$$

   下图1显示了使用中心，前向和后向差分方案的逻辑模型$\dot{\mathbf{P}}\left(\mathbf{t}\right)=\mathbf{P}\left(\mathbf{t}\right)\left(1-\mathbf{P}\left(\mathbf{t}\right)\right),\mathbf{P}\left(0\right)=0.5$的解。 该图显示，通过前向和后向差分方案获得的解与原始模型具有相同的行为，而中心有限差分方案与精确解无关。

   图1

2. 使用标准分母函数$\mathbf{h}$。 在点${\mathbf{t}}_{\mathbf{j}}$的一阶和二阶导数近似为：

   $$\begin{gathered} \frac{\mathbf{d}\mathbf{u}}{\mathbf{d}\mathbf{t}}\left({\mathbf{t}}_{\mathbf{i}}\right)\approx \frac{\mathbf{u}\left({\mathbf{t}}_{\mathbf{i}+1}\right)-\mathbf{u}\left({\mathbf{t}}_{\mathbf{i}}\right)}{\mathbf{h}}或\frac{\mathbf{d}\mathbf{u}}{\mathbf{d}\mathbf{t}}\left({\mathbf{t}}_{\mathbf{i}}\right)\approx \frac{\mathbf{u}\left({\mathbf{t}}_{\mathbf{i}}\right)-\mathbf{u}\left({\mathbf{t}}_{\mathbf{i}-1}\right)}{\mathbf{h}} \\ \frac{{\mathbf{d}}^2\mathbf{u}}{{\mathbf{d}\mathbf{t}}^2}\left({\mathbf{t}}_{\mathbf{i}}\right)\approx \frac{\mathbf{u}\left({\mathbf{t}}_{\mathbf{i}+1}\right)-2\mathbf{u}\left({\mathbf{t}}_{\mathbf{i}}\right)+\mathbf{u}\left({\mathbf{t}}_{\mathbf{i}+1}\right)}{{\mathbf{h}}^2} \end{gathered}$$

   在许多情况下，经典分母函数的选择不合适，并且可能导致数值不稳定。 请考虑下面的示例，一阶IVP：

   $$\frac{\mathbf{d}\mathbf{y}\left(\mathbf{t}\right)}{\mathbf{d}\mathbf{t}}=\mathbf{A}\mathbf{y}\left(\mathbf{t}\right)+\mathbf{t},\mathbf{y}\left(0\right)=0$$

   这个IVP确切解是

   $$\mathbf{y}\left(\mathbf{t}\right)=\frac{1}{{\mathbf{A}}^2}\left({\mathbf{e}}^{\mathbf{A}\mathbf{t}}-\left(\mathbf{A}\mathbf{t}+1\right)\right)$$

   间隔$\left[0,1\right]$ 将被划分为1000个子间隔，按照点$0={\mathbf{t}}_0<{\mathbf{t}}_1<\cdots <{\mathbf{t}}_{1000}=1$，每个子间隔的长度$\mathbf{h}=0.001$。 前向欧拉方法，隐式梯形法则和经典四阶龙格库塔法将用于在离散点处针对参数A的不同值解决给定问题。从理论上讲，欧拉方法的误差应为$\mathcal{O}\left(\mathbf{h}\right)=\mathcal{O}\left(10^{-3}\right)$ ，隐式梯形规则的误差应为$\mathcal{O}\left({\mathbf{h}}^2\right)=\mathcal{O}\left(10^{-6}\right)$ ，而四阶龙格库塔法误差应为$\mathcal{O}\left({\mathbf{h}}^4\right)=\mathcal{O}\left(10^{-12}\right)$ 。

   分别计算三种方法的不同A值误差的Python代码是：

   运行代码可获得以下结果：

   runfile(’/media/WindowsData1/PyFiles/
   numinstduetotrivstepfun.py’,
   wdir=’/media/WindowsData1/PyFiles’)
   --------------------------------------------------------
   A Euler Error Imp. Trapz Err RK4 Error
   --------------------------------------------------------
   1 1.35789622e-03 3.56321584e-07 2.30926389e-14
   3 1.00002555e-02 5.19127795e-06 4.50417481e-12
   5 7.35013397e-02 4.52531609e-05 1.53952406e-10
   7 5.39170234e-01 3.52721966e-04 3.11633030e-09
   9 3.94740396e+00 2.65343074e-03 4.88584107e-08
   11 2.88444591e+01 1.97153658e-02 6.58044428e-07
   13 2.10373074e+02 1.45836606e-01 8.01259330e-06
   15 1.53146894e+03 1.07702617e+00 9.07992144e-05
   17 1.11283670e+04 7.94936660e+00 9.75035611e-04
   19 8.07191089e+04 5.86610000e+01 1.00415309e-02
   21 5.84468245e+05 4.32849955e+02 1.00014394e-01
   23 4.22478467e+06 3.19388027e+03 9.69289817e-01
   25 3.04879507e+07 2.35668600e+04 9.18238944e+00
   27 2.19661624e+08 1.73896191e+05 8.53281952e+01
   29 1.58017839e+09 1.28317310e+06 7.79939406e+02
   --------------------------------------------------------
   

   可以看出，尽管步长$\mathbf{h}$保持恒定，但是三种方法的误差都随着$\mathbf{A}$的增加而增加。

3. 使用非线性项的局部逼近：标准有限差分方案使用诸如${\mathbf{x}}^{{\mathbf{j}}^2}或{\mathbf{x}}^{\mathbf{j}+1^2}$的局部逼近来表示诸如${\mathbf{x}}^2\left({\mathbf{t}}_{\mathbf{j}}\right)$ 的非线性项。 Mickens表明，非线性项的这种局部逼近可能导致数值不稳定。 Mickens给出的一个例子是指数衰减模型的解：

   $$\frac{\mathbf{d}\mathbf{P}\left(\mathbf{t}\right)}{\mathbf{d}\mathbf{t}}=-\mathbf{P}\left(\mathbf{t}\right),\mathbf{P}\left(0\right)=0.5$$

   指数衰减模型的前向欧拉离散化由下式给出：

   $${\mathbf{P}}^{\mathbf{n}+1}={\mathbf{P}}^{\mathbf{n}}-{\mathbf{h}\mathbf{P}}^{\mathbf{n}}=\left(1.0-\mathbf{h}\right){\mathbf{P}}^{\mathbf{n}}={\left(1-\mathbf{h}\right)}^{\mathbf{n}}{\mathbf{P}}^0$$

   对于步长$h$的不同值，所得结果如下图2所示。

   图2

   为指数衰减模型获得的解包括具有渐近稳定，周期性和不稳定动力学的解曲线，其中周期性和不稳定解不对应于指数衰减模型的任何真实解。 仅人口图衰减为零的第一个图就与微分方程的真实解一致。

   另一个例子是使用亨恩的方法来求解逻辑模型，

   $$\frac{\mathbf{d}\mathbf{x}\left(\mathbf{t}\right)}{\mathbf{d}\mathbf{t}}=\mathbf{x}\left(\mathbf{t}\right)\left(1-\mathbf{x}\left(\mathbf{t}\right)\right),\mathbf{x}\left(0\right)=0.5$$

   Python代码如下，

   下图3显示了针对亨恩的方法对$\mathbf{h}$的不同值进行逻辑模型的求解。

   图3

   从上图3可以看出，与步长$\mathbf{h}$的某些值相对应，逻辑模型的数值解可能具有周期性和混沌行为，这与模型的任何解都不对应。

非标准有限差分格式的构造规则

与标准有限差分方案一起设计的数值稳定性表明，由于标准模型中的某些固有误差，检索真实信息变得非常昂贵，导致这些标准方案无法准确显示。

考虑到导致标准有限差分方案中数值不稳定的原因，Micken为这种数值方案的设计建立了四个非标准构造规则。 这些规则为：

1. 因为使用与微分方程的阶不同的离散导数会导致数值不稳定，所以离散模型的阶应与微分方程的阶相同。 在此规则下，中心有限差分方案不能用作一阶微分方程离散模型中一阶导数的近似值。 可以使用前向或后向差分方案。

2. 非标准分母函数必须用于连续导数的离散表示。 $\frac{\mathbf{d}\mathbf{y}}{\mathbf{d}\mathbf{t}}$在$\mathbf{t}={\mathbf{t}}_{\mathbf{j}}$的非标准离散表示形式为：

   $$\frac{\mathbf{d}\mathbf{y}}{\mathbf{d}\mathbf{t}}\approx \frac{\mathbf{y}\left({\mathbf{t}}_{\mathbf{j}+1}\right)-\mathbf{\psi }\left(\mathbf{\lambda },\mathbf{h}\right)\mathbf{y}\left({\mathbf{t}}_{\mathbf{j}}\right)}{\mathbf{\varphi }\left(\mathbf{\lambda },\mathbf{h}\right)}$$

   其中$\mathbf{\lambda }$是模型参数的向量，$\mathbf{h}={\mathbf{t}}_{\mathbf{j}+1}-{\mathbf{t}}_{\mathbf{j}}$。 分子和分母函数$\mathbf{\psi }\left(\mathbf{\lambda },\mathbf{h}\right)$ 和$\mathbf{\varphi }\left(\mathbf{\lambda },\mathbf{h}\right)$ 应满足以下特性：

   $$当\mathbf{h}\to 0,\mathbf{\psi }\left(\mathbf{\lambda },\mathbf{h}\right)\to 1与\mathbf{\varphi }\left(\mathbf{\lambda },\mathbf{h}\right)\to \mathbf{h}$$

   举例说明使用非平凡分母如何导致数值不稳定，请考虑指数衰减模型

   $$\frac{\mathbf{d}\mathbf{P}\left(\mathbf{t}\right)}{\mathbf{d}\mathbf{t}}=-\mathbf{P}\left(\mathbf{t}\right),\mathbf{P}\left(0\right)=0.5$$

   在欧拉方法中，假设平凡的分母函数$\mathbf{h}$被分母函数代替

   $$\mathbf{\varphi }\left(\mathbf{h}\right)=\frac{1-{\mathbf{e}}^{-\mathbf{h}}}{\mathbf{h}}$$

   因此，所得离散模型为

   $${\mathbf{P}}^{\mathbf{j}+1}=\left(1-\mathbf{\varphi }\right){\mathbf{P}}^{\mathbf{j}},{\mathbf{P}}^0=0.5$$

   Python代码expdecnsden.py用于求解步长$\mathbf{h}$的不同值的指数衰减模型：

   执行该代码的结果如下图4所示，其中对于不同步长$\mathbf{h}$的值，使用非标准分母函数求解了指数衰减模型。

3. 非线性项应使用非局部近似来近似：像${\mathbf{y}}^2\left(\mathbf{t}\right)$ 的项应由${\mathbf{y}}^{\mathbf{j}}{\mathbf{y}}^{\mathbf{j}+1}$或${\mathbf{y}}^{\mathbf{j}-1}{\mathbf{y}}^{\mathbf{j}+1}$代替${\left({\mathbf{y}}^{\mathbf{j}}\right)}^2$来表示，其中${\left({\mathbf{y}}^{\mathbf{j}}\right)}^2$。 例如，对应于逻辑模型的离散模型

   $$\frac{\mathbf{d}\mathbf{y}}{\mathbf{d}\mathbf{t}}=\mathbf{y}\left(\mathbf{t}\right)\left(1-\mathbf{y}\left(\mathbf{t}\right)\right),\mathbf{y}\left(0\right)=0.5$$

   可能是

   $$\frac{{\mathbf{y}}^{\mathbf{j}+1}-{\mathbf{y}}^{\mathbf{j}}}{\mathbf{h}}={\mathbf{y}}^{\mathbf{j}}\left(1-{\mathbf{y}}^{\mathbf{j}+1}\right)\Rightarrow {\mathbf{y}}^{\mathbf{j}+1}=\frac{{\mathbf{y}}^{\mathbf{j}}\left(1+\mathbf{h}\right)}{1+{\mathbf{h}\mathbf{y}}^{\mathbf{j}}},{\mathbf{y}}^0=0.5$$

   Python代码nonlocapplog.py用于使用非局部逼近法求解逻辑模型。 该代码是

   下图5是通过执行python代码获得的。 它显示了使用不同步长值的逻辑模型的解决方案。 所有这些解都对应于微分方程模型的真实解。 因此，对于步长的所有值，离散模型不会遭受数值不稳定。

   图5

   上图5是通过执行python代码获得的。 它显示了使用不同步长值的逻辑模型的解决方案。 所有这些解都对应于微分方程模型的真实解。 因此，对于步长的所有值，离散模型不会遭受数值上的不稳定性。

4. 离散模型和微分方程的解之间的动态一致性。 这意味着微分方程解所具有的所有性质都必须由离散模型的解所具有。 特殊属性包括正性，单调性，有界性，极限环和其他周期解等。

非标准有限差分方案是基于非标准规则构造的微分方程的离散表示。

精确有限差分法

齐次线性常微分方程的精确有限差分

非线性方程的精确差分

具有线性和幂项的微分方程的精确有限差分

其他非标准有限差分法

详情参阅 - 亚图跨际