重心插值Newton-Raphson迭代法求解非线性分数阶Sobolev方程

Barycentric interpolation Newton-Raphson iterative method for solving nonlinear fractional Sobolev equation

摘要

为求解非线性分数阶Sobolev方程，采用了重心插值Newton-Raphson迭代法，构造了方程的Newton-Raphson迭代格式。对于方程中的分数阶导数项，采用分部积分法，克服Caputo分数阶导数的奇异性，使之转化为黎曼积分，并利用高斯正交公式近似计算分数阶导数，通过重心插值配点法离散方程，根据微分矩阵得到方程的矩阵方程，进而得到Jacobi矩阵，构造出非线性分数阶Sobolev方程的Newton-Raphson迭代格式。最后，通过三个非线性项不同的非线性分数阶Sobolev方程的数值算例结果表明，该方法对于求解非线性分数阶Sobolev方程的有效性。

Abstract

In order to solve the nonlinear fractional Sobolev equation, the barycentric interpolation Newton-Raphson iterative method is used to construct the Newton-Raphson iterative scheme of the equation. For the fractional derivative term in the equation, the integration by parts is used to overcome the singularity of the Caputo fractional derivative and transform it into Riemann integral. The fractional derivative is approximately calculated by Gauss quadrature formula, the barycentric interpolation collocation method is used to discretize the equation, and the matrix equation of the equation is obtained according to the differential matrix, and then the Jacobi matrix is obtained. The Newton-Raphson iterative scheme of the nonlinear fractional Sobolev equation is constructed. At last, the numerical results of three nonlinear fractional Sobolev equations with different nonlinear terms show that the method is effective for solving nonlinear fractional Sobolev equations.

Graphical abstract

关键词

重心插值Newton-Raphson迭代法 / 非线性分数阶Sobolev方程 / 矩阵方程 / Caputo分数阶导数

Key words

barycentric interpolation Newton-Raphson iterative method / nonlinear fractional Sobolev equation / matrix equation / Caputo fractional derivative

非线性分数阶偏微分方程在生物种群的扩散、流体流动、电磁波、动力系统的控制理论等^[1-3]应用科学中起着重要的作用。非线性分数阶Sobolev方程作为一类重要的非线性分数阶偏微分方程，包含了对空间和时间的三阶混合偏导数项，具有分数阶导数项和非线性项，可用于解决流体通过裂隙岩石的流动、黏土固结、波在介质中的运动以及湿度通过土壤的运动等问题^[4-7]。目前求解非线性分数阶Sobolev方程的研究较少，已有的数值方法有有限差分法^[8-9]、Crank-Nicolson方法^[10]、谱法^[11]等。

重心插值配点法是一种不需要划分计算网络、计算程序易于编写、数值稳定性好、精度高的配点型方法。重心插值在工程和物理方面应用广泛，如板的弯曲、振动频率等问题，也可应用于求解方程，如Kolmogorov-Petrovskii-Piskunov方程^[12]、双调和方程^[13]、Volterra积分方程^[14]、非线性热传导方程^[15]、电报方程^[16]、梁力振动方程^[17]等。重心插值Newton-Raphson迭代法是利用重心插值近似计算未知函数，将插值函数代入非线性微分方程并令微分方程在计算节点处精确成立，从而得到非线性代数方程组，再采用Newton-Raphson迭代法求解，得到非线性微分方程数值解的方法。重心插值Newton-Raphson迭代法已被广泛应用于非线性微分方程数值求解的研究中，如Sine-Gordon方程^[18]、非线性椭圆型方程^[19]、非线性梁方程^[20]、非线性扩散方程等。

本文的结构安排如下：第1节介绍了(1+1)维和(2+1)维非线性分数阶Sobolev方程。第2节介绍了重心型插值及Newton-Raphson迭代法的相关理论。第3节给出了非线性分数阶Sobolev方程矩阵方程，并构造出了方程的Newton-Raphson迭代格式。第4节给出了三个数值算例，对方法求解非线性分数阶Sobolev方程的有效性进行了验证。第5节给出了一些结论。

1 非线性分数阶Sobolev方程

本文研究以下形式的(1+1)维和(2+1)维非线性分数阶Sobolev方程。(1+1)维非线性分数阶Sobolev方程为

\left\{\begin{array}{l}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}-k_{\mathrm{1}}\frac{\partial \mathrm{\Delta }u\left(x,t\right)}{\partial t}-k_{\mathrm{2}}\mathrm{\Delta }u\left(x,t\right)\\ =G\left(u\right)+f\left(x,t\right),\left(x,t\right)\in \left[a,b\right]\times \left[\mathrm{0},T\right]\\ u\left(x,\mathrm{0}\right)=u_{\mathrm{0}}\left(x\right)\\ u\left(a,t\right)={\psi }_{\mathrm{1}}\left(t\right),u\left(b,t\right)={\psi }_{\mathrm{2}}\left(t\right)\end{array}\right.(1)

(2+1)维非线性分数阶Sobolev方程为

\left\{\begin{array}{l}\frac{{\partial }^{\alpha }u\left(x,y,t\right)}{\partial t^{\alpha }}-k_{\mathrm{1}}\frac{\partial \mathrm{\Delta }u\left(x,y,t\right)}{\partial t}-k_{\mathrm{2}}\mathrm{\Delta }u\left(x,y,t\right)\\ =G\left(u\right)+f\left(x,y,t\right)\\ u\left(x,y,\mathrm{0}\right)=u_{\mathrm{0}}\left(x,y\right)\\ u\left(a,y,t\right)=q_{\mathrm{1}}\left(y,t\right),u\left(b,y,t\right)=q_{\mathrm{2}}\left(y,t\right)\\ u\left(x,c,t\right)=p_{\mathrm{1}}\left(x,t\right),u\left(x,d,t\right)=p_{\mathrm{2}}\left(x,t\right)\end{array}\right.(2)

式中：\left(x,y,t\right)\in \Omega \times \left[\mathrm{0},T\right]，\Omega =\left[a,b\right]\times \left[c,d\right]，\mathrm{\Delta }表示Laplace算子，k_{\mathrm{1}}和k_{\mathrm{2}}表示两个大于0的常数，，G，f，u_{\mathrm{0}}，{\psi }_{\mathrm{1}}，{\psi }_{\mathrm{2}}，q_{\mathrm{1}}，q_{\mathrm{2}}，p_{\mathrm{1}}，p_{\mathrm{2}}均为已知函数，非线性项G\left(u\right)满足Lipschitz条件，即

\left|G\left(u_{\mathrm{1}}\right)-G\left(u_{\mathrm{2}}\right)\right|\le L\left|u_{\mathrm{1}}-u_{\mathrm{2}}\right|

其中：u_{\mathrm{1}},u_{\mathrm{2}}\in \mathbb{C}\left(\Omega \times \left[\mathrm{0},T\right]\right)，L为大于0的Lipschitz常数。方程中的分数阶导数部分表示为Caputo分数阶导数，形式如下：

\frac{{\partial }^{\alpha }\varphi \left(\mathit{x},t\right)}{\partial t^{\alpha }}=\frac{\mathrm{1}}{\mathrm{\Gamma }\left(\mathrm{1}-\alpha \right)}{\int }_{\mathrm{0}}^t\frac{\partial \varphi \left(\mathit{x},\tau \right)}{\partial \tau }\frac{\mathrm{1}}{{\left(t-\tau \right)}^{\alpha }}\mathrm{d}\tau

其中：\mathrm{\Gamma }\left(·\right)表示伽马函数，定义为

\mathrm{\Gamma }\left(\alpha \right)={\int }_{\mathrm{0}}^{+\mathrm{\infty }}{x}^{\alpha -\mathrm{1}}{\mathrm{e}}^{-x}\mathrm{d}x,x\in \left(\mathrm{0},+\mathrm{\infty }\right)

2 重心型插值及Newton-Raphson迭代法

2.1 重心型插值

本节以定义域为\left[a,b\right]的函数u\left(x\right)为例介绍重心型插值及其微分矩阵。设函数的定义域上有\left(m+\mathrm{1}\right)个等距节点或第二类切比雪夫节点，表示为x_{\mathrm{0}},

x_{\mathrm{1}},\cdots ,x_m，其中第二类切比雪夫节点的定义为x_i=\mathrm{c}\mathrm{o}\mathrm{s}\frac{i\mathrm{\pi }}{m}，i=\mathrm{0},\mathrm{1},\cdots ,m。此时，重心插值函数为

u={\displaystyle \sum _{i=\mathrm{0}}^m}{L}_i\left(x\right)u_i(3)

其中：u_i=u\left(x_i\right)；L_i(x)表示空间变量x方向上的基函数，定义为

L_i\left(x\right)=\frac{\frac{{\omega }_i}{x-x_i}}{{\displaystyle \sum _{h=\mathrm{0}}^m}\frac{{\omega }_h}{x-x_h}}(4)

当权函数{\omega }_i的定义为

{\omega }_i=\frac{\mathrm{1}}{{\displaystyle \prod _{h=\mathrm{0},h\ne i}^m}(x_i-x_h)},i=\mathrm{0},\mathrm{1},\cdots ,m(5)

时，重心插值式(3)为重心Lagrange插值。当权函数{\omega }_i的定义为

{\omega }_i={\sum _{r\in J_i}\left(-\mathrm{1}\right)}^r{\displaystyle \prod _{h=r,r\ne i}^{r+d_x}}\frac{\mathrm{1}}{x_i-x_h}(6)

时，重心插值式(3)为重心有理插值，其中：J_i=\left\{r\in J_x:i-d_x\le r\le i\right\}，\mathrm{0}\le d_x\le m，d_x为整数，J_x=\left\{\mathrm{0},\mathrm{1},\cdots ,m-d_x\right\}。

对式(3)求x的1阶导数，可得u\mathrm{\text{'}}\left(x\right)={\displaystyle \sum _{i=\mathrm{0}}^m}{L}_i^{\mathrm{\text{'}}}\left(x\right)u_i，其中

L_i^{\mathrm{\text{'}}}\left(x_{\vartheta }\right)=\left\{\begin{array}{l}\frac{{\omega }_i/{\omega }_{\vartheta }}{x_{\vartheta }-x_i},i\ne \vartheta \\ -{\displaystyle \sum _{i\ne \vartheta }^{}}{L}_i^{\mathrm{\text{'}}}\left(x_{\vartheta }\right),i=\vartheta \end{array}\right.(7)

且\vartheta =\mathrm{0},\mathrm{1},\cdots ,m，则L_i(x)在插值节点处的导数递推公式为

L_i^{\left(k_x\right)}\left(x_{\vartheta }\right)=\left\{\begin{array}{l}{k}_x\left[L_{\vartheta }^{\left(k_x-\mathrm{1}\right)}\left(x_{\vartheta }\right)L_i^{\mathrm{\text{'}}}\left(x_{\vartheta }\right)-\frac{L_i^{\left(k_x-\mathrm{1}\right)}\left(x_{\vartheta }\right)}{x_{\vartheta }-x_i}\right],i\ne \vartheta \\ -{\displaystyle \sum _{i=\mathrm{0},i\ne \vartheta }^m}{L}_i^{\left(k_x\right)}\left(x_{\vartheta }\right),i=\vartheta \end{array}\right.(8)

此时k_x\ge \mathrm{2}。

对于给定的插值节点x_{\mathrm{0}},x_{\mathrm{1}},\cdots ,x_m，代入到式(3)中，则k_x阶导数可表示为

u^{\left(k_x\right)}(x)={\displaystyle \sum _{i=\mathrm{0}}^m}{L}_i^{\left(k_x\right)}(x)u_i(9)

式(9)可用矩阵的形式表示为

{\mathit{u}}^{\left(k_x\right)}={\mathit{D}}^{\left(k_x\right)}\mathit{u}(10)

式中：{\mathit{u}}^{\left(k_x\right)}={\left[u_{\mathrm{0}}^{\left(k_x\right)},u_{\mathrm{1}}^{\left(k_x\right)},\cdots ,u_n^{\left(k_x\right)}\right]}^{\mathrm{T}}；{\mathit{D}}^{\left(k_x\right)}的元素表示为D_{\vartheta i}^{\left(k_x\right)}=L_i^{\left(k_x\right)}\left(x_{\vartheta }\right)；\mathit{u}={\left[u_{\mathrm{0}},u_{\mathrm{1}},\cdots ,u_n\right]}^{\mathrm{T}}。

2.2 Newton-Raphson迭代法

对于非线性方程组

\left\{\begin{array}{l}{f}_{\mathrm{1}}\left(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_n\right)=\mathrm{0}\\ f_{\mathrm{2}}\left(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_n\right)=\mathrm{0}\\ ⋮\\ f_n\left(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_n\right)=\mathrm{0}\end{array}\right.(11)

其矩阵形式为

\mathit{F}\left(\mathit{x}\right)=\mathrm{0}

其中：\mathit{F}\left(\mathit{x}\right)={\left[f_{\mathrm{1}}\left(\mathit{x}\right),f_{\mathrm{2}}\left(\mathit{x}\right),\cdots ,f_n\left(\mathit{x}\right)\right]}^{\mathrm{T}}，\mathit{x}=\left[x_{\mathrm{1}},x_{\mathrm{2}},\right.{\left.\cdots ,x_n\right]}^{\mathrm{T}}，则函数\mathit{F}\left(\mathit{x}\right)的Frechet导数的Jacobi矩阵表示为

\mathit{F}\mathrm{\text{'}}\left(\mathit{x}\right)=\left[\begin{array}{cccc}\frac{\partial f_{\mathrm{1}}\left(\mathit{x}\right)}{\partial x_{\mathrm{1}}}& \frac{\partial f_{\mathrm{1}}\left(\mathit{x}\right)}{\partial x_{\mathrm{2}}}& \cdots & \frac{\partial f_{\mathrm{1}}\left(\mathit{x}\right)}{\partial x_n}\\ \frac{\partial f_{\mathrm{2}}\left(\mathit{x}\right)}{\partial x_{\mathrm{1}}}& \frac{\partial f_{\mathrm{2}}\left(\mathit{x}\right)}{\partial x_{\mathrm{2}}}& \cdots & \frac{\partial f_{\mathrm{2}}\left(\mathit{x}\right)}{\partial x_n}\\ ⋮& ⋮& & ⋮\\ \frac{\partial f_n\left(\mathit{x}\right)}{\partial x_{\mathrm{1}}}& \frac{\partial f_n\left(\mathit{x}\right)}{\partial x_{\mathrm{2}}}& \cdots & \frac{\partial f_n\left(\mathit{x}\right)}{\partial x_n}\end{array}\right]

式(11)的Newton-Raphson迭代格式为

{\mathit{x}}_k={\mathit{x}}_{k-\mathrm{1}}-{\left[\mathit{F}\mathrm{\text{'}}\left({\mathit{x}}_{k-\mathrm{1}}\right)\right]}^{-\mathrm{1}}\mathit{F}\left({\mathit{x}}_{k-\mathrm{1}}\right),k=\mathrm{1},\mathrm{2},\cdots(12)

在实际应用Newton-Raphson迭代法求解非线性微分方程过程中，在计算{\mathit{x}}_k时，一般将式(12)转变为

\left[\mathit{F}\mathrm{\text{'}}\left({\mathit{x}}_{k-\mathrm{1}}\right)\right]{\mathit{x}}_k=\left[\mathit{F}\mathrm{\text{'}}\left({\mathit{x}}_{k-\mathrm{1}}\right)\right]{\mathit{x}}_{k-\mathrm{1}}-\mathit{F}\left({\mathit{x}}_{k-\mathrm{1}}\right),k=\mathrm{1},\mathrm{2},\cdots(13)

这样每次迭代只需计算式(13)即可。

3 矩阵方程及迭代格式

本节主要以(1+1)维非线性分数阶Sobolev方程为例，推导出方程的矩阵方程，并得到方程的Newton-Raphson迭代格式。

此时考虑双变量函数u\left(x,t\right)，定义域为\left[a,b\right]\times \left[\mathrm{0},T\right]，设定义域上有\left(m+\mathrm{1}\right)\times \left(n+\mathrm{1}\right)个张量型插值节点\left(x_i,t_k\right)，i=\mathrm{0},\mathrm{1},\cdots ,m，k=\mathrm{0},\mathrm{1},\cdots ,n。u\left(x,t\right)的重心插值函数为

u\left(x,t\right)={\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}{L}_i\left(x\right)T_k\left(t\right)u_{ik}(14)

式中：u_{ik}=u\left(x_i,t_k\right)，L_i(x)，T_k(t)为空间和时间变量方向上的插值基函数，定义为

L_i\left(x\right)=\frac{\frac{{\omega }_i}{x-x_i}}{{\displaystyle \sum _{h=\mathrm{0}}^m}\frac{{\omega }_h}{x-x_h}},T_k\left(t\right)=\frac{\frac{{\varphi }_k}{t-t_k}}{{\displaystyle \sum _{s=\mathrm{0}}^n}\frac{{\varphi }_s}{t-t_s}}(15)

其中：{\omega }_i，{\varphi }_k为权函数，定义为式(5)或式(6)。

式(14)在插值节点处的\left(k_x+k_t\right)阶导数可表示为

u^{\left(k_x+k_t\right)}(x,t)={\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}{L}_i^{\left(k_x\right)}(x){T_k^{\mathrm{\text{'}}}}^{\left(k_t\right)}(t)u_{ik}(16)

用矩阵的形式表示为

{\mathit{u}}^{\left(k_x+k_t\right)}=\left[{\mathit{D}}^{\left(k_x\mathrm{0}\right)}\otimes {\mathit{D}}^{\left(\mathrm{0}{k}_t\right)}\right]\mathit{u}(17)

其中：{\mathit{u}}^{\left(k_x+k_t\right)}={\left[u_{\mathrm{00}}^{\left(k_x+k_t\right)},\cdots ,u_{\mathrm{0}n}^{\left(k_x+k_t\right)},u_{\mathrm{10}}^{\left(k_x+k_t\right)},\cdots ,u_{mn}^{\left(k_x+k_t\right)}\right]}^{\mathrm{T}}，{\mathit{D}}^{\left(k_x\mathrm{0}\right)}的元素表示为D_{\vartheta i}^{\left(k_x\mathrm{0}\right)}=L_i^{\left(k_x\right)}\left(x_{\vartheta }\right)，{\mathit{D}}^{\left(\mathrm{0}{k}_t\right)}的元素为D_{\gamma k}^{\left(\mathrm{0}{k}_t\right)}=T_k^{\left(k_t\right)}\left(t_{\gamma }\right)，\mathit{u}={\left[u_{\mathrm{00}},\cdots ,u_{\mathrm{0}n},u_{\mathrm{10}},\cdots ,u_{mn}\right]}^{\mathrm{T}}。

对于(1+1)维非线性分数阶Sobolev方程式(1)，由于方程中分数阶导数项定义为\alpha阶的Caputo分数阶导数，对该项进行分部积分来克服奇异性，即

\begin{array}{l}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=\frac{\mathrm{1}}{\mathrm{\Gamma }\left(\mathrm{1}-\alpha \right)}{\int }_{\mathrm{0}}^t\frac{\partial u\left(x,\tau \right)}{\partial \tau }{\left(t-\tau \right)}^{-\alpha }\mathrm{d}\tau \\ =\frac{\mathrm{1}}{\mathrm{\Gamma }\left(\mathrm{1}-\alpha \right)\left(\mathrm{1}-\alpha \right)}\frac{\partial u\left(x,\mathrm{0}\right)}{\partial t}{t}^{\mathrm{1}-\alpha }\\ +\frac{\mathrm{1}}{\mathrm{\Gamma }\left(\mathrm{1}-\alpha \right)\left(\mathrm{1}-\alpha \right)}{\int }_{\mathrm{0}}^t\frac{{\partial }^{\mathrm{2}}u\left(x,\tau \right)}{\partial {\tau }^{\mathrm{2}}}\frac{\mathrm{d}\tau }{{\left(t-\tau \right)}^{\alpha -\mathrm{1}}}\\ ={\mathrm{\Gamma }}_{\alpha }^{\mathrm{1}}\left[\frac{\partial u\left(x,\mathrm{0}\right)}{\partial t}{t}^{\mathrm{1}-\alpha }+{\int }_{\mathrm{0}}^t\frac{{\partial }^{\mathrm{2}}u\left(x,\tau \right)}{\partial {\tau }^{\mathrm{2}}}\frac{\mathrm{d}\tau }{{\left(t-\tau \right)}^{\alpha -\mathrm{1}}}\right]\end{array}

其中{\mathrm{\Gamma }}_{\alpha }^{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{\Gamma }\left(\mathrm{1}-\alpha \right)\left(\mathrm{1}-\alpha \right)}，则(1+1)维非线性分数阶Sobolev方程可以表示为

\begin{array}{l}{\mathrm{\Gamma }}_{\alpha }^{\mathrm{1}}\left[\frac{\partial u\left(x,\mathrm{0}\right)}{\partial t}{t}^{\mathrm{1}-\alpha }+{\int }_{\mathrm{0}}^t\frac{{\partial }^{\mathrm{2}}u\left(x,\tau \right)}{\partial {\tau }^{\mathrm{2}}}\frac{\mathrm{d}\tau }{{\left(t-\tau \right)}^{\alpha -\mathrm{1}}}\right]\\ -k_{\mathrm{1}}{u}_{xxt}-k_{\mathrm{2}}{u}_{xx}=G\left(u\right)+f\left(x,t\right)\end{array}(18)

对于式(18)中的函数u\left(x,t\right)，对其使用双变量函数的重心插值式(14)，此时式(18)为

\begin{array}{l}{\mathrm{\Gamma }}_{\alpha }^{\mathrm{1}}{\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}\left[L_i\left(x\right)T_k^{\mathrm{\text{'}}}\left(\mathrm{0}\right)t^{\mathrm{1}-\alpha }\right]u_{ik}\\ +{\mathrm{\Gamma }}_{\alpha }^{\mathrm{1}}{\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}\left[L_i\left(x\right){\int }_{\mathrm{0}}^t\frac{T_k^{″}\left(\tau \right)}{{\left(t-\tau \right)}^{\alpha -\mathrm{1}}}\mathrm{d}\tau \right]u_{ik}\\ -k_{\mathrm{1}}{\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}{L}_i^{″}\left(x\right)T_k^{\mathrm{\text{'}}}\left(t\right)u_{ik}-k_{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}{L}_i^{″}\left(x\right)T_k\left(t\right)u_{ik}\\ =G\left[{\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}{L}_i^{}\left(x\right)T_k\left(t\right)u_{ik}\right]+f\left(x,t\right)\end{array}(19)

令x=x_{\vartheta }，t=t_{\gamma }，\vartheta =\mathrm{0},\mathrm{1},\cdots ,m，\gamma =\mathrm{0},\mathrm{1},\cdots ,n，则式(19)为

\begin{array}{l}{\mathrm{\Gamma }}_{\alpha }^{\mathrm{1}}{\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}\left[L_i\left(x_{\vartheta }\right)T_k^{\mathrm{\text{'}}}\left(\mathrm{0}\right)t_{\gamma }^{\mathrm{1}-\alpha }\right]u_{ik}\\ +{\mathrm{\Gamma }}_{\alpha }^{\mathrm{1}}{\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}\left[L_i\left(x_{\vartheta }\right){\int }_{\mathrm{0}}^{t_{\gamma }}\frac{T_k^{″}\left(\tau \right)}{{\left(t_{\gamma }-\tau \right)}^{\alpha -\mathrm{1}}}\mathrm{d}\tau \right]u_{ik}\\ -k_{\mathrm{1}}{\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}{L}_i^{″}\left(x_{\vartheta }\right)T_k^{\mathrm{\text{'}}}\left(t_{\gamma }\right)u_{ik}\\ -k_{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}{L}_i^{″}\left(x_{\vartheta }\right)T_k\left(t_{\gamma }\right)u_{ik}\\ =G\left[{\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}{L}_i^{}\left(x_{\vartheta }\right)T_k\left(t_{\gamma }\right)u_{ik}\right]+f\left(x_{\vartheta },t_{\gamma }\right)\end{array}(20)

对式(20)中积分项用权重为\rho \left(\tau \right)={\left(t_{\gamma }-\tau \right)}^{\mathrm{1}-\alpha }的高斯正交公式表示为

\begin{array}{l}{Q}_{\gamma k}=Q_k\left(t_{\gamma }\right)={\int }_{\mathrm{0}}^{t_{\gamma }}\frac{T_k^{″}\left(\tau \right)}{{\left(t_{\gamma }-\tau \right)}^{\alpha -\mathrm{1}}}\mathrm{d}\tau \\ ={\displaystyle \sum _{i=\mathrm{1}}^h}{T}_k^{″}\left({\tau }_i^{\gamma ,\alpha }\right)A_i^{\gamma ,\alpha }\end{array}(21)

式中：A_i^{\gamma ,\alpha }和{\tau }_i^{\gamma ,\alpha }是高斯权和高斯点；h为高斯点个数。

非线性方程组式(20)写成矩阵形式

\begin{array}{l}{\mathrm{\Gamma }}_{\alpha }^{\mathrm{1}}\left[{\mathit{T}}^{\mathrm{1},\alpha }\left({\mathit{I}}_{m+\mathrm{1}}\otimes {\mathit{D}}_{\mathrm{1}}^{\left(\mathrm{01}\right)}\right)+{\mathit{{\rm I}}}_{m+\mathrm{1}}\otimes \mathit{Q}\right]\mathit{U}\\ -k_{\mathrm{1}}{\mathit{D}}^{\left(\mathrm{20}\right)}\otimes {\mathit{D}}_{\mathrm{2}}^{\left(\mathrm{01}\right)}\mathit{U}\\ -k_{\mathrm{2}}{\mathit{D}}^{\left(\mathrm{20}\right)}\otimes {\mathit{I}}_{n+\mathrm{1}}\mathit{U}=G\left(\mathit{U}\right)+\mathit{F}\end{array}(22)

式中：{\mathit{I}}_{m+\mathrm{1}}，{\mathit{I}}_{n+\mathrm{1}}为单位矩阵。则式(22)简写成

\mathit{H}\mathit{U}=\mathit{F}+G\left(\mathit{U}\right)(23)

其中：{\mathit{D}}_{\mathrm{2}}^{\left(\mathrm{01}\right)}，{\mathit{D}}^{\left(\mathrm{20}\right)}分别由式(7)、式(8)所得，以及

\begin{array}{l}\mathit{H}={\mathrm{\Gamma }}_{\alpha }^{\mathrm{1}}\left[{\mathit{T}}^{\mathrm{1},\alpha }\left({\mathit{I}}_{m+\mathrm{1}}\otimes {\mathit{D}}_{\mathrm{1}}^{\left(\mathrm{01}\right)}\right)+{\mathit{I}}_{m+\mathrm{1}}\otimes \mathit{Q}\right]\\ -k_{\mathrm{1}}{\mathit{D}}^{\left(\mathrm{20}\right)}\otimes {\mathit{D}}_{\mathrm{2}}^{\left(\mathrm{01}\right)}-k_{\mathrm{2}}{\mathit{D}}^{\left(\mathrm{20}\right)}\otimes {\mathit{I}}_{n+\mathrm{1}}\\ \mathit{U}={\left[\begin{array}{cccc}\begin{array}{ccc}{u}_{\mathrm{00}},& \cdots ,& u_{\mathrm{0}n},\end{array}& u_{\mathrm{10}},& \cdots ,& u_{mn}\end{array}\right]}^{\mathrm{T}}\\ \mathit{F}={\left[\begin{array}{cccc}{f}_{\mathrm{00}},& \cdots ,& f_{\mathrm{0}n},& \begin{array}{ccc}{f}_{\mathrm{10}},& \cdots ,& f_{mn}\end{array}\end{array}\right]}^{\mathrm{T}}\end{array}

\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(t_{\gamma }^{\mathrm{1}-\alpha }\right)=\left[\begin{array}{cccc}{t}_{\mathrm{0}}^{\mathrm{1}-\alpha }& & & \\ & t_{\mathrm{1}}^{\mathrm{1}-\alpha }& & \\ & & \ddots & \\ & & & t_n^{\mathrm{1}-\alpha }\end{array}\right]

{\mathit{D}}_{\mathrm{1}}^{\left(\mathrm{01}\right)}={\left[\begin{array}{cccc}-{\displaystyle \sum _{k=\mathrm{1}}^n}{D}_{\mathrm{0}k}^{\left(\mathrm{01}\right)}& D_{\mathrm{01}}^{\left(\mathrm{01}\right)}& \cdots & D_{\mathrm{0}n}^{\left(\mathrm{01}\right)}\\ -{\displaystyle \sum _{k=\mathrm{1}}^n}{D}_{\mathrm{0}k}^{\left(\mathrm{01}\right)}& D_{\mathrm{01}}^{\left(\mathrm{01}\right)}& \cdots & D_{\mathrm{0}n}^{\left(\mathrm{01}\right)}\\ ⋮& ⋮& & ⋮\\ -{\displaystyle \sum _{k=\mathrm{1}}^n}{D}_{\mathrm{0}k}^{\left(\mathrm{01}\right)}& D_{\mathrm{01}}^{\left(\mathrm{01}\right)}& \cdots & D_{\mathrm{0}n}^{\left(\mathrm{01}\right)}\end{array}\right]}_{\left(n+\mathrm{1}\right)\times \left(n+\mathrm{1}\right)}

{\mathit{T}}^{\mathrm{1},\alpha }={\left[\begin{array}{cccc}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(t_{\gamma }^{\mathrm{1}-\alpha }\right)& & & \\ & \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(t_{\gamma }^{\mathrm{1}-\alpha }\right)& & \\ & & \ddots & \\ & & & \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(t_{\gamma }^{\mathrm{1}-\alpha }\right)\end{array}\right]}_{z\times z}

z=\left(m+\mathrm{1}\right)\times \left(n+\mathrm{1}\right),G\left(\mathit{U}\right)=\left[\begin{array}{c}G\left(u_{\mathrm{00}}\right)\\ ⋮\\ \begin{array}{l}G\left(u_{\mathrm{0}n}\right)\\ G\left(u_{\mathrm{10}}\right)\end{array}\\ ⋮\\ G\left(u_{mn}\right)\end{array}\right]

方程的矩阵方程式(23)可以表示为R\left(\mathit{U}\right)=\mathit{H}\mathit{U}-G\left(\mathit{U}\right)-\mathit{F}=\mathrm{0}，则其Jacobi导数矩阵为

R\mathrm{\text{'}}\left(\mathit{U}\right)=\mathit{H}-G\mathrm{\text{'}}\left(\mathit{U}\right)

式中：G\mathrm{\text{'}}\left(\mathit{U}\right)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[G\mathrm{\text{'}}\left(u_{\mathrm{00}}\right)\cdots G\mathrm{\text{'}}\left(u_{mn}\right)\right]，G\mathrm{\text{'}}\left(u\right)为函数G\left(u\right)关于u的一阶导函数。

则方程的Newton-Raphson迭代格式为{\mathit{U}}^{\left(\zeta \right)}={\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}-{\left[R\mathrm{\text{'}}\left({\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}\right)\right]}^{-\mathrm{1}}R\left[{\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}\right]，整理得

R\mathrm{\text{'}}\left[{\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}\right]{\mathit{U}}^{\left(\zeta \right)}=R\mathrm{\text{'}}\left[{\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}\right]{\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}-R\left[{\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}\right](24)

式中：{\mathit{U}}^{\left(\mathrm{0}\right)}为给定的迭代初值\widehat{u}构成的列向量；\zeta为迭代次数，\zeta =\mathrm{1},\mathrm{2},\cdots。

对于(2+1)维的非线性分数阶Sobolev方程式(2)，将方程的分数阶导数转化为黎曼积分，再对其使用三变量函数的重心插值配点法进行离散化处理，具体步骤同上述(1+1)维非线性分数阶Sobolev方程离散化过程类似。将离散后的非线性方程组写为矩阵的形式，有

\begin{array}{l}{\mathrm{\Gamma }}_{\alpha }^{\mathrm{1}}\left\{{\mathit{T}}^{\mathrm{1},\alpha }\left[{\mathit{I}}_{m+\mathrm{1}}\otimes {\mathit{I}}_{l+\mathrm{1}}\otimes {\mathit{D}}_{\mathrm{1}}^{\left(\mathrm{001}\right)}\right]+{\mathit{I}}_{m+\mathrm{1}}\otimes {\mathit{I}}_{l+\mathrm{1}}\otimes \mathit{Q}\right\}\mathit{U}\\ -k_{\mathrm{1}}\left[{\mathit{D}}^{\left(\mathrm{200}\right)}\otimes {\mathit{I}}_{l+\mathrm{1}}\otimes {\mathit{D}}_{\mathrm{2}}^{\left(\mathrm{001}\right)}+{\mathit{I}}_{m+\mathrm{1}}\otimes {\mathit{D}}^{\left(\mathrm{020}\right)}\otimes {\mathit{D}}_{\mathrm{2}}^{\left(\mathrm{001}\right)}\right]\mathit{U}\\ -k_{\mathrm{2}}\left[{\mathit{D}}^{\left(\mathrm{200}\right)}\otimes {\mathit{I}}_{l+\mathrm{1}}\otimes {\mathit{I}}_{n+\mathrm{1}}+{\mathit{I}}_{m+\mathrm{1}}\otimes {\mathit{D}}^{\left(\mathrm{020}\right)}\otimes {\mathit{I}}_{n+\mathrm{1}}\right]\mathit{U}\\ =\mathit{F}+G\left(\mathit{U}\right)\end{array}(25)

简写为\mathit{K}\mathit{U}=\mathit{F}+G\left(\mathit{U}\right)，且

\begin{array}{l}\mathit{K}={\mathrm{\Gamma }}_{\alpha }^{\mathrm{1}}\left\{{\mathit{T}}^{\mathrm{1},\alpha }\left[{\mathit{I}}_{m+\mathrm{1}}\otimes {\mathit{I}}_{l+\mathrm{1}}\otimes {\mathit{D}}_{\mathrm{1}}^{\left(\mathrm{001}\right)}\right]+{\mathit{I}}_{m+\mathrm{1}}\otimes {\mathit{I}}_{l+\mathrm{1}}\otimes \mathit{Q}\right\}\\ -k_{\mathrm{1}}\left[{\mathit{D}}^{\left(\mathrm{200}\right)}\otimes {\mathit{I}}_{l+\mathrm{1}}\otimes {\mathit{D}}_{\mathrm{2}}^{\left(\mathrm{001}\right)}+{\mathit{I}}_{m+\mathrm{1}}\otimes {\mathit{D}}^{\left(\mathrm{020}\right)}\otimes {\mathit{D}}_{\mathrm{2}}^{\left(\mathrm{001}\right)}\right]\\ -k_{\mathrm{2}}\left[{\mathit{D}}^{\left(\mathrm{200}\right)}\otimes {\mathit{I}}_{l+\mathrm{1}}\otimes {\mathit{I}}_{n+\mathrm{1}}+{\mathit{I}}_{m+\mathrm{1}}\otimes {\mathit{D}}^{\left(\mathrm{020}\right)}\otimes {\mathit{I}}_{n+\mathrm{1}}\right]\end{array}(26)

此时方程的矩阵方程可以表示为P\left(\mathit{U}\right)=\mathit{K}\mathit{U}-G\left(\mathit{U}\right)-\mathit{F}=\mathrm{0}，则其Jacobi导数矩阵为

P\mathrm{\text{'}}\left(\mathit{U}\right)=\mathit{K}-G\mathrm{\text{'}}\left(\mathit{U}\right)

式中：G\mathrm{\text{'}}\left(\mathit{U}\right)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[G\mathrm{\text{'}}\left(u_{\mathrm{000}}\right)\cdots G\mathrm{\text{'}}\left(u_{mln}\right)\right]，G\mathrm{\text{'}}\left(u\right)为函数G\left(u\right)关于u的一阶导函数。

则方程的Newton-Raphson迭代格式为

{\mathit{U}}^{\left(\zeta \right)}={\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}-{\left[P\mathrm{\text{'}}\left({\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}\right)\right]}^{-\mathrm{1}}P\left[{\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}\right]

整理得

P\mathrm{\text{'}}\left[{\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}\right]{\mathit{U}}^{\left(\zeta \right)}=P\mathrm{\text{'}}\left[{\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}\right]{\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}-P\left[{\mathit{U}}^{\left(\zeta -\mathrm{1}\right)}\right](27)

关于(1+1)维和(2+1)维非线性分数阶Sobolev方程初始和边界条件的离散和施加方法见文献[21]，此处不再赘述。

施加完初始和边界条件后，得到最终的矩阵方程，令\zeta =\mathrm{1}，给定一个迭代精度e，当满足

\begin{array}{l}{‖u^{\left(\zeta \right)}\left(x,t\right)-u^{\left(\zeta -\mathrm{1}\right)}\left(x,t\right)‖}_{\mathrm{2}}\\ =\sqrt[]{{\displaystyle \sum _{i=\mathrm{0}}^m}{\displaystyle \sum _{k=\mathrm{0}}^n}{\left[u^{\left(\zeta \right)}\left(x_i,t_k\right)-u^{\left(\zeta -\mathrm{1}\right)}\left(x_i,t_k\right)\right]}^{\mathrm{2}}}\le e\end{array}(28)

时，迭代终止，即可得到所求方程的数值解。式(28)中{‖·‖}_{\mathrm{2}}表示欧几里得范数。本文预先规定迭代精度为10^-6。

4 数值算例

为验证重心插值Newton-Raphson迭代法对求解非线性分数阶Sobolev方程的有效性，本节给出了三个非线性项不同且维数不同的非线性分数阶Sobolev方程的数值算例，所有算例均使用MATLAB(R2020b)版本实现。

绝对误差、相对误差的定义如下：

E_{\mathrm{a}}={‖\tilde{\mathit{u}}-\mathit{u}‖}_{\mathrm{2}},E_{\mathrm{r}}=\frac{{‖\tilde{\mathit{u}}-\mathit{u}‖}_{\mathrm{2}}}{{‖\mathit{u}‖}_{\mathrm{2}}}

其中：\tilde{\mathit{u}}，\mathit{u}分别表示由数值解和精确解构成的列向量。

算例1 考虑(1+1)维非线性分数阶Sobolev方程

\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}-\mathrm{\Delta }{u}_t\left(x,t\right)-\mathrm{\Delta }u\left(x,t\right)+u^{\mathrm{2}}=f\left(x,t\right)

其中：x,t\in \left[\mathrm{0},\mathrm{1}\right]，精确解为u\left(x,t\right)=t^{\mathrm{4}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }x\right)，方程对应的初始和边界条件可由精确解得到，且

\begin{array}{l}f\left(x,t\right)=\frac{\mathrm{24}{t}^{\left(\mathrm{4}-\alpha \right)}}{\mathrm{\Gamma }\left(\mathrm{5}-\alpha \right)}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }x\right)+{\mathrm{\pi }}^{\mathrm{2}}{t}^{\mathrm{4}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }x\right)\\ +\mathrm{4}{\mathrm{\pi }}^{\mathrm{2}}{t}^{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }x\right)+{\left[t^{\mathrm{4}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }x\right)\right]}^{\mathrm{2}}\end{array}

表1和表2给出了\alpha =\mathrm{0.5}时，重心Lagrange插值Newton-Raphson迭代法在不同节点类型的误差和收敛次数。通过结果可以看出，当选用切比雪夫节点时，重心Lagrange插值Newton-Raphson迭代法随着节点数量的增加，误差先减小后增大，且迭代次数一直维持在4次。当选用等距节点时，随着节点数量的增加，误差先减小后急剧增大，且表现出极大的数值不稳定性，当m=n=\mathrm{18}时，迭代次数增大到15次，当节点数量为24时，迭代次数达到1 000也无法满足式(28)，迭代终止。在切比雪夫节点上得到的误差小于在等距节点上得到的误差。

表3和表4给出了\alpha =\mathrm{0.5}，d_x=d_t=\mathrm{6}时，重心有理插值Newton-Raphson迭代法在不同节点类型的误差和收敛次数。通过结果可以看出，当选用切比雪夫节点时，重心有理插值Newton-Raphson迭代法随着节点数量的增加，误差减小，迭代次数一直维持在4次。当选用等距节点时，随着节点数量的增加，误差减小，迭代次数一直维持在4次，且在切比雪夫节点得到的误差小于在等距节点得到的误差。

表5给出了m=n=\mathrm{12}时，不同\alpha值下，重心Lagrange插值Newton-Raphson迭代法在不同节点类型的误差。通过结果可以看出，选取切比雪夫节点所得到的误差小于选取等距节点所得到的误差。

表6给出了m=n=\mathrm{12}，d_x=d_t=\mathrm{6}时，不同\alpha值下，重心有理插值Newton-Raphson迭代法在不同节点类型上的误差。通过结果可以看出，选取切比雪夫节点所得到的误差小于选取等距节点所得到的误差。

图1给出了重心Lagrange插值Newton-Raphson迭代法在m=n=\mathrm{12}，\alpha =\mathrm{0.5}，选取切比雪夫节点时的数值解和解析解。通过比较，可以得知解析解和数值解的吻合度较高。

图2给出了重心有理插值Newton-Raphson迭代法在m=n=\mathrm{12}，\alpha =\mathrm{0.5}，d_x=d_t=\mathrm{6},选取不同节点类型时的误差分布。可以看出，选取切比雪夫节点时得到的误差小于选取等距节点时得到的误差。

算例2 考虑(2+1)维非线性分数阶Sobolev方程

\begin{array}{l}\frac{{\partial }^{\alpha }u\left(x,y,t\right)}{\partial t^{\alpha }}-\mathrm{\Delta }{u}_t\left(x,y,t\right)-\mathrm{\Delta }u\left(x,y,t\right)\\ =f\left(x,y,t\right)+u(\mathrm{1}-u),\end{array}

其中：x,y,t\in \left[\mathrm{0},\mathrm{1}\right]，精确解为u\left(x,y,t\right)=\mathrm{s}\mathrm{i}\mathrm{n}\left(t\right){\mathrm{e}}^{x-y}，方程对应的初始和边界条件可由精确解得到，且

\begin{array}{l}f\left(x,y,t\right)\\ =\left[t^{\mathrm{1}-\alpha }{M}_{\mathrm{2},\mathrm{2}-\alpha }\left(-t^{\mathrm{2}}\right)-\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}\left(t\right)-\mathrm{3}\mathrm{s}\mathrm{i}\mathrm{n}\left(t\right)\right]{\mathrm{e}}^{x-y}\\ +\left[{\left(\mathrm{s}\mathrm{i}\mathrm{n}\left(t\right)\right)}^{\mathrm{2}}{\mathrm{e}}^{x-y}\right]{\mathrm{e}}^{x-y}\end{array}

f\left(x,y,t\right)表达式中的M_{a,b}表示包含的两个参数a,b\in {\mathbb{R}}^{+}的Mittag-Leffler函数，定义为

M_{a,b}\left(s\right)={\displaystyle \sum _{h=\mathrm{0}}^{\mathrm{\infty }}}\frac{s^h}{\mathrm{\Gamma }\left(ah+b\right)},\forall s\in \mathbb{C}

且规定M_{\mathrm{1},\mathrm{1}}\left(s\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left(s\right),\forall s\in \mathbb{C}。

表7和表8给出了\alpha =\mathrm{0.5}时，重心Lagrange插值Newton-Raphson迭代法在不同节点类型的误差和收敛次数。通过结果可以看出，当选用切比雪夫节点时，重心Lagrange插值Newton-Raphson迭代法随着节点数量的增加，误差先减小后增大，当m=l=n=\mathrm{6}时，迭代次数为3，之后迭代次数变为4；当选用等距节点时，随着节点数量的增加，误差先减小后急剧增大，当m=l=n=\mathrm{12}时，迭代次数增大到11，当节点数量为14时，迭代次数达到1 000也无法满足式(28)，迭代终止。在切比雪夫节点处得到的误差小于在等距节点处得到的误差。

表9和表10给出了\alpha =\mathrm{0.5}，d_x=d_y=d_t=\mathrm{6}时，重心有理插值Newton-Raphson迭代法在不同节点类型的误差和收敛次数。通过结果可以看出，当选用切比雪夫节点时，重心有理插值Newton-Raphson迭代法随着节点数量的增加，误差先减小增大，当m=l=n=\mathrm{6}时，迭代次数为3，之后迭代次数变为4；当选用等距节点时，随着节点数量的增加，误差减小，迭代次数均为4。在切比雪夫节点处得到的误差小于在等距节点处得到的误差。

表11给出了m=l=n=\mathrm{12}，不同\alpha值下，重心Lagrange插值Newton-Raphson迭代法在不同节点类型的误差。通过结果可以看出，选取切比雪夫节点所得到的误差小于选取等距节点所得到的误差。

表12给出了m=l=n=\mathrm{12}，d_x=d_y=d_t=\mathrm{6}时，不同\alpha值下，重心有理插值Newton-Raphson迭代法在不同节点类型上的误差。通过结果可以看出，选取切比雪夫节点所得到的误差小于选取等距节点所得到的误差。

图3给出了重心Lagrange插值Newton-Raphson迭代法在m=l=n=\mathrm{12}，\alpha =\mathrm{0.5}，t=\mathrm{1}时，选取切比雪夫节点的数值解和解析解。通过比较，可以得知解析解和数值解的吻合度较高。

图4给出了重心有理插值Newton-Raphson迭代法在m=l=n=\mathrm{10}，\alpha =\mathrm{0.5}，d_x=d_y=d_t=\mathrm{6}，t=\mathrm{1}时，选取不同节点类型时的误差分布。可以看出，选取切比雪夫节点得到的误差小于选取等距节点得到的误差。

算例3 考虑(2+1)维非线性分数阶Sobolev方程

\begin{array}{l}\frac{{\partial }^{\alpha }u\left(x,y,t\right)}{\partial t^{\alpha }}-\mathrm{2}\mathrm{\Delta }{u}_t\left(x,y,t\right)-\mathrm{\Delta }u\left(x,y,t\right)\\ =f\left(x,y,t\right)+u^{\mathrm{3}}\end{array}

其中：x,y,t\in \left[\mathrm{0},\mathrm{1}\right]，精确解为u\left(x,y,t\right)={\mathrm{e}}^{x-y-t}，方程对应的初始和边界条件可由精确解得到，且

f\left(x,y,t\right)=-t^{\mathrm{1}-\alpha }{M}_{\mathrm{1},\mathrm{2}-\alpha }\left(-t\right){\mathrm{e}}^{x-y}+\mathrm{2}{\mathrm{e}}^{x-y-t}-{\left({\mathrm{e}}^{x-y-t}\right)}^{\mathrm{3}}

表13和表14给出了\alpha =\mathrm{0.5}时，重心Lagrange插值Newton-Raphson迭代法在不同节点类型的误差和收敛次数。通过结果可以看出，当选用切比雪夫节点时，重心Lagrange插值Newton-Raphson迭代法随着节点数量的增加，误差先减小后增大，迭代次数始终保持在4；当选用等距节点时，随着节点数量的增加，误差先减小后急剧增大，当m=l=n=\mathrm{12}时，迭代次数增大到5，当节点数量为14时，迭代次数达到1 000也无法满足式(28)，迭代终止。在切比雪夫节点得到的误差小于在等距节点得到的误差。

表15和表16给出了\alpha =\mathrm{0.5}，d_x=d_y=d_t=\mathrm{6}时，重心有理插值Newton-Raphson迭代法在不同节点类型的误差和收敛次数。通过结果可以看出，当选用切比雪夫节点，重心有理插值Newton-Raphson迭代法随着节点数量的增加，误差先减小后增大，迭代次数始终保持在4；当选用等距节点时，随着节点数量的增加，误差减小，迭代次数始终保持在4。在切比雪夫节点得到的误差小于在等距节点得到的误差。

图5给出了重心Lagrange插值Newton-Raphson迭代法在m=l=n=\mathrm{12}，\alpha =\mathrm{0.5}，t=\mathrm{1}时，选取切比雪夫节点的数值解和解析解。通过比较，可以得知解析解和数值解的吻合度较高。

5 结论

本文使用重心插值Newton-Raphson迭代法求解了(1+1)维和(2+1)维非线性分数阶Sobolev方程。对Caputo分数阶导数部分进行分部积分，使其转化为黎曼积分，再利用重心插值配点法对方程进行了离散，得到对应迭代格式下的矩阵方程，构造了对应的Newton-Raphson迭代格式。数值算例的结果表明，选取切比雪夫节点所得到的误差小于选取等距节点所得到的误差，且数值解和解析解的吻合度较高，重心插值Newton-Raphson迭代法对于求解非线性分数阶Sobolev方程是有效的。

参考文献

[1] AL-SMADI M, ARQUB O A, HADID S. An attractive analytical technique for coupled system of fractional partial differential equations in shallow water waves with conformable derivative[J]. Communications in Theoretical Physics, 2020, 72(8): 085001.

[2] ZAYED E M E, AMER Y A, SHOHIB R M A. The fractional complex transformation for nonlinear fractional partial differential equations in the mathematical physics[J]. Journal of the Association of Arab Universities for Basic and Applied Sciences, 2016, 19: 59-69.

[3] SINGH J, KUMAR D, KİLİÇMAN A. Numerical solutions of nonlinear fractional partial differential equations arising in spatial diffusion of biological populations[J]. Abstract and Applied Analysis, 2014, 2014: 535793.

[4] CHIYANEH A B, DURU H. On adaptive mesh for the initial boundary value singularly perturbed delay Sobolev problems[J]. Numerical Methods for Partial Differential Equations, 2020, 36(2): 228-248.

[5] KUMBINARASAIAH S. Numerical solution for the (2+1)dimensional Sobolev and regularized long wave equations arise in fluid mechanics via wavelet technique[J]. Partial Differential Equations in Applied Mathematics, 2021, 3: 100016.

[6] HEYDARI M H, ZHAGHARIAN S, RAZZAGHI M. Discrete Chebyshev polynomials for the numerical solution of stochastic fractional two-dimensional Sobolev equation[J]. Communications in Nonlinear Science and Numerical Simulation, 2024, 130: 107742.

[7] GUO T, NIKAN O, QIU W L, et al. Localized meshless approaches based on theta method and BDF2 for nonlinear Sobolev equation arising from fluid dynamics[J]. Communications in Nonlinear Science and Numerical Simulation, 2023, 117: 106989.

[8] AZIN H, HABIBIRAD A, BAGHANI O. Legendre-finite difference method for solving fractional nonlinear Sobolev equation with Caputo derivative[J]. Journal of Computational Science, 2023, 74: 102177.

[9] NIU Y X, LIU Y, LI H, et al. Fast high-order compact difference scheme for the nonlinear distributed-order fractional Sobolev model appearing in porous media[J]. Mathematics and Computers in Simulation, 2023, 203: 387-407.

[10] QIN Y F, YANG X C, REN Y Z, et al. A Newton linearized Crank-Nicolson method for the nonlinear space fractional sobolev equation[J]. Journal of Function Spaces, 2021, 2021(1): 9979791.

[11] HAQ S, HUSSAIN M. Application of meshfree spectral method for the solution of multi-dimensional time-fractional Sobolev equations[J]. Engineering Analysis With Boundary Elements, 2019, 106: 201-216.

[12] LI J, CHENG Y L. Barycentric rational interpolation method for solving KPP equation[J]. Electronic Research Archive, 2023, 31(5): 3014-3029.

[13] LI J, CHENG Y L. Barycentric rational method for solving biharmonic equation by depression of order[J]. Numerical Methods for Partial Differential Equations, 2021, 37(3): 1993-2007.

[14] LI J, CHENG Y L. Linear barycentric rational collocation method for solving second-order volterra integro-differential equation[J]. Computational and Applied Mathematics, 2020, 39(2): 92.

[15] TORKAMAN S, LOGHMANI G B, HEYDARI M, et al. Novel numerical solutions of nonlinear heat transfer problems using the linear barycentric rational interpolation[J]. Heat Transfer-Asian Research, 2019, 48: 1318-1344.

[16] LI J, SU X N, QU J Z. Linear barycentric rational collocation method for solving telegraph equation[J]. Mathematical Methods in the Applied Sciences, 2021, 44(14): 11720-11737.

[17] LI J, SANG Y. Linear barycentric rational collocation method for beam force vibration equation[J]. Shock and Vibration, 2021, 2021(1): 5584274.

[18] LI J, QU J Z. Barycentric Lagrange interpolation collocation method for solving the sine-Gordon equation[J]. Wave Motion, 2023, 120: 103159.

[19] 王兆清, 李淑萍. 非线性问题的重心插值配点法[M]. 北京: 国防工业出版社, 2015.

[20] WANG Z Q, JIANG J, TANG B T, et al. Barycentric interpolation Newton-Raphson iterative method for solving nonlinear beam equations[J]. Applied Mechanics and Materials, 2014, 684: 41-48.

[21] 屈金铮, 李金, 苏晓宁. 求解非线性伪抛物方程的重心Lagrange插值配点法[J]. 山东大学学报(理学版), 2023, 58(4): 29-39.

基金资助

河北省自然科学基金项目(A2019209533)

山东省自然科学基金项目(ZR2022MA003)