两区域分数阶捕食者–猎物种群模型的参数反问题

期刊菜单

两区域分数阶捕食者–猎物种群模型的参数反问题
Parameter Inverse Problem of the Two-Patch Fractional Predator-Prey Population Model

DOI: 10.12677/aam.2024.136248, PDF, HTML, XML, 下载: 25 浏览: 50
作者: 施芳, 张应洪, 刘雪林：贵州师范大学数学科学学院，贵州贵阳
关键词: 分数阶捕食者–猎物种群模型；参数反问题；解的唯一性；有限差分法；同伦正则化算法；Fractional Predator-Prey Population Model； Parameter Inverse Problem； Uniqueness of the Solution； Finite Difference Method； Homotopy Regularization Algorithm

摘要: 本文考虑分数阶捕食者–猎物种群模型以及该模型的参数反问题。介绍正问题解的存在性与唯一性，提出了微分阶数和生长速率的反问题，并对反问题解的存在唯一性给出了证明。利用有限差分法求解正问题获得附加数据，然后根据同伦正则化算法进行参数的数值反演，数值结果表明该算法能够有效求解分数阶捕食者–猎物种群模型参数反问题。

Abstract: This paper considers the fractional predator-prey population model and the parameter inverse problem of the model. The existence and uniqueness of forward problem solutions are introduced, the inverse problems of differential order and growth rate are proposed, and the existence and uniqueness of inverse problem solutions is proved. The finite difference method is used to solve the positive problem to obtain additional data, and then numerically inverse the parameters according to the homotopy regularization algorithm. The numerical results show that the algorithm can effectively solve the parameter inverse problem of the fractional predator-prey population model.

文章引用：施芳, 张应洪, 刘雪林. 两区域分数阶捕食者–猎物种群模型的参数反问题[J]. 应用数学进展, 2024, 13(6): 2595-2604. https://doi.org/10.12677/aam.2024.136248

1. 引言

分数阶微分方程在数学、物理学和生物学中有着重要的研究意义。近年来，分数阶微分方程被广泛的应用到水文学[1] [2]、黏弹性[3]、斑图结构(如蟒蛇，长颈鹿) [4]和混沌动力学等多个领域。分数阶微分系统的理论研究也越来越多，其中包括非线性分数阶微分系统的可解性[5]-[10]。特别地，对分数阶耦合系统的研究有着重大意义，在研究某些复杂问题时会提出该模型，例如，两区域分数阶捕食者–猎物种群模型，该模型是根据捕食者与猎物种群之间的关系建立的，此模型已在文[11]中研究讨论过，证明了系统解的存在性、唯一性和解对初值的依赖。

另一方面，该模型建立之后，对模型中参数的确定是一类很重要的科学问题。目前，反演分数阶微分方程中的参数以及源项问题引起了人们的关注，见Tuan [12]，Zhang [13]，李功胜[14] [15]等。

本文通过差分格式对两区域分数阶捕食者–猎物种群模型正问题进行求解，得到附加数据，进而利用附加数据且根据同伦正则化算法来反演该模型中的微分阶数和猎物再生速率。本文内容安排如下：

第2节给出本篇文章中所要用到的相关理论知识，第3节介绍一般性的分数阶耦合方程组解的存在性、唯一性，第4节基于上述种群模型提出反问题和证明反问题解的唯一性，第5节根据同伦正则化算法反演微分阶数和猎物再生速率，第6节给出结论。

2. 预备知识

定义1 [16] Caputo分数阶导数：

${}_{a}^{C}D_{t}^{α} f (t) = \frac{1}{Γ (n - α)} \int_{a}^{t} \frac{f^{(n)} (τ)}{{(t - τ)}^{α - n + 1}} d τ, n - 1 < α \leq n,$ (2.1)

其中 $α$ 是一个正常数， $n \in N^{+}$ 。

定义2 [16] Caputo分数阶导数的L1插值逼近：对于 $α (0 < α < 1)$ 阶Caputo导数

$_{0}^{C} D_{t}^{α} f (t) = \frac{1}{Γ (1 - α)} \int_{0}^{t} \frac{f^{'} (s)}{{(t - s)}^{α}} d s,$

取正整数N。记 $τ = \frac{T}{N}$ ， $t_{k} = k τ$ ， $0 \leq k \leq N$ 及

$a_{l}^{(α)} = {(l + 1)}^{1 - α} - l^{1 - α}, l \geq 0,$

${_{0}^{C} D_{t}^{α} f (t) |}_{t = t_{n}} = \frac{1}{Γ (1 - α)} \int_{0}^{t_{n}} \frac{f^{'} (t)}{{(t_{n} - t)}^{α}} d t = \frac{1}{Γ (1 - α)} \sum_{k = 1}^{n} \int_{t_{k - 1}}^{t_{k}} \frac{f^{'} (t)}{{(t_{n} - t)}^{α}} d t .$

在区间 $[t_{k - 1}, t_{k}]$ 上对 $f (t)$ 作线性插值，并用其来近似上式的 $f (t)$ ，经过化简可以得到Caputo分数阶导数的L1插值逼近公式：

${}_{0}^{C}D_{t}^{α} f (t_{n}) \approx \frac{τ^{- α}}{Γ (2 - α)} [a_{0}^{(α)} f (t_{n}) - \sum_{k = 1}^{n - 1} (a_{n - k - 1}^{(α)} - a_{n - k}^{(α)}) f (t_{k}) - a_{n - 1}^{(α)} f (t_{0})] .$

为简化记号，记 $D^{α} : = {}_{0}^{C}D_{t}^{α}$ 。

本文考虑两区域分数阶捕食者–猎物种群模型[11]：

$\begin{array}{l} D^{α} U_{1} (t) = r U_{1} (t) (1 - \frac{U_{1}}{k_{2}}) - a U_{1} G_{1} + m_{12} U_{2} - m_{21} U_{1}, \\ D^{α} U_{2} (t) = r U_{2} (t) (1 - \frac{U_{2}}{k_{2}}) - a U_{2} G_{2} - m_{12} U_{2} + m_{21} U_{1}, \\ D^{α} G_{1} (t) = - μ_{1} G_{1} + a e U_{1} G_{1} + n_{12} G_{2} - n_{21} G_{1}, \\ D^{α} G_{2} (t) = - μ_{2} G_{2} + a e U_{2} G_{2} - n_{12} G_{2} + n_{21} G_{1}, \end{array}$ (2.2)

模型中所有的参数均为非负数，其中 $U_{i}$ 表示区域i中的猎物种群数， $G_{i}$ 表示区域i中捕食者种群数， $k_{1}$ ， $k_{2}$ 是猎物分别在区域1和区域2的承载能力，r是猎物的再生长速率，a是捕食者的狩猎率， $m_{i j}$ 是猎物从区域j到区域i的迁移速率， $μ_{1}$ ， $μ_{2}$ 表示区域1与区域2中捕食者的死亡率，e是猎物成为新捕食者的转换率， $n_{i j}$ 是捕食者从区域j到区域i的迁移率。

为了对模型(2.2)分析，给出以下分数阶耦合系统：

${\begin{cases} D^{α} y_{i} (t) = f_{i} (y_{1} (t), y_{2} (t), y_{3} (t), y_{4} (t)), 0 < t \leq T \\ y_{i} (0) = y_{i}^{0}, i = 1, 2, 3, 4. \end{cases}$ (2.3)

定义3 [11]对于模型(2.3)设：

$\overset{\lor}{Y}, \overset{\land}{Y} \in {(C^{1} [0, T])}^{4}$ 且 $\overset{\lor}{Y} = ({\overset{\lor}{y}}_{1}, {\overset{\lor}{y}}_{2}, {\overset{\lor}{y}}_{3}, {\overset{\lor}{y}}_{4})$ 和 $\overset{\land}{Y} = ({\overset{\land}{y}}_{1}, {\overset{\land}{y}}_{2}, {\overset{\land}{y}}_{3}, {\overset{\land}{y}}_{4})$ 满足

$- y_{i}^{0} \geq 0 \geq {\overset{\land}{y}}_{i} - y_{i}^{0}, i = 1, 2, 3, 4.$

并且 $f_{i}$ 是混合拟单调的，设：

$\frac{\partial f_{1}}{\partial y_{2}} \geq 0, \frac{\partial f_{1}}{\partial y_{3}} \leq 0, \frac{\partial f_{1}}{\partial y_{4}} \geq 0, \frac{\partial f_{2}}{\partial y_{1}} \geq 0, \frac{\partial f_{2}}{\partial y_{3}} \geq 0, \frac{\partial f_{2}}{\partial y_{4}} \leq 0,$

$\frac{\partial f_{3}}{\partial y_{1}} \geq 0, \frac{\partial f_{3}}{\partial y_{2}} \geq 0, \frac{\partial f_{3}}{\partial y_{4}} \geq 0, \frac{\partial f_{4}}{\partial y_{1}} \geq 0, \frac{\partial f_{4}}{\partial y_{2}} \geq 0, \frac{\partial f_{4}}{\partial y_{3}} \geq 0,$

且

$D^{α} {\overset{\lor}{y}}_{1} - f_{1} ({\overset{\lor}{y}}_{1}, {\overset{\lor}{y}}_{2}, {\overset{\lor}{y}}_{3}, {\overset{\lor}{y}}_{4}) \geq 0 \geq D^{α} {\overset{\land}{y}}_{1} - f_{1} ({\overset{\land}{y}}_{1}, {\overset{\land}{y}}_{2}, {\overset{\land}{y}}_{3}, {\overset{\land}{y}}_{4}),$

$D^{α} {\overset{\lor}{y}}_{2} - f_{2} ({\overset{\lor}{y}}_{1}, {\overset{\lor}{y}}_{2}, {\overset{\lor}{y}}_{3}, {\overset{\lor}{y}}_{4}) \geq 0 \geq D^{α} {\overset{\land}{y}}_{2} - f_{2} ({\overset{\land}{y}}_{1}, {\overset{\land}{y}}_{2}, {\overset{\land}{y}}_{3}, {\overset{\land}{y}}_{4}),$

$D^{α} {\overset{\lor}{y}}_{3} - f_{3} ({\overset{\lor}{y}}_{1}, {\overset{\lor}{y}}_{2}, {\overset{\lor}{y}}_{3}, {\overset{\lor}{y}}_{4}) \geq 0 \geq D^{α} {\overset{\land}{y}}_{3} - f_{3} ({\overset{\land}{y}}_{1}, {\overset{\land}{y}}_{2}, {\overset{\land}{y}}_{3}, {\overset{\land}{y}}_{4}),$

$D^{α} {\overset{\lor}{y}}_{4} - f_{4} ({\overset{\lor}{y}}_{1}, {\overset{\lor}{y}}_{2}, {\overset{\lor}{y}}_{3}, {\overset{\lor}{y}}_{4}) \geq 0 \geq D^{α} {\overset{\land}{y}}_{4} - f_{4} ({\overset{\land}{y}}_{1}, {\overset{\land}{y}}_{2}, {\overset{\land}{y}}_{3}, {\overset{\land}{y}}_{4}),$

$\overset{\lor}{Y}$ 和 $\overset{\land}{Y}$ 分别是模型(2.3)的上解和下解。

引理1 [11]设 $f_{1}, f_{2}, f_{3}, f_{4}$ 是混合拟单调，满足Lipschitz条件，也满足不等式

$- 4 L > - \frac{Γ (α + 1)}{T^{α}}$

如果上解和下解 $\overset{\lor}{Y}, \overset{\land}{Y} \in {(C^{1} [0, T])}^{4}$ ，且满足 $\overset{\land}{Y} \leq \overset{\lor}{Y}$ ，那么(2.3)式在 $[\overset{\land}{Y}, \overset{\lor}{Y}]$ 上有唯一解。

3. 正问题解的存在唯一性

这一节，将给出正问题模型，并应用上下解法以及相关定理来分析该模型解的存在唯一性。

对于两区域分数阶捕食者–猎物种群模型(2.2)与初值条件

$U_{1} (0) = U_{1}^{0}, U_{2} (0) = U_{2}^{0}, G_{1} (0) = G_{1}^{0}, G_{2} (0) = G_{2}^{0}$

一起构成一个求解 $U_{1} (t)$ ， $U_{2} (t)$ ， $G_{1} (t)$ ， $G_{2} (t)$ 的正问题。为了方便读者更好的理解，我们将简要证明正问题解的存在唯一性。

利用上下解法，显然有

$(\overset{\lor}{U_{1}}, \overset{\lor}{U_{2}}, \overset{\lor}{G_{1}}, \overset{\lor}{G_{2}}) \leq (k_{1}, k_{2}, k_{1}, k_{2})$ , $(\overset{\land}{U_{1}}, \overset{\land}{U_{2}}, \overset{\land}{G_{1}}, \overset{\land}{G_{2}}) \geq (0, 0, 0, 0)$

从而可以得到如下式子：

$\overset{\lor}{U_{i}} - U_{i}^{0} \geq 0 \geq \overset{\land}{U_{i}} - U_{i}^{0}, i = 1, 2,$ (3.1)

$\overset{\lor}{G_{i}} - G_{i}^{0} \geq 0 \geq \overset{\land}{G_{i}} - G_{i}^{0}, i = 1, 2,$ (3.2)

$D^{α} \overset{\lor}{U_{1}} - f_{1} (\overset{\lor}{U_{1}}, \overset{\lor}{U_{2}}, \overset{\land}{G_{1}}, \overset{\lor}{G_{2}}) \geq 0 \geq D^{α} \overset{\land}{U_{1}} - f_{1} (\overset{\land}{U_{1}}, \overset{\land}{U_{2}}, \overset{\lor}{G_{3}}, \overset{\land}{G_{4}}),$ (3.3)

$D^{α} \overset{\lor}{U_{2}} - f_{2} (\overset{\lor}{U_{1}}, \overset{\lor}{U_{2}}, \overset{\lor}{G_{1}}, \overset{\land}{G_{2}}) \geq 0 \geq D^{α} \overset{\land}{U_{2}} - f_{2} (\overset{\land}{U_{1}}, \overset{\land}{U_{2}}, \overset{\land}{G_{3}}, \overset{\lor}{G_{4}}),$ (3.4)

$D^{α} \overset{\lor}{G_{1}} - f_{3} (\overset{\lor}{U_{1}}, \overset{\lor}{U_{2}}, \overset{\lor}{G_{1}}, \overset{\lor}{G_{2}}) \geq 0 \geq D^{α} \overset{\land}{G_{1}} - f_{3} (\overset{\land}{U_{1}}, \overset{\land}{U_{2}}, \overset{\land}{G_{3}}, \overset{\land}{G_{4}}),$ (3.5)

$D^{α} \overset{\lor}{G_{2}} - f_{4} (\overset{\lor}{U_{1}}, \overset{\lor}{U_{2}}, \overset{\lor}{G_{1}}, \overset{\lor}{G_{2}}) \geq 0 \geq D^{α} \overset{\land}{G_{2}} - f_{1} (\overset{\land}{U_{1}}, \overset{\land}{U_{2}}, \overset{\land}{G_{3}}, \overset{\land}{G_{4}}),$ (3.6)

根据定义的上下解法，我们选择合适的值 $\overset{\land}{U_{i}}$ ， $\overset{\land}{G_{i}}$ ， $\overset{\lor}{U_{i}}$ ， $\overset{\lor}{G_{i}}$ $(i = 1, 2)$ 满足不等式(3.1)~(3.6)。

此外通过计算，有

$\frac{\partial f_{1}}{\partial U_{2}} = m_{12} \geq 0, \frac{\partial f_{1}}{\partial G_{1}} = - a U_{1} \leq 0, \frac{\partial f_{1}}{\partial G_{2}} = 0,$

$\frac{\partial f_{1}}{\partial U_{1}} = m_{21} \geq 0, \frac{\partial f_{2}}{\partial G_{1}} = 0, \frac{\partial f_{2}}{\partial G_{2}} = - a G_{2} \leq 0,$

$\frac{\partial f_{3}}{\partial U_{1}} = a e G_{1} \geq 0, \frac{\partial f_{3}}{\partial U_{2}} = 0, \frac{\partial f_{3}}{\partial G_{2}} = n_{12} \geq 0,$

$\frac{\partial f_{4}}{\partial U_{1}} = 0, \frac{\partial f_{4}}{\partial U_{2}} = a e G_{2}, \frac{\partial f_{4}}{\partial G_{1}} = n_{21} \geq 0.$

另外还有：

$\begin{array}{l} \frac{\partial f_{1}}{\partial U_{1}} = r - \frac{2 r}{k_{1}} U_{1} - a P_{1} - m_{21}, \\ \frac{\partial f_{2}}{\partial U_{2}} = r - \frac{2 r}{k_{2}} U_{2} - a P_{2} - m_{12}, \\ \frac{\partial f_{3}}{\partial G_{1}} = - μ_{1} + a e U_{1} - n_{21}, \\ \frac{\partial f_{4}}{\partial G_{2}} = - μ_{2} + a e U_{2} - n_{12} . \end{array}$

其中 $f_{i} (i = 1, 2, 3, 4)$ 是混合拟单调且Lipschitz连续。 $(\overset{\lor}{U_{1}}, \overset{\lor}{U_{2}}, \overset{\lor}{G_{1}}, \overset{\lor}{G_{2}})$ 和 $(\overset{\land}{U_{1}}, \overset{\land}{U_{2}}, \overset{\land}{G_{3}}, \overset{\land}{G_{4}})$ 分别是正问题的上解和下解，又由引理1，可以推出正问题的解具有唯一性。

4. 微分阶数与猎物再生长速率的反问题

4.1. 反问题的形成

分数阶微分方程具有记忆和遗传性质，其中的微分阶数对反常扩散现象具有一定的影响，一般情况下，直接确定微分阶数较为困难，包括正问题中的一些重要系数也是很难直接测量。对此，本文将以确定正问题中微分阶数 $α$ 和猎物再生速率r构造反问题。基于此，设 $t = T$ 时正问题的解为一个观测数据，记为：

$F (T) = (U_{1} (T), U_{2} (T), G_{1} (T), G_{2} (T)),$ (4.1)

从而由正问题及观测数据(4.1)式形成了确定微分阶数和猎物再生速率的参数反问题。

4.2. 反问题解的唯一性

给定一个集合为W，定义 $W = {(α, r) : 0 < α < 1, 0 < r < C}$ ，并且假设未知的参数 $(α, r) \in W$ ，其中C是一个正常数，对于任意给定的 $(α, r) \in W$ ，记 $U_{1}^{α, r}$ ， $U_{2}^{α, r}$ ， $G_{1}^{α, r}$ ， $G_{2}^{α, r}$ 是正问题的解。

定理1. 正问题中， $0 < α < 1$ ， $0 < r < C$ ，这里的C为正常数。 $U_{1}^{α, r}$ ， $U_{2}^{α, r}$ ， $G_{1}^{α, r}$ ， $G_{2}^{α, r}$ 是正问题对应于 $(α, r) \in W$ 的解，那么由观测数据(4.1)式可以唯一确定 $α$ 和r。

证明. 为方便书写，我们考虑一般性的分数阶耦合方程组(2.3)结合附加数据形成的反问题解的唯一性。对于 $(α_{j}, r_{j}) \in W (j = 1, 2)$ ，设 $y_{i}^{α_{1}, r_{1}} (t) = y_{i}^{α_{2}, r_{2}} (t)$ ，下面将证明 $α_{1} = α_{2}$ ， $r_{1} = r_{2}$ 。这里利用反证法先证明 $α_{1} = α_{2}$ ，记 $r : = r_{1} = r_{2}$ ，不妨设 $α_{1} > α_{2}$ 。

对模型(2.3)作Laplace变换[17]可得：

$s^{α} {\tilde{y}}_{i} (s) - s^{α - 1} y_{i}^{0} = {\tilde{f}}_{i} (y_{i} (s)), i = 1, 2, 3, 4;$

将 $(α_{1}, r)$ 与 $(α_{2}, r)$ 对应正问题作Laplace变换：

$s^{α_{1}} {\tilde{y}}_{i}^{α_{1}, r} (s) - s^{α_{1} - 1} y_{i}^{0} = {\tilde{f}}_{i}^{α_{1}, r} (y_{i} (s)),$

$s^{α_{2}} {\tilde{y}}_{i}^{α_{2}, r} (s) - s^{α_{2} - 1} y_{i}^{0} = {\tilde{f}}_{i}^{α_{2}, r} (y_{i} (s)) .$

其中 $y_{i}^{α_{1}, r} (s)$ 与 $y_{i}^{α_{2}, r} (s)$ 为正问题对应于 $(α_{1}, r)$ 与 $(α_{2}, r)$ 的解

继而有：

${\tilde{y}}_{i}^{α_{1}, r} (s) = \frac{s^{α_{1} - 1} y_{i}^{0} + {\tilde{f}}_{i}^{α_{1}, r} (y_{i} (s))}{s^{α_{1}}},$ (4.2)

${\tilde{y}}_{i}^{α_{2}, r} (s) = \frac{s^{α_{2} - 1} y_{i}^{0} + {\tilde{f}}_{i}^{α_{2}, r} (y_{i} (s))}{s^{α_{2}}},$ (4.3)

由假设条件可知，右端项 ${\tilde{f}}_{i}^{α_{1}, r} (y_{i} (s)) = {\tilde{f}}_{i}^{α_{2}, r} (y_{i} (s)), i = 1, 2, 3, 4$ ，用(4.2)~(4.3)式有：

$y_{i}^{α_{1}, r} (s) - y_{i}^{α_{2}, r} (s) = \frac{{\tilde{f}}_{i}^{α_{1}, r} (y_{i} (s))}{s^{α_{1}}} - \frac{{\tilde{f}}_{i}^{α_{2}, r} (y_{i} (s))}{s^{α_{2}}}$

因当 $0 < s < 1$ 时， $s^{α_{1}} < s^{α_{2}}$ ，即是 $y_{i}^{α_{1}, r} (s) < y_{i}^{α_{2}, r} (s)$ ，与假设条件矛盾，因此 $α_{1} > α_{2}$ 不成立，故 $α_{1} = α_{2}$ 。

接下来将证明 $r_{1} = r_{2}$ ，不妨设 $r_{1} > r_{2}$ ，类似于上文的做法，对模型(3.1)作Laplace变换且化简可得：

$r_{1} {\tilde{U}}_{i}^{α, r_{1}} (s) = {\tilde{g}}_{i}^{α, r_{_{1}}} ({\tilde{U}}_{1}^{α, r_{1}} (s), {\tilde{U}}_{2}^{α, r_{1}} (s), {\tilde{G}}_{1}^{α, r_{1}} (s), {\tilde{G}}_{2}^{α, r_{1}} (s)),$ (4.4)

$r_{2} {\tilde{U}}_{i}^{α, r_{2}} (s) = {\tilde{g}}_{i}^{α, r_{_{2}}} ({\tilde{U}}_{1}^{α, r_{2}} (s), {\tilde{U}}_{2}^{α, r_{2}} (s), {\tilde{G}}_{1}^{α, r_{2}} (s), {\tilde{G}}_{2}^{α, r_{2}} (s)),$ (4.5)

此时 $i = 1, 2$ 。由原来的四个方程变为两个方程，同理根据假设知 ${\tilde{g}}_{i}^{α, r_{_{1}}} = {\tilde{g}}_{i}^{α, r_{_{2}}}$ ，再由(4.4)~(4.5)式可得：

${\tilde{U}}_{i}^{α, r_{1}} (s) - {\tilde{U}}_{i}^{α, r_{2}} (s) = \frac{{\tilde{g}}_{i}^{α, r_{_{1}}}}{r_{1}} - \frac{{\tilde{g}}_{i}^{α, r_{_{2}}}}{r_{2}}, (i = 1, 2) .$

当 $r_{1} > r_{2}$ 时，即是 $U_{i}^{α, r_{1}} (s) < U_{i}^{α, r_{2}} (s)$ ，与假设条件矛盾，因此 $r_{1} > r_{2}$ 不成立，故 $r_{1} = r_{2}$ 。综上所述，反问题的解存在唯一性，证毕。

5. 求正问题的解及反演参数

反演参数时，需要附加数据，进而需要获得正问题在 $t = T$ 时刻的解，由于正问题的精确解很难获得，所以使用有限差分法来求其数值解。往下将介绍同伦正则化算法，该算法已在文[18] [19]出现，接着用该算法来反演算例中的参数。

5.1. 求正问题的数值解

在数值实验中，我们对 $0 < t \leq T$ 进行N等分剖分，记 $t_{n} = n τ, (n = 0, 1, 2, \dots, N)$ ，时间步长为 $τ = \frac{T}{N}$ ，

令 $U_{1}^{n} = U_{1} (t_{n})$ ， $U_{2}^{n} = U_{2} (t_{n})$ ， $G_{1}^{n} = G_{1} (t_{n})$ ， $G_{2}^{n} = G_{2} (t_{n})$ 。在正问题中Caputo分数阶导数可以用L1公式来近似代替，右端非线性项用显示Euler法逼近，这样可以得到正问题的差分格式：

$\begin{array}{l} \frac{τ^{- α}}{Γ (2 - α)} [U_{1} (t_{n}) - \sum_{k = 1}^{n - 1} (a_{n - k - 1} - a_{n - k}) U_{1} (t_{k}) - a_{n - 1} U_{1} (t_{0})] \\ = r U_{1}^{n - 1} (1 - \frac{U_{1}^{n - 1}}{k_{1}}) - a U_{1}^{n - 1} G_{1}^{n - 1} + m_{12} U_{2}^{n - 1} - m_{21} U_{1}^{n - 1}, \end{array}$

$\begin{array}{l} \frac{τ^{- α}}{Γ (2 - α)} [U_{2} (t_{n}) - \sum_{k = 1}^{n - 1} (a_{n - k - 1} - a_{n - k}) U_{2} (t_{k}) - a_{n - 1} U_{2} (t_{0})] \\ = r U_{2}^{n - 1} (1 - \frac{U_{2}^{n - 1}}{k_{1}}) - a U_{2}^{n - 1} G_{2}^{n - 1} - m_{12} U_{2}^{n - 1} + m_{21} U_{1}^{n - 1}, \end{array}$

$\begin{array}{l} \frac{τ^{- α}}{Γ (2 - α)} [G_{1} (t_{n}) - \sum_{k = 1}^{n - 1} (a_{n - k - 1} - a_{n - k}) G_{1} (t_{k}) - a_{n - 1} G_{1} (t_{0})] \\ = - μ_{1} G_{1}^{n - 1} + a e U_{1}^{n - 1} G_{1}^{n - 1} + n_{12} G_{2}^{n - 1} - n_{21} G_{1}^{n - 1}, \end{array}$

$\begin{array}{l} \frac{τ^{- α}}{Γ (2 - α)} [G_{2} (t_{n}) - \sum_{k = 1}^{n - 1} (a_{n - k - 1} - a_{n - k}) G_{2} (t_{k}) - a_{n - 1} G_{2} (t_{0})] \\ = - μ_{2} G_{2}^{n - 1} + a e U_{2}^{n - 1} G_{2}^{n - 1} - n_{12} G_{2}^{n - 1} + n_{21} G_{1}^{n - 1}, \end{array}$

令 $θ = τ^{α} Γ (2 - α)$ ，并结合初值条件，可得：

$U_{1} (t_{n}) = θ [r U_{1}^{n - 1} (1 - \frac{U_{1}^{n - 1}}{k_{1}}) - a U_{1}^{n - 1} G_{1}^{n - 1} + m_{12} U_{2}^{n - 1} - m_{21} U_{1}^{n - 1}] + \sum_{k = 1}^{n - 1} (a_{n - k - 1} - a_{n - k}) U_{1} (t_{k}) + a_{n - 1} U_{1}^{0},$

$U_{2} (t_{n}) = θ [r U_{2}^{n - 1} (1 - \frac{U_{2}^{n - 1}}{k_{1}}) - a U_{2}^{n - 1} G_{2}^{n - 1} - m_{12} U_{2}^{n - 1} + m_{21} U_{1}^{n - 1}] + \sum_{k = 1}^{n - 1} (a_{n - k - 1} - a_{n - k}) U_{2} (t_{k}) + a_{n - 1} U_{2}^{0},$

$G_{1} (t_{n}) = θ [- μ_{1} G_{1}^{n - 1} + a e U_{1}^{n - 1} G_{1}^{n - 1} + n_{12} G_{2}^{n - 1} - n_{21} G_{1}^{n - 1}] + \sum_{k = 1}^{n - 1} (a_{n - k - 1} - a_{n - k}) G_{1} (t_{k}) + a_{n - 1} G_{1}^{0},$

$G_{2} (t_{n}) = θ [- μ_{2} G_{2}^{n - 1} + a e U_{2}^{n - 1} G_{2}^{n - 1} - n_{12} G_{2}^{n - 1} + n_{21} G_{1}^{n - 1}] + \sum_{k = 1}^{n - 1} (a_{n - k - 1} - a_{n - k}) G_{2} (t_{k}) + a_{n - 1} G_{2}^{0},$

其中 $a_{k} = {(k + 1)}^{1 - α} + k^{1 - α}, k = 1, 2, \dots$ 。

这样根据初值条件可以求出t在某个时刻的数值解，换言之，可以得出反问题中需要的观测数据。

5.2. 反问题参数反演

这一节，将应用同伦正则化算法[18] [19]来对反问题中微分阶数 $α$ 和猎物再生速率r进行数值反演，首先简单的描述该算法。

给出一些记号 $z = (α, r) \in W$ ，其中 $W = {(α, r) : 0 < α < 1, 0 < r < C}$ ，对所有的 $z_{j} \in W$ ，记 $F (z_{j}) (t)$ 为正问题的解。对反问题的求解实际就是对下面这个算子方程进行求解：

$F (z) (T) = ς (T),$ (5.1)

其中 $F (z) (T)$ 为 $z = (α, r) \in W$ 的计算输出， $ς (T)$ 是由真解对应的正问题解得到的观测数据。对(5.1)定义一个同伦：

$H (z, φ) = (1 - φ) {‖ F (z) (T) - ς (T) ‖}_{*}^{2} + φ {‖ z ‖}^{2},$ (5.2)

这样反问题的求解就转化为求下面极小问题的解：

$\min {(1 - φ) {‖ F (z) (T) - ς (T) ‖}_{L_{2}}^{2} + φ {‖ z ‖}^{2}}$ (5.3)

其中 $φ > 0$ 为同伦参数， $z = (α, r) \in W$ 。根据最佳摄动算法，上述问题又可转变为对应给定的初始值进行反复迭代优化求解，得出最佳摄动量。给定一个 $z_{j}$ ，通过以下方式进行迭代

$z_{j + 1} = z_{j} + δ_{z_{j}}, j = 0, 1, 2, \dots,$ (5.4)

这里的 $δ_{z_{j}}$ 满足下列极小问题：

$P (δ_{z_{j}}) = (1 - φ_{j}) {‖ F (z_{j} + δ_{z_{j}}) (T) - ς (T) ‖}_{L_{2}}^{2} + φ_{j} {‖ δ_{z_{j}} ‖}^{2},$ (5.5)

对 $F (z_{j} + δ_{z_{j}})$ 在 $z_{j}$ 处作Taylor展开，并略去高阶项可得：

$F (z_{j} + δ_{z_{j}}) (T) \approx F (z_{j}) (T) + F^{'} (z_{j}) (T) \cdot δ_{z_{j}} .$

再结合(5.5)式，有：

$P (δ_{z_{j}}) = (1 - φ_{j}) {‖ F^{'} (z_{j}) (T) \cdot δ_{z_{j}} - [ς (T) - F (z_{j} + δ_{z_{j}}) (T)] ‖}_{L_{2}}^{2} + φ_{j} {‖ δ_{z_{j}} ‖}^{2} .$ (5.6)

进一步可将(5.6)写为：

$P (δ_{z_{j}}) = (1 - φ_{j}) {‖ A δ_{z_{j}} - [η - ξ] ‖}^{2} + φ_{j} {‖ δ_{z_{j}} ‖}^{2},$ (5.7)

其中

$A = {(a_{m n})}_{M \times 2}, a_{m n} = \frac{F (z_{j} + ε e_{i}) (t) - F (z_{j}) (t)}{ε}, m = 1, 2, \dots, M, n = 1, 2;$ (5.8)

并且记

$ξ = {(F (z_{j}) (T))}^{T}, η = ς (T) = {(F (z) (T))}^{T}$ (5.9)

这里的 $ξ$ 表示对应于 $z_{j}$ 正问题在 $t = T$ 时的数值解， $η$ 表示对应真解z正问题在 $t = T$ 时的数值解。

在上文中， $e_{i}$ 是二维的单位向量， $ε$ 是数值微分步长，并且对(5.7)式关于 $δ_{z_{j}}$ 进行求导可得：

$((1 - φ_{j}) A^{T} A + φ_{j} I) δ_{z_{j}} = (1 - φ_{j}) A^{T} (η - ξ),$ (5.10)

这里 $φ_{j} \in (0, 1)$ 是依赖于迭代次数的同伦参数，选择同伦参数是运用同伦正则化算法的核心，参考文[18] [19]，定义

$φ_{j} = \frac{1}{1 + \exp (β (j - j_{0}))},$ (5.11)

这里的 $β$ 是调节参数，j为迭代次数， $j_{0}$ 为预估迭代次数。在实际的计算中，观测数据 $η$ 是有随机扰动的，设为 $η^{γ} = η + γ σ$ ，其中， $γ$ 为扰动水平且是正数， $σ$ 是随机向量且每个分量都在 $[- 1, 1]$ 之间。

通过以上算法我们就可以得到 $δ_{z_{j}}$ ，根据迭代式(5.4)可以反复得到最佳摄动量 $δ_{z_{j}}$ ，当 $‖ δ_{z_{j}} ‖ \leq e p s$ 时，停止迭代，这样可以逐渐得到反问题的最优解。整个迭代步骤如下：

第一步：利用真解带入正问题的差分格式进行计算，从而获得观测数据 $η$ ，并且给定控制精度eps，数值微分步长 $ε$ 。

第二步：根据给定的初始迭代 $z_{j} (j = 0, 1, 2, \dots)$ ，按照(5.8)~(5.9)求出矩阵A和 $ξ$ 。

第三步：同伦参数 $φ_{j}$ 依据(5.11)来选取，再利用(5.10)式计算出最佳摄动量 $δ_{z_{j}}$ 。

第四步：判断是否满足 $‖ δ_{z_{j}} ‖ \leq e p s$ ，如果满足，停止迭代，输出反演解 $z_{j + 1}$ 。如果不满足，则转到第二步继续迭代，此时由 $z_{j + 1}$ 代替 $z_{j}$ 。

5.3. 数值算例

算例1. 在数值实验中，取微分阶数 $α = 0.6$ ，猎物再生速率 $r = 2$ ，从而可知真解 $z = (0.6, 2)$ ，继而我们可以利用上文中的差分格式对正问题进行求解来获得附加数据，也即是得到扰动数据 $η^{γ}$ ，最后利用同伦正则化算法进行反演。

正问题中 $N = 1600$ ， $T = 1$ ， $a = 0.5$ ， $m_{12} = 2$ ， $m_{21} = 3$ ， $μ_{1} = 3$ ， $μ_{2} = 4$ ， $e = 0.3$ ， $k_{1} = 100$ ， $k_{2} = 150$ ， $n_{12} = 2.5$ ， $n_{21} = 2.8$ ，初值 $U_{1}^{0} = 20$ ， $U_{2}^{0} = 30$ ， $G_{1}^{0} = 10$ ， $G_{2}^{0} = 8$ ；在反演计算中，预估迭代次数 $j_{0} = 8$ ，调节参数 $β = 0.9$ ，数值微分步长 $ε = 0.01$ ，初始迭代 $z_{0} = (0.3, 0.1)$ ，控制精度 $e p s = 1 e - 6$ 。将根据不同的扰动水平得到反演解，并把反演解十五次的平均值、反演解平均值与真解的相对误差和迭代次数的平均值列在表1中，定义其相对误差为：

$\bar{e r r} = \frac{{‖ \bar{z} - z ‖}_{2}}{{‖ z ‖}_{2}}$

其中， $\bar{z}$ 是连续反演十五次解的的平均值， $\bar{j}$ 代表连续迭代次数的平均值。当 $γ = 0$ 时，表示数据没有扰动。

Table 1. Error results of true and inverse solutions of different γ

表1. 不同γ下，真解与反演解的误差结果

$γ$	$\bar{z}$	$\bar{e r r}$	$\bar{j}$
5%	$(0 .58341,1 .96972)$	2.952e−1	15.4
1%	$(0.60659, 2.02706)$	5.083e−2	15.7
0.1%	$(0.60013, 2.00062)$	3.108e−4	16
0.01%	$(0.60007, 2.00025)$	1.389e−5	16
0	$(0.60000, 2.00000)$	7.394e−11	16

算例2. 在该算例中，微分阶数 $α = 0. 7$ ，猎物再生长速率为 $r = 8$ ，也即是反问题的真解 $z = (0. 7, 8)$ ，初始迭代 $z_{0} = (0.3, 5)$ ，正问题模型中参数 $m_{12} = 0$ ， $n_{12} = 0$ ，控制精度 $e p s = 1 e - 5$ ，其余参数与算例1中一样，反演结果位于表2。

Tablr 2. Error results of true and inverse solutions of different γ

表2. 不同γ下，真解与反演解的误差结果

$γ$	$\bar{z}$	$\bar{e r r}$	$\bar{j}$
5%	$(0.69477, 7.99187)$	2.128e−3	13.45
1%	$(0.69930, 7.99635)$	4.047e−4	13.05
0.1%	$(0.70184, 8.00730)$	3.993e−5	13
0.01%	$(0.70002, 8.00005)$	7.289e−6	13
0	$(0.70000, 8.00000)$	2.809e−8	13

从上面两个算例的结果可以看出，当数据没有扰动的情况下，反演解和真解基本一致，从而反映了该算法的有效性。同样，在数据有扰动的情况下，也可以看出，随着数据的扰动水平下降时，反演解逐渐靠近真解。

6. 总结

本文考虑两区域分数阶捕食者-猎物种群模型，通过有限差分法求解正问题来获得附加数据。进一步，对模型中的微分阶数和猎物再生速率，应用同伦正则算法去反演求解，通过数值算例来验证该算法的有效性。关于本文中反问题数值解的误差分析理论是我们未来要做的工作。

参考文献

[1]	Benson, D.A., Wheatcraft, S.W. and Meerschaert, M.M. (2000) Application of a Fractional Advection-Dispersion Equation. Water Resources Research, 36, 1403-1412. https://doi.org/10.1029/2000wr900031
[2]	Benson, D.A., Wheatcraft, S.W. and Meerschaert, M.M. (2000) The Fractional-Order Governing Equation of Lévy Motion. Water Resources Research, 36, 1413-1423. https://doi.org/10.1029/2000wr900032
[3]	Ahmad, B. and Alsaedi, A. (2010) Existence and Uniqueness of Solutions for Coupled Systems of Higher-Order Nonlinear Fractional Differential Equations. Fixed Point Theory and Applications, 2010, Article No. 364560. https://doi.org/10.1155/2010/364560
[4]	Turing, A.M. (1952) The Chemical Basis of Morphogebesis. Philosophical Transactions of the Royal Society of London B, 237, 37-72.
[5]	Bai, Z. (2010) On Positive Solutions of a Nonlocal Fractional Boundary Value Problem. Nonlinear Analysis: Theory, Methods & Applications, 72, 916-924. https://doi.org/10.1016/j.na.2009.07.033
[6]	Daftardar-Gejji, V. (2005) Positive Solutions of a System of Non-Autonomous Fractional Differential Equations. Journal of Mathematical Analysis and Applications, 302, 56-64. https://doi.org/10.1016/j.jmaa.2004.08.007
[7]	Delbosco, D. and Rodino, L. (1996) Existence and Uniqueness for a Nonlinear Fractional Differential Equation. Journal of Mathematical Analysis and Applications, 204, 609-625. https://doi.org/10.1006/jmaa.1996.0456
[8]	Diethelm, K. (2010) The Analysis of Fractional Differential Equations. Springer-Verlag.
[9]	Li, C. and Deng, W. (2007) Remarks on Fractional Derivatives. Applied Mathematics and Computation, 187, 777-784. https://doi.org/10.1016/j.amc.2006.08.163
[10]	Li, C.F., Luo, X.N. and Zhou, Y. (2010) Existence of Positive Solutions of the Boundary Value Problem for Nonlinear Fractional Differential Equations. Computers & Mathematics with Applications, 59, 1363-1375. https://doi.org/10.1016/j.camwa.2009.06.029
[11]	Zhou, X. and Xu, C. (2017) Analysis of a System of Autonomous Fractional Differential Equations. International Journal of Biomathematics, 10, 1750094. https://doi.org/10.1142/s1793524517500942
[12]	Tuan, V.K. (2011) Inverse Problem for Fractional Diffusion Equation. Fractional Calculus and Applied Analysis, 14, 31-55. https://doi.org/10.2478/s13540-011-0004-x
[13]	Zhang, Y. and Xu, X. (2011) Inverse Source Problem for a Fractional Diffusion Equation. Inverse Problems, 27, 035010. https://doi.org/10.1088/0266-5611/27/3/035010
[14]	Li, G., Gu, W. and Jia, X. (2012) Numerical Inversions for Space-Dependent Diffusion Coefficient in the Time Fractional Diffusion Equation. Journal of Inverse and Ill-Posed Problems, 20, 339-366. https://doi.org/10.1515/jip-2011-0012
[15]	Li, G., Zhang, D., Jia, X. and Yamamoto, M. (2013) Simultaneous Inversion for the Space-Dependent Diffusion Coefficient and the Fractional Order in the Time-Fractional Diffusion Equation. Inverse Problems, 29, 065014. https://doi.org/10.1088/0266-5611/29/6/065014
[16]	孙志忠, 高广花. 分数阶微分方程的差分方法[M]. 第2版. 北京: 科学出版社, 2020.
[17]	杨冰, 李功胜. 一个分数阶生长-抑制系统的参数反问题[J]. 计算数学, 2023, 45(2): 215-229.
[18]	贾现正, 张大利, 李功胜, 池光胜, 李慧玲. 空间-时间分数阶变系数对流扩散方程微分阶数的数值反演[J]. 计算数学, 2014, 36(2): 113-132.
[19]	Zhang, D., Li, G., Jia, X. and Li, H. (2013) Simultaneous Inversion for Space-Dependent Diffusion Coefficient and Source Magnitude in the Time Fractional Diffusion Equation. Journal of Mathematics Research, 5, 65-78. https://doi.org/10.5539/jmr.v5n2p65

为你推荐

友情链接