基于复合Heun方法的不确定延迟微分方程的最小二乘估计

Least Squares Estimation of Uncertain Delay Differential Equations Based on the Composite Heun Method

摘要: 文章基于复合Heun方法，研究了不确定延迟微分方程的最小二乘估计问题。利用复合Heun方法，建立了不确定延迟微分方程的差分形式，从而将漂移项中参数的估计问题转化为优化问题，扩散项中参数的估计问题转化为方程的求解。数值算例表明，复合Heun方法的最小二乘估计计算结果优于传统的欧拉方法。最后，利用本文提出的参数估计方法对美洲鹤种群进行了建模分析。

Abstract: Based on the composite Heun method, the article investigates the least squares estimation problem for uncertain delay differential equations. By utilizing the composite Heun method, a difference form is established for uncertain delay differential equations. Then the parameter estimation problem in the drift term is transformed into an optimization problem, and the parameter estimation problem in the diffusion term is converted to solve a system of equations. Numerical examples demonstrate that the least squares estimation using the composite Heun method outperforms the conventional Euler method. Finally, a modeling analysis of the population dynamics of whooping cranes is conducted by applying the parameter estimation method developed in this study.

1. 引言

为了研究广泛存在的不确定现象，基于规范性、对偶性、次可加性和乘积公理，Liu [1]创建了不确定理论。此后，为了描述动态不确定系统，Liu [2]引入了不确定过程的概念，不确定过程本质上是随时间变化的不确定变量。为了处理白噪声，Liu [2]进一步定义了Liu过程，即其是一个平稳的独立增量过程，且增量是正态不确定变量。在此基础上，Liu [3]引入了不确定微分方程。Barbacioru [4]提出了不确定延迟微分方程，并针对一类特殊的不确定延迟微分方程，证明了局部存在唯一的解。Ge和Zhu [5]根据巴拿赫不动点定理，证明了不确定延迟微分方程解的存在唯一性定理。Wang和Ning [6]研究了不确定延迟微分方程的均值稳定性、矩稳定性以及它们之间的相互联系。Jia和Sheng [7]探讨了不确定延迟微分方程的分布稳定性，并建立了分布稳定性的充分条件。

基于观测数据，如何估计不确定微分方程中的未知参数是一个重要问题。Yao和Liu [8]使用矩估计方法来估计不确定微分方程的参数。Liu和Jia [9]基于矩估计方法，分析了不确定延迟微分方程的参数估计问题。Gao、Gao和Yang [10]利用泰勒展开式，得到未知延迟时间的近似值，从而确定参数和延迟时间的估计值。Chen、Yao和Sheng [11]提出了利用最小二乘估计方法，估计不确定微分方程中的未知参数，并根据这一原理，得到了一些特殊类型的不确定微分方程的参数估计值，最后通过解析例子和数值实验说明了方法的有效性。Zhang、Sheng和Wang [12]利用最小二乘估计方法，研究了高阶不确定微分方程的参数估计问题。Liu和Liu [13]提出了极大似然估计方法，并用来估计不确定微分方程的未知参数。

在对不确定微分方程进行参数估计时，需要用差分法对方程进行离散化处理，上述文献的参数估计方法主要基于欧拉差分方法。为了提高精度，Zhou、Tan和Wang [14]提出了复合Heun差分方法，并用于不确定微分方程的参数估计。本文首先利用复合Heun方法，建立了不确定延迟微分方程的差分形式，再基于最小二乘法估计方法估计不确定延迟微分方程中的未知参数，并通过数值例子比较复合Heun方法和欧拉方法的有效性。

2. 复合Heun方法

考虑如下的不确定微分方程

$\text{d}{X}_t=f\left(t,X_t\right)\text{d}t+g\left(t,X_t\right)\text{d}{C}_t$ (1)

其中$f$ 和$g$ 分别是漂移项和扩散项，且是满足利普希茨条件和线性增长条件的连续可微函数，$C_t$ 是Liu过程。对于初始值$X_0$ ，不确定微分方程有唯一解。Zhou、Tan和Wang [14]提出的不确定微分方程(1)的复合Heun差分格式为

$X_{t_{n+1}}=X_{t_n}+{\Delta }_{t_n}\left(\left(1-\theta \right)f\left(t_n,X_{t_n}\right)+\theta f\left(t_{n+1},X_{t_n}+f\left(t_n,X_{t_n}\right){\Delta }_{t_n}\right)\right)+g\left(t_n,X_{t_n}\right)\left(C_{t_{n+1}}-C_{t_n}\right)$ (2)

其中$\theta$ 是松弛系数。显然，复合Heun方法是一种预测校正方法，它的近似度要比欧拉近似更精确。当$\theta =0$ 时，复合Heun方法变成了欧拉方法。

进一步地，本文将复合Heun方法应用到近似不确定延迟微分方程。考虑如下的不确定延迟微分方程

$\text{d}{X}_t=f\left(t,X_t,X_{t-\tau }\right)\text{d}t+g\left(t,X_t,X_{t-\tau }\right)\text{d}{C}_t$ (3)

$f$ 和$g$ 分别是漂移项和扩散项，且是满足利普希茨条件和线性增长条件的连续可微函数，$C_t$ 是Liu过程。类似地，根据复合Heun方法，方程(3)的复合Heun差分格式为

$\begin{array}{l}{X}_{t_{n+1}}=X_{t_n}+{\Delta }_{t_n}\left(\left(1-\theta \right)f\left(t_n,X_{t_n},X_{t_n-\tau }\right)+\theta f\left(t_{n+1},X_{t_n}+hf\left(t_n,X_{t_n},X_{t_n-\tau }\right),X_{t_{n+1}-\tau }\right)\right)\\ +g\left(t_n,X_{t_n},X_{t_n-\tau }\right)\left(C_{t_{n+1}}-C_{t_n}\right)\end{array}$

然后将复合Heun方法同最小二乘估计法相结合，利用最小噪声原理获得不确定延迟微分方程中漂移项和扩散项中的未知参数。

3. 最小二乘估计

考虑如下的不确定延迟微分方程

$\{\begin{array}{l}\text{d}{X}_t=f\left(t,X_t,X_{t-\tau };\alpha \right)\text{d}t+g\left(t,X_t,X_{t-\tau };\beta \right)\text{d}{C}_t,\text{ 0}\le t\le T\\ X_t=\phi \left(t\right),-\tau \le 0\le t\end{array}$ (4)

其中$f$ 和$g$ 分别是漂移项和扩散项，且满足利普希茨条件和线性增长条件的连续可微函数，$\alpha$ 和$\beta$ 均是包含未知参数的向量，$\tau$ 是给定的延迟时间，$C_t$ 是Liu过程，$\phi \left(t\right)$ 是给定的已知函数。下面利用$X_t$ 的观测数据对未知向量$\alpha$ 和$\beta$ 进行估计。将区间$\left[0,T\right]$ 进行$N$ 等分，记步长$h=T/N$ ，节点$t_i=ih\left(i=0,1,\cdots ,N\right)$ ，即$0=t_0<t_1<\cdots <t_N\le T$ 。假设存在正整数$m$ ，使得$\tau =mh$ 成立。根据Zhou、Tan和Wang [14]提出的不确定微分方程的复合Heun方法，方程(4)的复合Heun差分格式是

$\{\begin{array}{l}{X}_{t_{n+1}}=X_{t_n}+h\left(\left(1-\theta \right)f\left(t_n,X_{t_n},X_{t_{n-m}};\alpha \right)+\theta f\left(t_{n+1},X_{t_n}+hf\left(t_n,X_{t_n},X_{t_{n-m}};\alpha \right),X_{t_{n-m+1}};\alpha \right)\right)\\ +g\left(t_n,X_{t_n},X_{t_{n-m}};\beta \right)\left(C_{t_{n+1}}-C_{t_n}\right),\text{ 0}\le n\le N-1\\ X_{t_{n-m}}=\phi \left(t_{n-m}\right),n-m\le 0\end{array}$ (5)

其中$\theta$ 是松弛系数。

定理2.1 假设$X_t$ 有$N+1$ 个观测数据$x_{t_0},x_{t_1},\cdots ,x_{t_N}$ ，方程(5)中未知向量$\alpha$ 的最小二乘估计值$\widehat{\alpha }$ 为以下优化问题的解

$\underset{\alpha }{\min }{{\displaystyle \sum _{n=0}^{N-1}\left(x_{t_{n+1}}-x_{t_n}-h\left(\left(1-\theta \right)f\left(t_n,x_{t_n},x_{t_{n-m}};\alpha \right)+\theta f\left(t_{n+1},x_{t_n}+hf\left(t_n,x_{t_n},x_{t_{n-m}};\alpha \right),x_{t_{n-m+1}};\alpha \right)\right)\right)}}^2$ (6)

未知向量$\beta$ 的最小二乘估计值$\widehat{\beta }$ 为以下方程的解

$\begin{array}{l}{\displaystyle \sum _{n=0}^{N-1}g{\left(t_n,x_{t_n},x_{t_{n-m}};\beta \right)}^2{h}^2=}{\displaystyle \sum _{n=0}^{N-1}(x_{t_{n+1}}-x_{t_n}-h(\left(1-\theta \right)f\left(t_n,x_{t_n},x_{t_{n-m}};\widehat{\alpha }\right)}\\ {+\theta f\left(t_{n+1},x_{t_n}+hf\left(t_n,x_{t_n},x_{t_{n-m}};\widehat{\alpha }\right),x_{t_{n-m+1}};\widehat{\alpha }\right)))}^2\end{array}$ (7)

证明：由复合Heun差分格式(5)中的第一个方程，可以得到

$\begin{array}{l}g\left(t_n,X_{t_n},X_{t_{n-m}};\beta \right)\left(C_{t_{n+1}}-C_{t_n}\right)=X_{t_{n+1}}-X_{t_n}-h(\left(1-\theta \right)f\left(t_n,X_{t_n},X_{t_{n-m}};\alpha \right)\\ +\theta f\left(t_{n+1},X_{t_n}+hf\left(t_n,X_{t_n},X_{t_{n-m}};\alpha \right),X_{t_{n-m+1}};\alpha \right))\end{array}$ (8)

上式等号左边项可以视为噪声项，应该尽可能的小，因此可以通过最小化等号的右边项来找到$\alpha$ 的最小二乘估计。对于$X_t$ 的观测值$x_{t_0},x_{t_1},\cdots ,x_{t_N}$ ，$\alpha$ 的最小二乘估计值就是下面优化问题的解

$\underset{\alpha }{\min }{{\displaystyle \sum _{n=0}^{N-1}\left(x_{t_{n+1}}-x_{t_n}-h\left(\left(1-\theta \right)f\left(t_n,x_{t_n},x_{t_{n-m}};\alpha \right)+\theta f\left(t_{n+1},x_{t_n}+hf\left(t_n,x_{t_n},x_{t_{n-m}};\alpha \right),x_{t_{n-m+1}};\alpha \right)\right)\right)}}^2$ (9)

令$\widehat{\alpha }$ 代表$\alpha$ 的估计值，那么根据式(8)，未知向量$\beta$ 的最小二乘估计值$\widehat{\beta }$ 满足下面的方程

$\begin{array}{l}E\left[{\displaystyle \sum _{n=0}^{N-1}g{\left(t_n,x_{t_n},x_{t_{n-m}};\beta \right)}^2{\left(C_{t_{n+1}}-C_{t_n}\right)}^2}\right]={\displaystyle \sum _{n=0}^{N-1}(x_{t_{n+1}}-x_{t_n}-h(\left(1-\theta \right)f\left(t_n,x_{t_n},x_{t_{n-m}};\widehat{\alpha }\right)}\\ {+\theta f\left(t_{n+1},x_{t_n}+hf\left(t_n,x_{t_n},x_{t_{n-m}};\widehat{\alpha }\right),x_{t_{n-m+1}};\widehat{\alpha }\right)))}^2\end{array}$

由于$C_t$ 是一个独立增量过程，所以增量$C_t-C_s$ 是正态不确定变量，其期望值为0，方差为${\left(t-s\right)}^2$ 。因此，得到

$\begin{array}{l}E\left[{\displaystyle \sum _{n=0}^{N-1}g{\left(t_n,x_{t_n},x_{t_{n-m}};\beta \right)}^2{\left(C_{t_{n+1}}-C_{t_n}\right)}^2}\right]\\ \text{=}{\displaystyle \sum _{n=0}^{N-1}E\left[g{\left(t_n,x_{t_n},x_{t_{n-m}};\beta \right)}^2{\left(C_{t_{n+1}}-C_{t_n}\right)}^2\right]}={\displaystyle \sum _{n=0}^{N-1}g{\left(t_n,x_{t_n},x_{t_{n-m}};\beta \right)}^{2}E{\left[\left(C_{t_{n+1}}-C_{t_n}\right)\right]}^2}\\ ={\displaystyle \sum _{n=0}^{N-1}g{\left(t_n,x_{t_n},x_{t_{n-m}};\beta \right)}^2{\left(t_{n+1}-t_n\right)}^2}={\displaystyle \sum _{n=0}^{N-1}g{\left(t_n,x_{t_n},x_{t_{n-m}};\beta \right)}^2{h}^2}\end{array}$

即最小二乘估计值$\widehat{\beta }$ 满足下面的方程

$\begin{array}{l}{\displaystyle \sum _{n=0}^{N-1}g{\left(t_n,x_{t_n},x_{t_{n-m}};\beta \right)}^2{h}^2}={\displaystyle \sum _{n=0}^{N-1}(x_{t_{n+1}}-x_{t_n}-h(\left(1-\theta \right)f\left(t_n,x_{t_n},x_{t_{n-m}};\widehat{\alpha }\right)}\\ {+\theta f\left(t_{n+1},x_{t_n}+hf\left(t_n,x_{t_n},x_{t_{n-m}};\widehat{\alpha }\right),x_{t_{n-m+1}};\widehat{\alpha }\right)))}^2\end{array}$ (10)

例1不确定延迟微分方程

$\{\begin{array}{l}\text{d}{X}_t=\alpha X_{t-1}\text{d}t+\beta \text{d}{C}_t,0\le t\le T\\ X_t=1,-1\le t\le 0\end{array}$ (11)

其中$\alpha$ 和$\beta$ 是未知参数，并且$\alpha ,\beta$ 均大于0。利用$X_t$ 在$t_0,t_1,\cdots ,t_N$ 时刻观测数据，确定$\alpha$ 和$\beta$ 的最小二乘估计值。

解：基于复合Heun方法，得到微分方程的差分格式

$X_{t_{n+1}}=X_{t_n}+h\left(\left(1-\theta \right)\alpha X_{t_{n-1}}+\theta \alpha X_{t_n}\right)+\beta \left(C_{t_{n+1}}-C_{t_n}\right)$

即

$\beta \left(C_{t_{n+1}}-C_{t_n}\right)=X_{t_{n+1}}-X_{t_n}-h\left(\left(1-\theta \right)\alpha X_{t_{n-1}}+\theta \alpha X_{t_n}\right)$

通过求解下面的优化问题

$\underset{\alpha }{\min }{\displaystyle \sum _{n=0}^{N-1}{\left(x_{t_{n+1}}-x_{t_n}-\alpha h\left(\left(1-\theta \right)x_{t_{n-1}}+\theta x_{t_n}\right)\right)}^2}$

得到$\alpha$ 的最小二乘估计值为

$\widehat{\alpha }={\displaystyle \sum _{n=0}^{N-1}\left(x_{t_{n+1}}-x_{t_n}\right)\left(\left(1-\theta \right)x_{t_{n-1}}+\theta x_{t_n}\right)}/{\displaystyle \sum _{n=0}^{N-1}{\left(\left(1-\theta \right)x_{t_{n-1}}+\theta x_{t_n}\right)}^{2}h}$

未知参数$\beta$ 的最小二乘估计值$\widehat{\beta }$ 满足

${\displaystyle \sum _{n=0}^{N-1}{\beta }^2{h}^2}={\displaystyle \sum _{n=0}^{N-1}{\left(x_{t_{n+1}}-x_{t_n}-\widehat{\alpha }h\left(\left(1-\theta \right)x_{t_{n-1}}+\theta x_{t_n}\right)\right)}^2}$

求解得到

$\widehat{\beta }=\sqrt{{\displaystyle \sum _{n=0}^{N-1}{\left(x_{t_{n+1}}-x_{t_n}-\widehat{\alpha }h\left(\left(1-\theta \right)x_{t_{n-1}}+\theta x_{t_n}\right)\right)}^2/N{h}^2}}$

例2不确定延迟微分方程

$\{\begin{array}{l}\text{d}{X}_t=\alpha X_t\text{d}t+\beta X_{t-1}\text{d}{C}_t,0\le t\le T\\ X_t=1,-1\le t\le 0\end{array}$ (12)

其中$\alpha$ 和$\beta$ 是未知参数，并且$\alpha ,\beta$ 均大于0。利用$X_t$ 在$t_0,t_1,\cdots ,t_N$ 时刻观测数据，确定$\alpha$ 和$\beta$ 的最小二乘估计值。

解：基于复合Heun方法，得到微分方程的差分格式

$X_{t_{n+1}}=X_{t_n}+h\left(\left(1-\theta \right)\alpha X_{t_n}+\theta \alpha \left(X_{t_n}+h\alpha X_{t_n}\right)\right)+\beta X_{t_{n-1}}\left(C_{t_{n+1}}-C_{t_n}\right)$

即

$\beta X_{t_{n-1}}\left(C_{t_{n+1}}-C_{t_n}\right)=X_{t_{n+1}}-X_{t_n}-h\left(\left(1-\theta \right)\alpha X_{t_n}+\theta \alpha \left(X_{t_n}+h\alpha X_{t_n}\right)\right)$

通过求解下面的优化问题

$\underset{\alpha }{\min }{\displaystyle \sum _{n=0}^{N-1}{\left(x_{t_{n+1}}-x_{t_n}-\alpha h\left(1+\theta h\alpha \right)x_{t_n}\right)}^2}$

得到$\alpha$ 的最小二乘估计值$\widehat{\alpha }$ 。未知参数$\beta$ 的最小二乘估计值$\widehat{\beta }$ 满足

${\displaystyle \sum _{n=0}^{N-1}{\beta }^2{h}^2{x}_{t_{n-1}}^2}={\displaystyle \sum _{n=0}^{N-1}{\left(x_{t_{n+1}}-x_{t_n}-\widehat{\alpha }h\left(1+\theta h\alpha \right)x_{t_n}\right)}^2}$

求解得到

$\widehat{\beta }=\sqrt{{\displaystyle \sum _{n=0}^{N-1}{\left(x_{t_{n+1}}-x_{t_n}-\widehat{\alpha }h\left(1+\theta h\alpha \right)x_{t_n}\right)}^2}/{\displaystyle \sum _{n=0}^{N-1}{h}^2}{x}_{t_{n-1}}^2}$

4. 数值实验

本节中，通过两个数值算例，验证本文提出算法的有效性。给出如下算法步骤：

步骤1 对于不确定变量$\eta ~ℒ\left(0,1\right)$ ，即$\eta$ 的不确定分布为

$\Psi \left(x\right)=\{\begin{array}{l}0,�x\le 0�\\ x,�0<x<1�\\ 1,�x\ge 1�\end{array}$

得到其样本${\eta }_n\left(n=1,2,\cdots ,N\right)$ 。

步骤2 基于不确定标准正态分布的逆不确定分布，计算得到

${\xi }_n=\frac{\sqrt{3}}{\pi }\ln \frac{{\eta }_n}{1-{\eta }_n},n=1,2,\cdots ,N$

那么${\xi }_1,{\xi }_2,\cdots ,{\xi }_N$ 可视为标准正态不确定变量$\xi$ 的样本。特别地，${\xi }_1,{\xi }_2,\cdots ,{\xi }_N$ 可视为$\frac{C_{t_{n+1}}-C_{t_n}}{h}$ 的样本。

步骤3 对于给定的初值$x_{t_0}=\phi \left(t_0\right),x_{t_{0-m}}=\phi \left(t_0-mh\right)$ ，利用方程(1)的欧拉方法，计算得到$X_t$ 的观测值$x_{t_0},x_{t_1},\cdots ,x_{t_N}$ 。

步骤4根据步骤3得到$X_t$ 的观测数据，利用基于复合Heun方法的最小二乘估计算法，计算未知参数的估计值。

为了评估估计的准确性，定义偏差函数如下

$\text{Bias}\left(\mu ,\widehat{\mu }\right)=\left|\mu -\widehat{\mu }\right|$

其中$\widehat{\mu }$ 是$\mu$ 的估计值。显然，偏差越小，估计的效果越好。

例3不确定过程$X_t$ 满足如下的不确定延迟微分方程

$\{\begin{array}{l}\text{d}{X}_t=\alpha X_{t-1}\text{d}t+\beta \text{d}{C}_t,\text{ 0}\le t\le 2\\ X_t=1,-1\le t\le 0\end{array}$ (13)

其中$\alpha$ 和$\beta$ 是未知参数，且$\alpha ,\beta$ 均大于0。为比较复合Heun方法和欧拉方法，设定参数$\alpha =0.5,\beta =0.2$ ，选取步长$h=0.1$ 。利用上述算法，计算得到$X_t$ 的观测数据，如表1所示。取$\theta =0.9$ ，利用复合Heun方法计算得到$\alpha$ 和$\beta$ 的最小二乘估计值分别为

$\widehat{\alpha }=0.5258,\widehat{\beta }=0.2687$

因此，得到估计方程

$\{\begin{array}{l}\text{d}{X}_t=0.5258X_{t-1}\text{d}t+0.2687\text{d}{C}_t,\text{ 0}\le t\le 2\\ X_t=1,-1\le t\le 0\end{array}$ (14)

Table 1. Sample data in Example 3

表1. 例3的样本数据

$n$	0	1	2	3	4	5	6
$t_n$	0	0.1	0.2	0.3	0.4	0.5	0.6
$x_{t_n}$	1	1.1000	1.2068	1.2814	1.3775	1.4498	1.5143
$n$	7	8	9	10	11	12	13
$t_n$	0.7	0.8	0.9	1.0	1.1	1.2	1.3
$x_{t_n}$	1.6162	1.6910	1.7586	1.8385	1.9361	2.0306	2.1399
$n$	14	15	16	17	18	19	20
$t_n$	1.4	1.5	1.6	1.7	1.8	1.9	2.0
$x_{t_n}$	2.3009	2.4437	2.5603	2.7187	2.8243	2.9281	3.0602

而利用欧拉方法计算得到$\alpha$ 和$\beta$ 的估计值为

$\widehat{\alpha }=0.5503,\widehat{\beta }=0.2718$

相应的估计方程为

$\{\begin{array}{l}\text{d}{X}_t=0.5503X_{t-1}\text{d}t+0.2718\text{d}{C}_t,0\le t\le 2\\ X_t=1,-1\le t\le 0\end{array}$ (15)

Figure 1. Sample data and α-paths of $X_t$ in Example 3

图1. 例3的样本数据和$X_t$ 的α轨道

由图1，可以发现，$X_t$ 的观测数据均落在估计方程(14)和(15)的0.01-轨道与0.99-轨道之间的区域，这说明两种方法均是可行的。

比较两种方法得到的估计值，从表2可以发现，基于复合Heun方法的最小二乘估计方法优于基于欧拉方法的最小二乘估计方法。

Table 2. The biases of the least squares estimation based on the composite Heun method and Euler method

表2. 基于复合Heun方法和欧拉方法的最小二乘估计的偏差

	Bias $\left(\alpha ,\widehat{\alpha }\right)$	Bias $\left(\beta ,\widehat{\beta }\right)$
复合Heun方法	0.0258	0.0687
欧拉方法	0.0503	0.0718

例4设不确定过程$X_t$ 满足如下的不确定延迟微分方程

$\{\begin{array}{l}\text{d}{X}_t=\alpha X_{t-1}\text{d}t+\beta X_{t-1}\text{d}{C}_t,0\le t\le 4\\ X_t=1,-1\le t\le 0\end{array}$ (16)

其中$\alpha$ 和$\beta$ 是未知参数，并且$\alpha ,\beta$ 均大于0。设定参数$\alpha =-0.2,\beta =0.3$ ，选取步长$h=0.2$ 。利用初始条件，计算得到$X_t$ 的20组观测值，如表3所示。

Table 3. Sample data in Example 4

表3. 例4的样本数据

$n$	0	1	2	3	4	5	6
$t_n$	0	0.2	0.4	0.6	0.8	1.0	1.2
$x_{t_n}$	1	0.9000	0.8491	0.7737	0.7732	0.7132	0.7347
$n$	7	8	9	10	11	12	13
$t_n$	1.4	1.6	1.8	2.0	2.2	2.4	2.6
$x_{t_n}$	0.6859	0.6214	0.6159	0.6547	0.6751	0.6674	0.6732
$n$	14	15	16	17	18	19	20
$t_n$	2.8	3.0	3.2	3.4	3.6	3.8	4.0
$x_{t_n}$	0.6878	0.7071	0.6521	0.6359	0.5854	0.5552	0.5566

基于复合Heun方法，$\alpha$ 的最小二乘估计值$\widehat{\alpha }$ 是如下优化问题的解

$\underset{\alpha }{\min }{\displaystyle \sum _{n=0}^{N-1}{\left(x_{t_{n+1}}-x_{t_n}-\alpha h\left(\left(1-\theta \right)x_{t_{n-1}}+\theta x_{t_n}\right)\right)}^2}$

取$\theta =0.8$ ，得到

$\widehat{\alpha }=-0.1767$

参数$\beta$ 的最小二乘估计值$\widehat{\beta }$ 满足方程

${\displaystyle \sum _{n=0}^{N-1}{\beta }^2{h}^2{x}_{t_{n-1}}^2}={\displaystyle \sum _{n=0}^{N-1}{\left(x_{t_{n+1}}-x_{t_n}-\widehat{\alpha }h\left(\left(1-\theta \right)x_{t_{n-1}}+\theta x_{t_n}\right)\right)}^2}$

求得

$\widehat{\beta }=0.1075$

因此，估计方程为

$\{\begin{array}{l}\text{d}{X}_t=-0.1767X_{t-1}\text{d}t+0.1075\text{d}{C}_t,0\le t\le 4\\ X_t=1,-1\le t\le 0\end{array}$ (17)

类似地，利用欧拉方法，计算得到$\alpha$ 和$\beta$ 的估计值为

$\widehat{\alpha }=-0.1724,\widehat{\beta }=0.1074$

相应的估计方程为

$\{\begin{array}{l}\text{d}{X}_t=-0.1724X_{t-1}\text{d}t+0.1074\text{d}{C}_t,0\le t\le 4\\ X_t=1,-1\le t\le 0\end{array}$ (18)

Table 4. The biases of the least squares estimation based on the composite Heun method and Euler method

表4. 基于复合Heun方法和欧拉方法的最小二乘估计的偏差

	Bias $\left(\alpha ,\widehat{\alpha }\right)$	Bias $\left(\beta ,\widehat{\beta }\right)$
复合Heun方法	0.0233	0.1925
欧拉方法	0.0276	0.1926

由图2，可以发现，$X_t$ 的观测数据均落在估计方程(17)和(18)的0.01-轨道与0.90-轨道之间的区域，这说明两种方法均是可行的。

Figure 2. Sample data and α-paths of $X_t$ in Example 4

图2. 例4的样本数据和$X_t$ 的α轨道

比较两种方法得到的估计值，从表4可以发现，基于复合Heun方法的最小二乘估计方法优于基于欧拉方法的最小二乘估计方法。

5. 实证分析

在本节中，利用1939年至2010年间美国德克萨斯州阿兰萨斯国家野生动物保护区周围越冬地的美洲鹤种群数量[15]，建立不确定延迟微分方程模型。假设美洲鹤种群数量满足Gao、Gao和Yang [10]提出的不确定延迟微分方程

$\{\begin{array}{l}\text{d}{X}_t=\alpha X_t\left(1-\frac{X_{t-1}}{K}\right)\text{d}t+\beta X_t\left(1-\frac{X_{t-1}}{K}\right)\text{d}{C}_t,t\in \left[0,71\right]\\ X_t=18,t\in \left[-1,0\right]\end{array}$ (19)

其中$\alpha ,\beta ,K$ 是未知参数，且$t_{n+1}-t_n=1$ 。美洲鹤种群数量$X_t$ 的观察值为$x_{t_0},x_{t_1},\cdots ,x_{t_{71}}$ ，如表5所示。

Table 5. Whooping crane population from 1939 to 2010

表5. 1939年至2010年美洲鹤种群数量

22	26	16	19	21	18	22	25	31	30	34	31
25	21	24	21	28	24	26	32	33	36	39	32
33	42	44	43	48	50	56	57	59	51	49	49
57	69	72	75	76	78	73	73	75	86	97	110
134	138	146	146	132	136	143	133	158	160	182	183
188	180	176	185	194	217	220	237	266	270	264	283

基于复合Heun方法，利用最小二乘估计法计算未知参数。估计值$\widehat{\alpha }$ 和$\widehat{K}$ 是以下优化问题的解

$\underset{\alpha ,K}{\min }{\displaystyle \sum _{n=0}^{N-1}{\left(x_{t_{n+1}}-x_{t_n}-\left(1-\theta \right)\alpha x_{t_n}\left(1-\frac{x_{t_{n-1}}}{K}\right)-\theta \alpha \left(x_{t_n}+\alpha x_{t_n}\left(1-\frac{x_{t_{n-1}}}{K}\right)\right)\left(1-\frac{x_{t_n}}{K}\right)\right)}^2}$

取$\theta =0.7$ ，得到

$\widehat{\alpha }=0.0590,K=499.9994$

估计值$\widehat{\beta }$ 满足

$\begin{array}{l}{{\displaystyle \sum _{n=0}^{N-1}{\beta }^2{x}_{t_n}^2\left(1-\frac{x_{t_{n-1}}}{K}\right)}}^2={\displaystyle \sum _{n=0}^{N-1}(x_{t_{n+1}}-x_{t_n}-\left(1-\theta \right)\widehat{\alpha }{x}_{t_n}\left(1-\frac{x_{t_{n-1}}}{K}\right)}\\ -{\theta \alpha \left(x_{t_n}+\alpha x_{t_n}\left(1-\frac{x_{t_{n-1}}}{K}\right)\right)\left(1-\frac{x_{t_n}}{K}\right))}^2\end{array}$

求解可得

$\widehat{\beta }=0.1048$

因此估计方程为

$\{\begin{array}{l}\text{d}{X}_t=0.0590X_t\left(1-\frac{X_{t-1}}{499.9994}\right)\text{d}t+0.1048X_t\left(1-\frac{X_{t-1}}{499.9994}\right)\text{d}{C}_t,t\in \left[0,71\right]\\ X_t=18,t\in \left[-1,0\right]\end{array}$

从图3，可以发现$X_t$ 的观测数据均落在估计方程0.01-轨道与0.99-轨道之间的区域，这说明估计方程是合理的。

Figure 3. Whooping crane aerial survey data from 1939 to 2010 and α-path of $X_t$

图3. 1939年至2010年美洲鹤的航测数据及$X_t$ 的α轨道

6. 结论

不确定微分方程可以用于描述不确定性动态系统，在实际应用中，参数估计是其重要的研究方向。以往的研究大多是基于欧拉方法对不确定微分方程进行离散，得到方程的差分形式，从而对未知参数进行估计。为了提高差分方程的精度，本文利用复合Heun方法，将其应用于不确定延迟微分方程的最小二乘估计。数值计算结果表明，与欧拉方法相比，复合Heun方法具有更高的精度。

NOTES

^*通讯作者。