0 引言

Logistic方程是一类重要的生物数学模型, 其描述了资源总量约束条件下的种群增长行为。单位时间种群密度的增长变化满足^[1-2]

$ \frac{\partial u}{\partial t}=D \Delta\left(\frac{u(x, t)}{P(x)}\right)+r u(x, t)\left(1-\frac{u(x, t)}{K(x)}\right), $

(1)

其中：u(x, t)表示空间x处、t时刻的种群密度；D>0为扩散系数；Δ表示关于空间变量x的Laplace算子；r>0为内禀增长率；P(x)为空间依赖的扩散策略；K(x)为空间依赖的环境容纳量。P(x)、K(x)均为取正值的非常数函数。文献[3-6]对这类考虑人均可用资源分布的单物种增长扩散或两物种竞争-共存扩散模式，进行了定性分析与数值模拟研究。

近20年来，分数阶微积分与分数阶微分方程研究受到广泛关注，尤其是在特殊环境中具有记忆特性的物质会表现出反常扩散的行为，此时应用分数阶扩散模型来描述可能更加符合实际情况^[7-8]。在生存资源总量有限的约束下，考虑种群增长扩散的历史因素与记忆性，引入种群密度函数关于时间的分数阶导数来表示种群的一种次增长率，可得到一类分数阶Logistic非线性扩散模型

$ R \partial_{t}^{\alpha} u=D \Delta\left(\frac{u(x, t)}{P(x)}\right)+r u(x, t)\left(1-\frac{u(x, t)}{K(x)}\right), $

(2)

其中：R>0表示时间尺度因子；∂_t^αu表示u(x, t)关于时间t>0的α阶Caputo分数阶导数，定义为^[7-8]

$ \partial_{t}^{\alpha} u=\frac{1}{\Gamma(1-\alpha)} \int_{0}^{t} \frac{\partial u(x, s)}{\partial s} \frac{\mathrm{d} s}{(t-s)^{\alpha}}, $

其中：α∈(0, 1)为微分阶数；Γ(·)为Gamma函数；其他各量与式(1)中的意义相同。

非线性的时间分数阶Logistic模型可以描述资源量有限但空间分布不均匀下单物种的次扩散过程及规律。对于一般非线性分数阶扩散方程的定解问题，包括整数阶的非线性扩散问题，已有的研究多是关于数值解法的^[9-14]。对于非线性分数阶或整数阶扩散方程的数值求解，主要是有限差分或有限元方法，且对于所建立的非线性格式，大多是基于能量估计，应用离散Gronwall不等式来证明差分解的稳定性与最优误差估计。本文基于模型(2)利用变量替换及分数阶导数的离散建立有限差分格式，在容纳量适当大的条件下，应用估计系数矩阵谱半径的方法证明格式的稳定性与收敛性。这种证明方法对于非线性分数阶方程差分法的理论分析尚不多见。

另一方面，对于分数阶Logistic模型(2)，环境容纳量函数K(x)是一个关键量，虽然可以估计其最大值，但一般难以获得准确的分布值，这导致相应的反问题研究。虽然关于时间分数阶扩散方程的反问题有不少研究^[15-17]，但关于非线性分数阶微分方程相关反问题的研究很少。文献[18]较早研究了时间分数阶扩散中非线性源项反演的一种条件唯一性，文献[19]针对一类具有非线性边值的时间分数阶扩散方程证明了微分阶数反演的一种Lipschitz稳定性，文献[20]应用变分型正则化方法研究了时间分数阶二维非线性扩散中的边界反演问题。对于扩散方程的参数反演估计研究，从最佳摄动量算法^[21]、贝叶斯方法^[22]到粒子群优化算法^[23]等，一直是数学物理反问题研究的重点。

本文主要探讨模型(2)在一维空间的差分解与容纳量函数的反演估计。与已有研究不同的是，本文在正问题方面对非线性差分格式的理论分析方法相对简便，主要是利用容纳量函数与Logistic源项的先验性质得到系数矩阵的谱半径估计，进而得到差分解的稳定性与收敛性。在参数反演方面，由于原方程及边界条件的复杂性，需要经过变量替换进行求解，但变换后容纳量函数不仅作为方程的系数，而且出现在初始条件中。因而，反问题转化为关于系数与初值的联合反演问题。对此，本文将应用同伦正则化算法^[24-25]进行数值反演。这种反演方法适合多个不同类型参数的联合反演，其主要思想是基于正问题的数值解，通过最优化观测数据，输入、输出数据结合正则化策略，得到稳定的迭代反演格式及最优反演解。

1 正问题与差分格式 1.1 正问题

设L>0, T>0, 记Ω=(0, L), Ω_T=Ω×(0, T), 在Ω_T上考虑一维时间分数阶非线性Logistic次增长扩散。设R=1，且P(x)=K(x)$\not \equiv$constant, 则方程(2)变为

$ \begin{aligned} & \partial_{t}^{\alpha} u(x, t)=D \frac{\partial^{2}}{\partial x^{2}}(u(x, t) / K(x))+ \\ & r u(x, t)(1-u(x, t) / K(x)) 。\end{aligned} $

(3)

给定初边值条件

$ u(x, 0)=u_{0}(x), x \in \varOmega, $

(4)

$ \begin{aligned} & \left.\frac{\partial(u(x, t) / K(x))}{\partial x}\right|_{x=0}= \\ & \left.\frac{\partial(u(x, t) / K(x))}{\partial x}\right|_{x=L}= \\ & 0, 0<t \leqslant T, \end{aligned} $

(5)

则由式(3)~(5)构成求解u=u(x, t)的定解问题。假设环境容纳量函数K(x)满足

$ K(x) \in C(\bar{\varOmega}), K(x)>1, \forall x \in \bar{\varOmega}, $

(6)

该条件说明环境容纳量在空间区域是连续分布的，且其最小值大于1，这是符合实际情况的自然条件。

由于方程(3)的右端非线性项是Logistic型的，其满足Lipschitz条件，则在初值分布非负连续条件下，正问题(3)~(5)存在唯一的非负有界解。由于方程及边界条件的复杂性，下面通过变量替换给出该问题的有限差分解。

1.2 差分格式的建立

注意到模型(3)~(5)中u(x, t)/K(x)的一致性，令v(x, t)=u(x, t)/K(x), 则方程(3)化为

$ \begin{aligned} & K(x) \partial_{t}^{\alpha} v(x, t)=D v_{x x}(x, t)+ \\ & r K(x) v(x, t)(1-v(x, t)), \end{aligned} $

(7)

相应的初边值条件变为

$ v(x, 0)=\frac{u_{0}(x)}{K(x)}, x \in \varOmega, $

(8)

$ \left.v_{x}\right|_{x=0}=\left.v_{x}\right|_{x=L}=0, 0<t \leqslant T, $

(9)

下面给出求解正问题(7)~(9)的差分格式。根据问题的物理属性，存在常数C>0, 使得0≤v(x, t)≤C成立, v(x, t)$\not \equiv$0。此外，为叙述方便，记方程(7)的右端非线性项为

$ f(v):=r K(x) v(x, t)(1-v(x, t)) 。$

对于M, N>0, 区域剖分为网格Ω_h={(x_i, t_n)x_i=ih, i=0, 1, …, M; t_n=nτ, n=0, 1, …, N}, 其中：空间步长h=L/M；时间步长τ=T/N。则在(x_i, t_n+1)处方程(7)为

$ K\left(x_{i}\right) \partial_{t}^{\alpha} v\left(x_{i}, t_{n+1}\right)=D v_{x x}\left(x_{i}, t_{n+1}\right)+f\left(x_{i}, t_{n+1}\right) 。$

(10)

根据Caputo分数阶导数的定义，直接离散时间分数阶导数为

$ \begin{aligned} & \partial_{t}^{\alpha} v\left(x_{i}, t_{n+1}\right)=\frac{1}{\Gamma(1-\alpha)} \int_{0}^{t_{n+1}}\left(t_{n+1}-\eta\right)^{-\alpha} \frac{\partial v\left(x_{i}, \eta\right)}{\partial \eta} \mathrm{d} \eta= \\ & \frac{\tau^{-\alpha}}{\Gamma(2-\alpha)} \sum\limits_{k=0}^{n} d_{k}\left[v\left(x_{i}, t_{n-k+1}\right)-v\left(x_{i}, t_{n-k}\right)\right]+O(\tau), \end{aligned} $

(11)

其中：d_k=(k+1)^1-α-k^1-α, k=0, 1, …, n。基于通常的中心差分，二阶导数离散为

$ \begin{aligned} & v_{x x}\left(x_{i}, t_{n+1}\right)=\frac{1}{h^{2}}\left(v\left(x_{i-1}, t_{n+1}\right)-2 v\left(x_{i}, t_{n+1}\right)+\right. \\ & \left.v\left(x_{i+1}, t_{n+1}\right)\right)+O\left(h^{2}\right):=\delta_{x}^{2} v\left(x_{i}, t_{n+1}\right)+O\left(h^{2}\right) 。\end{aligned} $

(12)

对于非线性项f(v), 利用泰勒展开，有

$ \begin{aligned} & f\left(v\left(x_{i}, t_{n+1}\right)\right)=2 f\left(v\left(x_{i}, t_{n}\right)\right)- \\ & f\left(v\left(x_{i}, t_{n-1}\right)\right)+O\left(\tau^{2}\right) 。\end{aligned} $

(13)

将式(11)~(13)代入式(10)，利用边值条件v₀^j=v₁^j, v_M^j=v_M-1^j, j=1, 2, …, N, 并记v_i^j≈v(x_i, t_j), v_i⁰=v₀(x_i), K_i=K(x_i), i=0, 1, …, M, 可得差分格式：

当n=0时，

$ K_{i} \frac{\tau^{-\alpha}}{\Gamma(2-\alpha)}\left(v_{i}^{1}-v_{i}^{0}\right)-D \delta_{x}^{2} v_{i}^{1}=f\left(v_{i}^{0}\right); $

(14)

当n≥1时，

$ \begin{aligned} & K_{i} \frac{\tau^{-\alpha}}{\Gamma(2-\alpha)} \sum\limits_{k=0}^{n} d_{k}\left(v_{i}^{n-k+1}-v_{i}^{n-k}\right)-D \delta_{x}^{2} v_{i}^{n+1}= \\ & 2 f\left(v_{i}^{n}\right)-f\left(v_{i}^{n-1}\right)。\end{aligned} $

(15)

再记

$ \begin{gathered} p=\frac{D \tau^{\alpha} \Gamma(2-\alpha)}{h^{2}}, q=\tau^{\alpha} \Gamma(2-\alpha), \\ \boldsymbol{K}^{*}=\left(K_{1}, K_{2}, \cdots, K_{M-1}\right)^{\mathrm{T}}, \\ \boldsymbol{V}^{n}=\left(v_{1}^{n}, v_{2}^{n}, \cdots, v_{M-1}^{n}\right)^{\mathrm{T}}, \end{gathered} $

$ \begin{gathered} \boldsymbol{W}^{n}=\left(\left(v_{1}^{n}\right)^{2}, \left(v_{2}^{n}\right)^{2}, \cdots, \left(v_{M-1}^{n}\right)^{2}\right)^{\mathrm{T}}, \\ \boldsymbol{H}^{n}:=\boldsymbol{K}^{*} \cdot \boldsymbol{V}^{n}=\left(K_{1} v_{1}^{n}, K_{2} v_{2}^{n}, \cdots, K_{M-1} v_{M-1}^{n}\right)^{\mathrm{T}} , \\ \boldsymbol{I}^{n}:=\boldsymbol{K}^{*} \cdot \boldsymbol{W}^{n}=\left(K_{1}\left(v_{1}^{n}\right)^{2}, K_{2}\left(v_{2}^{n}\right)^{2}, \cdots, K_{M-1}\left(v_{M-1}^{n}\right)^{2}\right)^{\mathrm{T}} 。\end{gathered} $

注意到f(v)=rK(x)v(x, t)(1-v(x, t)), 经过整理可得矩阵形式的差分格式

$ \left\{\begin{aligned} \boldsymbol{A} \boldsymbol{V}^{1}= & (1+q r) \boldsymbol{H}^{0}-q r \boldsymbol{I}^{0}, \\ \boldsymbol{A} \boldsymbol{V}^{n+1}= & \sum\limits_{j=1}^{n}\left(d_{j-1}-d_{j}\right) \boldsymbol{H}^{n+1-j}+d_{n} \boldsymbol{H}^{0}+2 q r \boldsymbol{H}^{n}- \\ & 2 q r \boldsymbol{I}^{n}-q r \boldsymbol{H}^{n-1}+q r \boldsymbol{I}^{n-1}, \end{aligned}\right. $

(16)

其中：n=1, 2, …, N-1；系数矩阵

$ \boldsymbol{A}=\left(\begin{array}{cccccc} K_{1}+p & -p & 0 & \cdots & 0 & 0 \\ -p & K_{2}+2 p & -p & \cdots & 0 & 0 \\ 0 & -p & K_{3}+2 p & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & K_{M-2}+2 p & -p \\ 0 & 0 & 0 & \cdots & -p & K_{M-1}+p \end{array}\right)。$

(17)

注意到K_i>1，i=1, 2，…, M-1，p>0, 容易验证由(17)式定义的系数矩阵A严格对角占优，即知差分格式(16)唯一可解。对于系数矩阵A，仍有条件K_i>1，i=1, 2，…, M-1，以及不等式

$ \begin{aligned} & a_{i i}=K_{i}+\sum\limits_{j=1, j \neq i}^{M-1}\left|a_{i j}\right|>1+\sum\limits_{j=1, j \neq i}^{M-1}\left|a_{i j}\right|, \\ & i=1, 2, \cdots, M-1, \end{aligned} $

(18)

可知系数矩阵A的谱半径是严格大于1的。因而, 可应用谱估计的方法，证明差分格式(16)的稳定性与收敛性。

2 差分格式的稳定性与收敛性 2.1 稳定性

在下面的叙述中，C及C_j(j=1, 2, …)表示任意正常数。

定理1   设容纳量函数K(x)满足条件(6)，则差分格式(16)是无条件稳定的。

证明   对于差分格式(16)，设λ为系数矩阵A的任一特征值，则由Gerschgorin圆盘定理，可知

$ \left|\lambda-a_{i i}\right| \leqslant \sum\limits_{s=1, s \neq i}^{M-1}\left|a_{i s}\right|, i=1, 2, \cdots, M-1 。$

由条件(6)及式(18)，有

$ |\lambda| \geqslant a_{i i}-\sum\limits_{s=1, s \neq i}^{M-1}\left|a_{i s}\right|>1, i=1, 2, \cdots, M-1 。$

因而必有

$ \rho(\boldsymbol{A})=\max |\lambda|>1, $

得ρ(A^-1)＜1。根据矩阵范数的性质，可知存在A^-1的某种范数‖·‖, 使得‖A^-1‖＜1成立。

记初始数据的扰动为$\boldsymbol{E}^{0}=\tilde{\boldsymbol{V}}^{0}-\boldsymbol{V}^{0}$, 第n层的扰动为$\boldsymbol{E}^{n}=\tilde{\boldsymbol{V}}^{n}-\boldsymbol{V}^{n}(n \in \mathbf{N})$, 并记$\left\|\boldsymbol{K}^{*}\right\|_{\infty}= \max\limits _{x \in \varOmega}|K(x)|$。由差分格式(16)，可得

$ \left\{\begin{aligned} \boldsymbol{A} \boldsymbol{E}^{1}= & (1+q r)\left(\tilde{\boldsymbol{H}}^{0}-\boldsymbol{H}^{0}\right)+q r\left(\boldsymbol{I}^{0}-\tilde{\boldsymbol{I}}^{0}\right), \\ \boldsymbol{A} \boldsymbol{E}^{n+1}= & \sum\limits_{j=1}^{n}\left(d_{j-1}-d_{j}\right)\left(\tilde{\boldsymbol{H}}^{n+1-j}-\boldsymbol{H}^{n+1-j}\right)+ \\ & d_{n}\left(\tilde{\boldsymbol{H}}^{0}-\boldsymbol{H}^{0}\right)+2 q r\left(\tilde{\boldsymbol{H}}^{n}-\boldsymbol{H}^{n}\right)- \\ & 2 q r\left(\tilde{\boldsymbol{I}}^{n}-\boldsymbol{I}^{n}\right)-q r\left(\tilde{\boldsymbol{H}}^{n-1}-\boldsymbol{H}^{n-1}\right)+ \\ & q r\left(\tilde{\boldsymbol{I}}^{n-1}-\boldsymbol{I}^{n-1}\right), n \in \mathbf{N}。\end{aligned}\right. $

(19)

由式(19)的第1式，利用‖A^-1‖＜1及$\left\|\tilde{\boldsymbol{H}}^{0}-\boldsymbol{H}^{0}\right\| \leqslant\left\|\boldsymbol{K}^{*}\right\|_{\infty}\left\|\boldsymbol{E}^{0}\right\|$可得

$ \begin{aligned} & \left\|\boldsymbol{E}^{1}\right\| \leqslant(1+q r)\left\|\boldsymbol{A}^{-1}\right\|\left\|\tilde{\boldsymbol{H}}^{0}-\boldsymbol{H}^{0}\right\|+ \\ & q r\left\|\boldsymbol{A}^{-1}\right\|\left\|\tilde{\boldsymbol{I}}^{0}-\boldsymbol{I}^{0}\right\| \leqslant \\ & (1+q r)\left\|\boldsymbol{K}^{*}\right\|_{\infty}\left\|\boldsymbol{E}^{0}\right\|+q r\left\|\tilde{\boldsymbol{I}}^{0}-\boldsymbol{I}^{0}\right\| 。\end{aligned} $

(20)

注意到解的有界性，即有$\left\|\tilde{\boldsymbol{V}}^{n}\right\|, \left\|\boldsymbol{V}^{n}\right\| \leqslant C, n \in \mathbf{N} $。

因而

$ \begin{aligned} & \left\|\tilde{\boldsymbol{I}}^{0}-\boldsymbol{I}^{0}\right\| \leqslant\left\|\boldsymbol{K}^{*}\right\|_{\infty}\left\|\left(\tilde{\boldsymbol{V}}^{0}\right)^{2}-\left(\boldsymbol{V}^{0}\right)^{2}\right\| \leqslant \\ & \left\|\boldsymbol{K}^{*}\right\|_{\infty} \cdot 2 C\left\|\tilde{\boldsymbol{V}}^{0}-\boldsymbol{V}^{0}\right\|=C_{1}\left\|\boldsymbol{E}^{0}\right\|, \end{aligned} $

(21)

其中$C_{1}=2 C\left\|\boldsymbol{K}^{*}\right\|_{\infty}$。结合式(20)，有

$ \begin{aligned} & \left\|\boldsymbol{E}^{1}\right\| \leqslant(1+q r)\left\|\boldsymbol{K}^{*}\right\|_{\infty}\left\|\boldsymbol{E}^{0}\right\|+ \\ & q r C_{1}\left\|\boldsymbol{E}^{0}\right\| \leqslant C_{2}\left\|\boldsymbol{E}^{0}\right\|, \end{aligned} $

(22)

其中C₂=max{(1+qr)‖K^*‖_∞, qrC₁}。

假设对于j=2, 3，…, n, ‖E^j‖≤C‖E⁰‖成立, C为与n无关的正常数。j=1, 2, …, n时，由H^j、I^j与V^j的关系，可知

$ \begin{aligned} & \left\|\tilde{\boldsymbol{H}}^{j}-\boldsymbol{H}^{j}\right\| \leqslant\left\|\boldsymbol{K}^{*}\right\|_{\infty}\left\|\boldsymbol{E}^{j}\right\|, \left\|\tilde{\boldsymbol{I}}^{j}-\boldsymbol{I}^{j}\right\| \leqslant \\ & C_{1}\left\|\boldsymbol{E}^{j}\right\|, j=1, 2, \cdots, n。\end{aligned} $

(23)

由于$d_{j-1}-d_{j}>0(j=1, 2, \cdots, n)$，及$\sum\limits_{j=1}^{n}\left(d_{j-1}-d_{j}\right)+ d_{n}=d_{0}=1$, 由式(19)的第2式，利用‖A^-1‖＜1, 可知j=n+1时，有

$ \begin{aligned} & \left\|\boldsymbol{E}^{n+1}\right\| \leqslant\left\|\boldsymbol{K}^*\right\|_{\infty}\left\{\sum_{j=1}^n\left(d_{j-1}-d_j\right)\left\|\boldsymbol{E}^{n+1-j}\right\|\right\}+ \\ & \left\|\boldsymbol{K}^*\right\|_{\infty} d_n\left\|\boldsymbol{E}^0\right\|+2 q r\left\|\boldsymbol{K}^*\right\|_{\infty}\left\|\boldsymbol{E}^n\right\|+ \\ & 2 q r C_1\left\|\boldsymbol{E}^n\right\|+q r\left\|\boldsymbol{K}^*\right\|_{\infty}\left\|\boldsymbol{E}^{n-1}\right\|+q r C_1\left\|\boldsymbol{E}^{n-1}\right\| \leqslant \\ & C\left\|\boldsymbol{K}^*\right\|_{\infty}\left\|\boldsymbol{E}^0\right\|+\left(2 q r\left\|\boldsymbol{K}^*\right\|_{\infty}+2 q r C_1+\right.\\ & \left.q r\left\|\boldsymbol{K}^*\right\|_{\infty}+q r C_1\right)\left\|\boldsymbol{E}^0\right\|=C_3\left\|\boldsymbol{E}^0\right\|, \end{aligned} $

(24)

其中C₃>0为依赖于q、r、‖K^*‖_∞，但与n无关的正常数。由归纳法即知差分格式(16)是稳定的，证毕。

2.2 收敛性

定理2   对于任意给定的T∈(0, +∞), 差分格式(16)收敛，且当h, τ→0时，有

$ \left\|\boldsymbol{e}^{n}\right\| \leqslant \frac{C_{4} T^{\alpha}}{1-\alpha}\left(\tau+\tau^{2}+h^{2}\right)=O\left(\tau+h^{2}\right), $

(25)

其中：eⁿ=V^{exa, n}-Vⁿ(n∈N)表示第n层的解误差，V^{exa, n}表示第n层的精确解向量；C₄为与n无关的正常数。

证明   由差分格式(16)，注意到e⁰=0及H^{exa, 0}-H⁰=0, I^{exa, 0}-I⁰=0, 得

$ \left\{\begin{aligned} \boldsymbol{A} \boldsymbol{e}^{1}= & \boldsymbol{R}^{1}, \\ \boldsymbol{A} \boldsymbol{e}^{n+1}= & \sum\limits_{j=1}^{n}\left(d_{j-1}-d_{j}\right)\left(\boldsymbol{H}^{\text {exa }, n+1-j}-\boldsymbol{H}^{n+1-j}\right)+ \\ & 2 q r\left(\boldsymbol{H}^{\text {exa }, n}-\boldsymbol{H}^{n}\right)-2 q r\left(\boldsymbol{I}^{\text {exa }, n}-\boldsymbol{I}^{n}\right)- \\ & q r\left(\boldsymbol{H}^{\text {exa }, n-1}-\boldsymbol{H}^{n-1}\right)+q r\left(\boldsymbol{I}^{\text {exa }, n-1}-\boldsymbol{I}^{n-1}\right)+ \\ & \boldsymbol{R}^{n+1}, n \in \mathbf{N}。\end{aligned}\right. $

(26)

其中Rⁿ表示第n层的截断误差。由式(11)~(13)可知，存在与n无关的常数C>0,

$ \left\|\boldsymbol{R}^{n}\right\| \leqslant C \tau^{\alpha}\left(\tau+\tau^{2}+h^{2}\right), n \in \mathbf{N}。$

(27)

利用‖A^-1‖＜1及式(27)的估计，由式(26)的第1式，可得

$ \left\|\boldsymbol{e}^{1}\right\| \leqslant\left\|\boldsymbol{A}^{-1} \boldsymbol{R}^{1}\right\|<\left\|\boldsymbol{R}^{1}\right\| \leqslant C \tau^{\alpha}\left(\tau+\tau^{2}+h^{2}\right) 。$

(28)

注意到解的有界性，类似于式(23)的估计，可知对于j=1, 2，…, n,

$ \begin{gathered} \left\|\boldsymbol{H}^{\text {exa }, j}-\boldsymbol{H}^{j}\right\| \leqslant\left\|\boldsymbol{K}^{*}\right\|_{\infty}\left\|\boldsymbol{e}^{j}\right\|, \\ \left\|\boldsymbol{I}^{\text {exa }, j}-\boldsymbol{I}^{j}\right\| \leqslant C_{1}\left\|\boldsymbol{K}^{*}\right\|_{\infty}\left\|\boldsymbol{e}^{j}\right\| 。\end{gathered} $

(29)

假设对于j=2, 3，…, n,

$ \left\|\boldsymbol{e}^{j}\right\| \leqslant C \tau^{\alpha}\left(\tau+\tau^{2}+h^{2}\right), $

(30)

则当j=n+1时，由式(26)的第2式及式(29)，再次利用‖A^-1‖＜1, 并注意到d_n＜d_n-1＜…＜d₀=1, 及式(30)，得

$ \begin{aligned} & \left\|\boldsymbol{e}^{n+1}\right\| \leqslant\left\|\boldsymbol{K}^{*}\right\|_{\infty} \sum\limits_{j=1}^{n}\left|d_{j-1}-d_{j}\right|\left\|\boldsymbol{e}^{n+1-j}\right\|+ \\ & 2 q r\left\|\boldsymbol{K}^{*}\right\|_{\infty}\left\|\boldsymbol{e}^{n}\right\|+2 q r C_{1}\left\|\boldsymbol{K}^{*}\right\|_{\infty}\left\|\boldsymbol{e}^{n}\right\|+ \\ & q r\left\|\boldsymbol{K}^{*}\right\|_{\infty}\left\|\boldsymbol{e}^{n-1}\right\|+q r C_{1}\left\|\boldsymbol{K}^{*}\right\|_{\infty}\left\|\boldsymbol{e}^{n-1}\right\|+ \\ & \left\|\boldsymbol{R}^{n+1}\right\| \leqslant C \tau^{\alpha}\left(\tau+\tau^{2}+h^{2}\right)\left(\left\|\boldsymbol{K}^{*}\right\|_{\infty}+\right. \\ & \left.3 q r\left\|\boldsymbol{K}^{*}\right\|_{\infty}+3 C_{1} q r\left\|\boldsymbol{K}^{*}\right\|_{\infty}+1\right) \leqslant\\ & C_{4}\left(\tau+\tau^{2}+h^{2}\right) \frac{\tau^{\alpha}}{d_{n}}, \end{aligned} $

(31)

其中C₄为不依赖于n的正常数。

进一步，由于nτ≤T, 得

$ \left\|\boldsymbol{e}^{n+1}\right\| \leqslant C_4\left(\tau+\tau^2+h^2\right) T^\alpha \frac{1}{n^\alpha d_n} \text { 。} $

注意到d_n=(n+1)^1-α-n^1-α, 以及

$ \lim\limits _{n \rightarrow \infty} \frac{1}{n^{\alpha}\left[(n+1)^{1-\alpha}-n^{1-\alpha}\right]}=\frac{1}{1-\alpha}, $

得

$ \left\|\boldsymbol{e}^{n+1}\right\| \leqslant \frac{C_{4} T^{\alpha}}{1-\alpha}\left(\tau+\tau^{2}+h^{2}\right) 。$

(32)

由式(28)、(30)与(32)，根据归纳法即知∀n∈N, 估计式(25)成立，差分格式(16)收敛。定理证毕。

2.3 数值试验

算例1   考虑带强制项的方程

$ \begin{aligned} & \partial_{t}^{\alpha} u(x, t)=D \frac{\partial^{2}}{\partial x^{2}}(u(x, t) / K(x))+ \\ & r u(x, t)(1-u(x, t) / K(x))+\gamma(x, t), \end{aligned} $

其中γ(x, t)为强制项。不妨设微分阶数α=0.3, 其他模型参数D=0.01, r=0.001, 环境容纳量K(x)=120+e^-x。设正问题的真解为

$ u(x, t)=100 x^{2}(1-x)^{2}\left(1+t^{2 \alpha}\right), $

则强制项为

$ \begin{aligned} & \gamma(x, t)=\frac{100 x^{2}(1-x)^{2} \Gamma(2 \alpha+1)}{\Gamma(\alpha+1)} t^{\alpha}- \\ & \frac{\partial^{2}}{\partial x^{2}}\left(\frac{x^{2}(1-x)^{2}}{120+e^{-x}}\right)\left(1+t^{2 \alpha}\right)-0.1 x^{2}(1-x)^{2} \\ & \left(1-\frac{100 x^{2}(1-x)^{2}\left(1+t^{2 \alpha}\right)}{120+e^{-x}}\right)\left(1+t^{2 \alpha}\right) 。\end{aligned} $

应用差分格式(16)，取网格比τ=h²进行计算，计算结果如表 1所示，其中Err表示真解与数值解在t=0.5时的相对误差，定义为Err=‖u(x, 0.5)-u^*(x, 0.5)‖₂/‖u(x, 0.5)‖₂, 这里u^*(x, 0.5)表示t=0.5时的差分解，Rat表示收敛速率。

从表 1的结果可以看出，随着网格加密，数值解与真解的误差逐步减小。图 1绘制了网格逐渐加密(τ=h²)时的真解与数值解的图像。可以看到，随着网格加密，数值解收敛到真解。

3 反问题与数值反演 3.1 确定环境容纳量的反问题

方程(3)中的环境容纳量K(x)是一个重要的分布参数，利用适当的观测数据对此参数进行数值反演是本节的主要工作。如果考虑基于式(3)~(5)确定K(x)的反演问题，不仅方程是复杂的，而且边界条件也含有K(x), 目前看并没有合适的反演方法。值得注意的是，经过变量替换后，得到的方程(7)及初边值条件(8)~(9)是相对简便的模型。虽然方程(7)中仍然含有K(x), 但边界条件简化了, 且K(x)出现在初始条件中，这相比变换前的模型，至少从数值方法上容易实现对K(x)的反演。

给定空间内某一点处的观测数据作为附加数据，有

$ v\left(x_{0}, t\right)=\theta(t), 0<t \leqslant T, $

(33)

其中：x₀∈Ω为固定的观测点。实际问题中，测量数据总带有误差。设θ^ε(t)为相应的扰动数据，满足条件

$ \left\|\theta^{\varepsilon}(t)-\theta(t)\right\|_{L^{2}(0, T)} \leqslant \varepsilon, $

这里ε>0称为扰动水平。由式(7)~(9)和扰动数据θ^ε(t), 形成了数据扰动条件下确定环境容纳量K(x)的反问题。对于这类反演问题的数值求解, 同伦正则化算法是一种有效方法^[24-25]。

3.2 数值反演

本节计算中，同伦参数为

$ \lambda_{n}=\frac{1}{1+\exp (0.8(n-5))}, $

其中n是迭代次数。此外，若无特别说明，计算区域Ω_T=(0, 1)×(0, 1), 附加数据在x₀=0.5处取得。正问题求解中, 初值u(x, 0)=x²(1-x)², 空间离散水平M=100, 时间离散水平N=100。具体反演模拟中，利用环境容纳量的真解K(x)得到附加数据(33)，进而在多项式逼近空间Φ^S=span{1, x, …, x^S-1}中, 应用同伦正则化算法对K(x)进行反演重建。

算例2   首先考察逼近空间的维数对反演算法的影响。设容纳量的真解K(x)=10+exp(-x), 微分阶数α=0.5, 扩散系数D=0.1, 增长率r=0.01。根据exp(-x)在多项式逼近空间中的展开式，真解K(x)在空间Φ^S中表示为

$ \boldsymbol{a}=\left(11, -1, 1 / 2, -1 / 6, \cdots, (-1)^{S-1} \frac{1}{(S-1) !}\right) 。$

由于截断误差的存在, 理论上逼近空间维数越高, 反演结果应该越精确。应用不带扰动的数据, 分别在逼近空间维数S=3, 4, 5, 6时，对K(x)进行反演重建。初始迭代取a₀=(1, 0), 其中0表示相应维数的零向量。在不同维数逼近空间中的反演解与真解的图像如图 2所示。

由图 2看出, 当逼近空间维数增大时, 反演解收敛于真解, 说明了反演算法有限维逼近收敛性。

算例3   考察数据扰动对反演算法的影响。仍取容纳量K(x)=10+exp(-x), 微分阶数、扩散系数等参数同于例2，在逼近空间span{1, x, x², x³, x⁴}中进行反演计算。此时, 容纳量真解表示为a=(11, -1, 1/2, -1/6, 1/24)。取初始迭代a₀=(1, 0), 不同扰动水平下的平均反演结果列于表 2, 其中ε表示扰动水平, a^inv表示连续10次反演解的平均值，Err表示平均反演解与真解的误差。

从表 2的计算结果可以看出，随着扰动水平的降低，反演解与真解的误差逐步减小，说明了反演算法的一种数值稳定性。

算例4   设容纳量的真解为K(x)=20+sin(x), 扩散系数、增长率、ε、a^inv和Err同算例3，微分阶数α=0.75。根据sin(x)的展开式, K(x)在逼近空间中近似表示为a=(20, 1, 0, -1/6, 0)。仍取初始迭代a₀=(1, 0), 不同扰动水平下的平均反演结果列于表 3。

从表 3可以看出, 完全类似于算例3的结果, 随着扰动水平的降低, 反演解与真解的误差逐步减小, 反演解收敛到真解。

4 结论

本文主要研究了一个非线性分数阶Logistic增长扩散模型的差分解与空间依赖容纳量函数的反演重建问题。在环境容纳量适当大的先验条件下，利用谱半径估计方法证明了差分格式的稳定性与收敛性。对于重建容纳量函数的反问题，利用一个内点的观测数据，应用同伦正则化算法进行了有效的数值反演。下一步，对于这类时间分数阶Logistic次增长扩散模型，我们将考虑基于终值观测数据确定容纳量函数的反问题，同时开展容纳量反演重建的理论研究等。