关于非局部扩散模型的一种快速预处理算法

摘要

变系数非局部扩散模型可以被一种快速配置法进行有效的数值离散。离散后得到一个系数矩阵具有 Toeplitz 结构且稠密的线性方程组。由于系数矩阵是非对称的，该线性方程组可以用广义极小残量法（GMRES）方法求解。为了提高 GMRES 方法的收敛率，构造了系数矩阵的 Toeplitz 及循环预处理子，并提出了预处理 GMRES 方法求解该线性方程组。数值算例也表明了该预处理算法的有效性。

关键词

非局部扩散模型; 快速配置法; Toeplitz矩阵; 预处理; GMRES方法

对于自然界那些用整数阶扩散方程不太好描述的现象，如反常扩散和多孔介质流等${}^{[\mathrm{1}-\mathrm{3}]}$，除了可以用分数阶扩散方程建模以外，非局部扩散模型也是一种可供选择的建模方法。非局部扩散理论对间断性问题及其他连续体变形问题给出了一个恰当的描述，而裂缝、断裂等奇异性问题在经典理论方程下无法给出恰当描述。非局部理论被引入到连续介质力学当中，2000年Silling提出了一种叫做peridynamic的非局部扩散模型${}^{[\mathrm{4}-\mathrm{7}]}$，在非局部扩散模型框架下，物体内部不再有接触力，并且可以保证建立的方程都保持同一积分形式，对间断或其他奇异性连续体的变形问题提供了新的研究思路。

然而，利用一般的数值方法数值离散非局部扩散模型，通常得到一个稠密的刚度矩阵，利用直接法求解的计算量和存储量非常大，从而阻碍了其广泛应用。最近Wang C等提出了一种快速的配置法数值离散变系数非局部扩散模型${}^{[\mathrm{8}]}$，分析了刚度矩阵的结构，并利用直接法和共轭梯度平方法求解了变系数非局部扩散模型经过快速配置方法数值离散得到的线性方程组。因为直接法的计算量太大，共轭梯度平方法不稳定，所以这两个方法不实用。

由于变系数非局部扩散模型经过快速配置方法数值离散得到的线性方程组系数矩阵为非对称的，因此可以用间接法中的GMRES方法${}^{[\mathrm{9}]}$求解。在实际求解过程中，系数矩阵的坏条件数会导致GMRES方法收敛的很慢甚至不收敛，而提高收敛速度最常用的方法是预处理技术${}^{[\mathrm{10}-\mathrm{13}]}$。一个好的预处理子应该至少保持某种特殊结构，并使得预处理后的系数矩阵的条件数很小。在文献［8］的研究基础上，本文采用了Toeplitz及循环预处理GMRES方法求解了变系数非局部扩散模型经过快速配置方法数值离散得到的线性方程组。

1 问题的来源

本节给出了所需要求解的线性方程组，即变系数非局部扩散模型经过快速配置法数值离散后所得到的离散方程。

1.1　非局部扩散模型及其快速配置法

变系数非局部扩散模型的一类体积约束非局部Dirichlet边值问题，可以表示如下${}^{[\mathrm{5}]}$：

$$\left\{\begin{array}{ll}\begin{array}{l}{\int }_{x-\delta }^{x+\delta }(\alpha ( x)+\alpha (y))\\ \hspace{1em}\sigma (x-y)(u(x)-u(y))\\ \hspace{1em}\mathrm{d}y=f(x)\end{array}& \\ \hspace{1em}\hspace{1em}\hspace{1em}x\in (a,b)& \\ u(x)=g(x)& \\ x\in (a-\delta ,a]\bigcup [b,b+\delta ]& \end{array}\right.$$

（1）

式中： $\delta >\mathrm{0}$表示非局部扩散现象扩散的水平方向的范围； $\alpha (\cdot )$用来表示可变的扩散系数，并假设其有正的上下界且光滑。核被定义为

$$\sigma (x)=\frac{\mathrm{1}}{{\left|x\right|}^{\mathrm{1}+\gamma }}$$

（2）

指数$\gamma$$<\mathrm{1}$表示核强度，因此积分算子是弱奇异的。取$N$为正整数，网格步长$h=\frac{b-a}{N}$，定义$K=⌊\frac{\delta }{h}⌋$为$\frac{\delta }{h}$的取整。在区间$(a-\delta ,b+\delta )$上，定义一个均匀剖分：$x_i=a+ih,-K\le i\le N+K\mathrm{。}$令

$$\psi (\xi )=\left\{\begin{array}{ll}\mathrm{1}-\left|\xi \right|& \xi \in [-\mathrm{1,1}]\\ \mathrm{0}& \xi \notin [-\mathrm{1,1}]\end{array}\right.$$

又令${\varphi }_i(x)=\psi (\frac{(x-x_i)}{h})$是以$x=x_i,i=\mathrm{1,2},\mathrm{3},\dots ,N$为中心的Hat函数。选取试验函数$u$为

$$u(x)={\displaystyle \sum _{j=-K}^{N+K}}{u}_j{\varphi }_j(x)={\displaystyle \sum _{j=\mathrm{1}}^{N-\mathrm{1}}}{u}_j{\varphi }_j(x)+({\displaystyle \sum _{j=\mathrm{0}}^{-K}}+$$

$$\begin{array}{l}\\ \hspace{1em}\hspace{1em}\hspace{1em}{\displaystyle \sum _{j=N}^{N+K}})g(x_j){\varphi }_j(x)\end{array}$$

（3）

将试验函数式（3）代入式（1），并使控制方程取配置点$\{x_i{\}}_{i=\mathrm{1}}^{N-\mathrm{1}}$处的值，得

$$\begin{array}{l}{\int }_{x_i-\delta }^{x_i+\delta }\frac{\alpha (x_i)+\alpha (y)}{{\left|x_i-y\right|}^{\mathrm{1}+\gamma }}(u(x_i)\hspace{1em}-\\ {\displaystyle \sum _{j=-K}^{N+K}}{u}_j{\varphi }_j(y))\mathrm{d}y=f(x_i)\mathrm{ }\mathrm{ }\\ \hspace{1em}\hspace{1em}\mathrm{1}\le i\le N-\mathrm{1}\end{array}$$

（4）

在节点值$u_j$处，方程（4）可表示如下：

$$\begin{array}{ll}{u}_i{\int }_{x_i-\delta }^{x_i+\delta }\frac{(\alpha (x_i)+\alpha (y))(\mathrm{1}-{\varphi }_i(y))}{{\left|x_i-y\right|}^{\mathrm{1}+\gamma }}\mathrm{d}y-& \\ \hspace{1em}{\displaystyle \sum _{j=\mathrm{1},j\ne i}^{N-\mathrm{1}}}{u}_j{\int }_{x_i-\delta }^{x_i+\delta }\frac{(\alpha (x_i)+\alpha (y)){\varphi }_j(y)}{{\left|x_i-y\right|}^{\mathrm{1}+\gamma }}\mathrm{d}y& =\\ \begin{array}{l}\hspace{1em}\hspace{1em}f(x_i)+({\displaystyle \sum _{j=\mathrm{0}}^{-K}}+{\displaystyle \sum _{j=N}^{N+K}})\mathrm{ }\\ g(x_j){\int }_{x_i-\delta }^{x_i+\delta }\frac{(\alpha (x_i)+\alpha (y)){\varphi }_j(y)}{{\left|x_i-y\right|}^{\mathrm{1}+\gamma }}\mathrm{d}y\end{array}& \\ \hspace{1em}\hspace{1em}\mathrm{1}\le i\le N-\mathrm{1}& \end{array}$$

（5）

数值格式（5）可以表示如下：

$$\mathbf{A}\mathbf{u}=\mathit{f}$$

（6）

这里， $\mathit{u}={\left[u_{\mathrm{1}},u_{\mathrm{2}},\cdots ,u_{\mathit{N}-\mathrm{1}}\right]}^{\mathrm{{\rm T}}}$其中${\left\{u_{\mathit{j}}\right\}}_{j=\mathrm{1}}^{N-\mathrm{1}}$由式（3）给出，刚度矩阵${\mathit{A}}_{}={\left[A_{\mathit{i},\mathit{j}}\right]}_{i,j=\mathrm{1}}^{N-\mathrm{1}}$和右端项${\mathit{f}}_{}={\left[f_{\mathrm{1}},f_{\mathrm{2}},\cdots ,f_{\mathit{N}-\mathrm{1}}\right]}^{\mathrm{{\rm T}}}$被定义如下：

$$\begin{array}{l}{A}_{\mathit{i},\mathit{j}}=-{\int }_{x_i-\delta }^{x_i+\delta }\frac{(\alpha (x_i)+\alpha (y)){\varphi }_j(y)}{{\left|x_i-y\right|}^{\mathrm{1}+\gamma }}\mathrm{d}y\hspace{1em}\hspace{1em}j\ne i\\ A_{\mathit{i},\mathit{i}}={\int }_{x_i-\delta }^{x_i+\delta }\frac{(\alpha (x_i)+\alpha (y))(\mathrm{1}-{\varphi }_i(y))}{{\left|x_i-y\right|}^{\mathrm{1}+\gamma }}\mathrm{d}y\end{array}$$

（7）

由式（7）可以得知，刚度矩阵$\mathit{A}$的每一行都有非零元且，随着增加而增加。换句话说，刚度矩阵$\mathit{A}$是一个渐进稠密的矩阵，直接法求解需要的储存空间和的计算量。

1.2　刚度矩阵的元素估计

参考文献［8］首先估计了式（7）中的刚度矩阵$\mathit{A}$的对角线元素，对于有

$$\begin{array}{l}{A}_{\mathit{i},\mathit{i}}=\alpha (x_i){\int }_{x_i-\delta }^{x_i+\delta }\frac{\mathrm{1}-{\varphi }_i(y)}{{\left|x_i-y\right|}^{\mathrm{1}+\gamma }}\mathrm{d}y+\\ \hspace{1em}\hspace{1em}{\int }_{x_i-\delta }^{x_i+\delta }\frac{\alpha (y)(\mathrm{1}-{\varphi }_i(y))}{{\left|x_i-y\right|}^{\mathrm{1}+\gamma }}\mathrm{d}y\end{array}$$

（8）

通过计算得到了刚度矩阵对角线元素近似值为$$\begin{array}{l}{A}_{\mathit{i},\mathit{i}}\approx \mathrm{4}\alpha (x_i)(\frac{h^{-\gamma }-{\delta }^{-\gamma }}{\gamma }+\frac{h^{-\gamma }}{\mathrm{1}-\gamma })+\\ \hspace{1em}\hspace{1em}{\alpha }^{\text{'}}(x_i)(\frac{{\delta }^{\mathrm{2}-\gamma }-h^{\mathrm{2}-\gamma }}{\mathrm{2}-\gamma }+\frac{h^{\mathrm{2}-\gamma }}{\mathrm{3}-\gamma })\end{array}$$

（9）

然后计算式（7）中的刚度矩阵$\mathit{A}$的非对角线元素，有

$$\begin{array}{l}{A}_{\mathit{i},\mathit{j}}=-\alpha (x_i){\int }_{x_i-\delta }^{x_i+\delta }\frac{{\varphi }_j(y)}{{\left|x_i-y\right|}^{\mathrm{1}+\gamma }}\mathrm{d}y-\\ \hspace{1em}\hspace{1em}{\int }_{x_i-\delta }^{x_i+\delta }\frac{\alpha (y){\varphi }_j(y)}{{\left|x_i-y\right|}^{\mathrm{1}+\gamma }}\mathrm{d}y\hspace{1em}\hspace{1em}\mathrm{ }j\ne i\end{array}$$

（10）

由于的支集为，则有

$$A_{\mathit{i},\mathit{j}}=\mathrm{0}\mathrm{ }\mathrm{ }\left|i-j\right|\ge \mathrm{ }K+\mathrm{2}$$

（11）

因此，刚度矩阵$\mathit{A}$是一个带状矩阵，带状的宽度为，故只需计算$A_{\mathit{i},\mathit{j}}$在的不同情形下的数值。在每一区间上，通过近似，由参考文献［8］得到如下结果：

情形1：.的支集为，则有

$$\begin{array}{l}{A}_{\mathit{i},\mathit{j}}\approx -\alpha (x_{j-\frac{\mathrm{1}}{\mathrm{2}}})l_{j-i}-\alpha (x_{j+\frac{\mathrm{1}}{\mathrm{2}}})r_{j-i}-\alpha (x_i)(l_{j-i}+\\ \hspace{1em}\hspace{1em}{r}_{j-i})\end{array}$$

（12）

式（12）中的和定义如下，其中，即

$$\begin{array}{l}{l}_k=\left\{\begin{array}{ll}\begin{array}{l}\frac{h^{-\gamma }}{\gamma }(\mathrm{1}-k)[{(k-\mathrm{1})}^{-\gamma }-k^{-\gamma }]+\\ \hspace{1em}\frac{h^{-\gamma }}{\mathrm{1}-\gamma }[k^{\mathrm{1}-\gamma }-{(k-\mathrm{1})}^{\mathrm{1}-\gamma }]\end{array}& \\ \hspace{1em}\hspace{1em}\gamma \ne \mathrm{0}& \\ \mathrm{1}+(\mathrm{1}-k)\mathrm{l}\mathrm{n}(\frac{k}{k-\mathrm{1}})& \\ \hspace{1em}\hspace{1em}\gamma =\mathrm{0}& \end{array}\right.\\ r_k=\left\{\begin{array}{ll}\begin{array}{l}\frac{h^{-\gamma }}{\gamma }(k+\mathrm{1})[k^{-\gamma }-{(k+\mathrm{1})}^{-\gamma }]+\\ \hspace{1em}\frac{h^{-\gamma }}{\mathrm{1}-\gamma }[k^{\mathrm{1}-\gamma }-{(k+\mathrm{1})}^{\mathrm{1}-\gamma }]\end{array}& \\ \hspace{1em}\hspace{1em}\gamma \ne \mathrm{0}& \\ (k+\mathrm{1})\mathrm{l}\mathrm{n}(\frac{k+\mathrm{1}}{k})& \\ \hspace{1em}\gamma =\mathrm{0}& \end{array}\right.\end{array}$$

情形2： ，则有

$$\begin{array}{l}{A}_{\mathit{i},\mathit{j}}\approx -\alpha (x_{j+\frac{\mathrm{1}}{\mathrm{2}}})l_{i-j}-\alpha (x_{j-\frac{\mathrm{1}}{\mathrm{2}}})r_{i-j}-\\ \hspace{1em}\hspace{1em}\alpha (x_i)(l_{i-j}+r_{i-j})\end{array}$$

（13）

情形3：，的支集，则有

$$\begin{array}{l}{A}_{\mathit{i},\mathit{i}+\mathit{K}}\approx -\alpha (x_{i+K-\frac{\mathrm{1}}{\mathrm{2}}})l_K-\alpha (x_{i+K+\frac{\mathrm{1}}{\mathrm{2}}})r_K-\\ \hspace{1em}\hspace{1em}\alpha (x_i)(l_K+r_K)\end{array}$$

（14）

其中

$$\begin{array}{l}\begin{array}{l}{l}_k=\\ \left\{\begin{array}{ll}\begin{array}{l}\frac{h^{-\gamma }}{\gamma }(\mathrm{1}-K)[{(K-\mathrm{1})}^{-\gamma }-K^{-\gamma }]+\\ \hspace{1em}\frac{h^{-\gamma }}{\mathrm{1}-\gamma }[K^{\mathrm{1}-\gamma }-{(K-\mathrm{1})}^{\mathrm{1}-\gamma }]\end{array}& \\ \hspace{1em}\hspace{1em}\gamma \ne \mathrm{0}& \\ \mathrm{1}+(\mathrm{1}-K)\mathrm{l}\mathrm{n}(\frac{K}{K-\mathrm{1}})& \\ \hspace{1em}\gamma =\mathrm{0}& \end{array}\right.\end{array}\\ r_k=\left\{\begin{array}{ll}\begin{array}{l}(K+\mathrm{1})\frac{{(Kh)}^{-\gamma }-{\delta }^{-\gamma }}{\gamma }-\\ \hspace{1em}\frac{{\delta }^{\mathrm{1}-\gamma }-{(Kh)}^{\mathrm{1}-\gamma }}{h(\mathrm{1}-\gamma )}\end{array}& \\ \hspace{1em}\hspace{1em}\gamma \ne \mathrm{0}& \\ (K+\mathrm{1})\mathrm{l}\mathrm{n}(\frac{\delta }{Kh})-\frac{\delta -Kh}{h}& \\ \hspace{1em}\gamma =\mathrm{0}& \end{array}\right.\end{array}$$

情形4：，利用对称性，得

$$\begin{array}{l}{A}_{\mathit{i},\mathit{i}-\mathit{K}}=-\alpha (x_{i-K+\frac{\mathrm{1}}{\mathrm{2}}})l_K-\alpha (x_{i-K-\frac{\mathrm{1}}{\mathrm{2}}})r_K-\\ \hspace{1em}\hspace{1em}\alpha (x_i)(l_K+r_K)\end{array}$$

（15）

情形5：，的支集，有

$$A_{\mathit{i},\mathit{i}+\mathit{K}+\mathrm{1}}\approx -\alpha (x_{i+K+\frac{\mathrm{1}}{\mathrm{2}}})l_{K+\mathrm{1}}-\alpha (x_i)l_{K+\mathrm{1}}$$

（16）

其中

$$l_{K+\mathrm{1}}=\left\{\begin{array}{ll}\frac{K\delta {}^{-\gamma }}{\gamma }-\frac{K{}^{\mathrm{1}-\gamma }{h}^{-\gamma }}{\gamma }+\frac{{\delta }^{\mathrm{1}-\gamma }-{(Kh)}^{\mathrm{1}-\gamma }}{h(\mathrm{1}-\gamma )}& \\ \hspace{1em}\hspace{1em}\gamma \ne \mathrm{0}& \\ -K\mathrm{l}\mathrm{n}(\frac{\delta }{Kh})+\frac{\delta -Kh}{h}& \\ \hspace{1em}\hspace{1em}\gamma =\mathrm{0}& \end{array}\right.$$

情形六：.利用对称性，得

$$A_{\mathit{i},\mathit{i}-\mathit{K}-\mathrm{1}}\approx -\alpha (x_{i-K-\frac{\mathrm{1}}{\mathrm{2}}})l_{K+\mathrm{1}}-\alpha (x_i)l_{K+\mathrm{1}}$$

（17）

1.3　刚度矩阵的结构分析

从参考文献［8］可知，刚度矩阵$\mathit{A}$可以分解为

$$\begin{array}{l}\mathit{A}={d_{\mathrm{1}}}^{h,\delta }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathit{\alpha })+{d_{\mathrm{2}}}^{h,\delta }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{\alpha }}^{″})+T\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{\alpha }}^{-})+\\ \hspace{1em}{T}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{\alpha }}^{+})+\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathit{\alpha })(\mathit{T}+{\mathit{T}}^{\mathrm{T}})\end{array}$$

（18）

其中，${\mathit{T}}^{}={\mathit{L}}^{\mathrm{{\rm T}}}+{\mathit{R}}^{}$，

$$\begin{array}{ll}\mathit{\alpha }=[\alpha (x_{\mathrm{1}}),\alpha (x_{\mathrm{2}}),\cdots ,\alpha (x_{N-\mathrm{1}}{)]}^{\mathrm{T}}& \\ {\mathit{\alpha }}^{-}=[\alpha (x\frac{\mathrm{1}}{\mathrm{2}}),\alpha (x\frac{\mathrm{3}}{\mathrm{2}}),\cdots ,\alpha (x_{N-\frac{\mathrm{3}}{\mathrm{2}}}{)]}^{\mathrm{T}}& \\ {\mathit{\alpha }}^{\text{'}}=[{\alpha }^{\text{'}}(x_{\mathrm{1}}),{\alpha }^{\text{'}}(x_{\mathrm{2}}),\cdots ,{\alpha }^{\text{'}}(x_{N-\mathrm{1}}{)]}^{\mathrm{T}}& \\ {\mathit{\alpha }}^{+}=[\alpha (x_{\frac{\mathrm{3}}{\mathrm{2}}}),\alpha (x_{\frac{\mathrm{5}}{\mathrm{2}}}),\cdots ,\alpha (x_{N-\frac{\mathrm{1}}{\mathrm{2}}}{)]}^{\mathrm{T}}& \\ {d_{\mathrm{1}}}^{h,\delta }=\mathrm{4}(\frac{h^{-\gamma }-{\delta }^{-\gamma }}{\gamma }+\frac{h^{-\gamma }}{\mathrm{1}-\gamma })& \\ {d_{\mathrm{2}}}^{h,\delta }=(\frac{{\delta }^{\mathrm{2}-\gamma }-h^{\mathrm{2}-\gamma }}{\mathrm{2}-\gamma }+\frac{h^{\mathrm{2}-\gamma }}{\mathrm{3}-\gamma })& \end{array}$$

（19）

$$\begin{array}{l}\mathit{L}=\\ \left(\begin{array}{cccccccccc}\mathrm{0}& \mathrm{0}& \cdots & \cdots & \mathrm{0}& \mathrm{0}& \cdots & \cdots & \mathrm{0}& \mathrm{0}\\ -l_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \cdots & \mathrm{0}& \mathrm{0}& \ddots & \ddots & \ddots & \mathrm{0}\\ ⋮& -l_{\mathrm{1}}& \mathrm{0}& \cdots & \ddots & \mathrm{0}& \ddots & \ddots & \ddots & ⋮\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & \ddots & \ddots & ⋮\\ -l_{K+\mathrm{1}}& -l_K& \ddots & \cdots & \mathrm{0}& \mathrm{0}& \cdots & \ddots & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -l_{K+\mathrm{1}}& -l_K& \cdots & -l_{\mathrm{1}}& \mathrm{0}& \cdots & \cdots & \mathrm{0}& \mathrm{0}\\ ⋮& \ddots & \ddots & \ddots & ⋮& ⋮& \ddots & \ddots & ⋮& ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & ⋮& \ddots & \ddots & ⋮& ⋮\\ \mathrm{0}& \ddots & \ddots & \ddots & -l_{K+\mathrm{1}}& -l_K& \cdots & \cdots & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \cdots & \cdots & \mathrm{0}& -l_{K+\mathrm{1}}& \cdots & \cdots & -l_{\mathrm{1}}& \mathrm{0}\end{array}\right)\end{array}$$

（20）

$$\begin{array}{l}\mathit{R}=\\ \left(\begin{array}{cccccccccc}\mathrm{0}& \mathrm{0}& \cdots & \cdots & \mathrm{0}& \mathrm{0}& \cdots & \cdots & \mathrm{0}& \mathrm{0}\\ -r_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \cdots & \mathrm{0}& \mathrm{0}& \ddots & \ddots & \ddots & \mathrm{0}\\ ⋮& -r_{\mathrm{1}}& \mathrm{0}& \cdots & \ddots & \mathrm{0}& \ddots & \ddots & \ddots & ⋮\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & \ddots & \ddots & ⋮\\ -r_K& -r_{K-\mathrm{1}}& \ddots & \cdots & \mathrm{0}& \mathrm{0}& \cdots & \ddots & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -r_K& -r_{K-\mathrm{1}}& \cdots & -r_{\mathrm{1}}& \mathrm{0}& \cdots & \cdots & \mathrm{0}& \mathrm{0}\\ ⋮& \ddots & \ddots & \ddots & ⋮& ⋮& \ddots & \ddots & ⋮& ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & ⋮& \ddots & \ddots & ⋮& ⋮\\ \mathrm{0}& \ddots & \ddots & \ddots & -r_K& -r_{K-\mathrm{1}}& \cdots & \cdots & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \cdots & \cdots & \mathrm{0}& -r_K& \cdots & \cdots & -r_{\mathrm{1}}& \mathrm{0}\end{array}\right)\end{array}$$

（21）

由式（18）可知： 刚度矩阵$\mathit{A}$为Toeplitz-like矩

阵${}^{[\mathrm{17}-\mathrm{18}]}$，且存储$\mathit{A}$需要的存储量。刚度矩阵$\mathit{A}$与任意向量作乘法只需要的计算量。

2 预处理GMRES算法

本文利用快速配置法数值求解变系数非局部扩散模型最终归结为线性方程组式（6）的求解，其中系数矩阵$\mathit{A}$由式（18）所定义，右端项$\mathit{f}$由式（7）所定义。由于系数矩阵$\mathit{A}$为非对称矩阵，所以可以利用GMRES方法求解线性方程组式（6）。为了提高GMRES方法的收敛率，本节构造了系数矩阵$\mathit{A}$的Toeplitz及循环预处理子，提出了预处理GMRES算法。首先，本文给出几个基本定义：

定义1  ${}^{[\mathrm{14}-\mathrm{15}]}$称具有如下结构，矩阵${\mathit{T}}_n$为Toeplitz阵。

$${\mathit{T}}_n=\left[\begin{array}{ccccc}{t}_{\mathrm{0}}& t_{-\mathrm{1}}& \cdots & t_{-(n-\mathrm{2})}& t_{-(n-\mathrm{1})}\\ t_{\mathrm{1}}& t_{\mathrm{0}}& t_{-\mathrm{1}}& \cdots & t_{-(n-\mathrm{2})}\\ t_{\mathrm{2}}& t_{\mathrm{1}}& t_{\mathrm{0}}& \ddots & ⋮\\ ⋮& ⋮& \ddots & \ddots & t_{-\mathrm{1}}\\ t_{n-\mathrm{1}}& t_{n-\mathrm{2}}& \cdots & t_{\mathrm{1}}& t_{\mathrm{0}}\end{array}\right]$$

即.设表示定义在区间内所有以为周期的实值连续函数。对所有的，令

$$t_k=\frac{\mathrm{1}}{\mathrm{2}\pi }{\int }_{-\pi }^{\pi }f(x)e^{-ikx}\mathrm{d}xk=\mathrm{0},\pm \mathrm{1},\pm \mathrm{2}\cdots$$

（22）

为的Fourier函数，则称函数是Toeplitz矩阵${\mathit{T}}_n$生成的函数。

定义2  称具有如下结构，矩阵${\mathit{C}}_n$为循环矩阵。

$${\mathit{C}}_n=\left[\begin{array}{ccccc}{c}_{\mathrm{0}}& c_{n-\mathrm{1}}& \cdots & c_{\mathrm{2}}& c_{\mathrm{1}}\\ c_{\mathrm{1}}& t_{\mathrm{0}}& c_{n-\mathrm{1}}& \cdots & c_{\mathrm{2}}\\ c_{\mathrm{2}}& c_{\mathrm{1}}& c_{\mathrm{0}}& \ddots & ⋮\\ ⋮& ⋮& \ddots & \ddots & c_{n-\mathrm{1}}\\ c_{n-\mathrm{1}}& c_{n-\mathrm{2}}& \cdots & c_{\mathrm{1}}& c_{\mathrm{0}}\end{array}\right]$$

即且.循环矩阵可以通过Fourier矩阵实现对角化，即

$${\mathit{C}}_n={{\mathit{F}}_n}^{\mathrm{*}}{\wedge }_n{\mathit{F}}_n$$

（23）

其中，Fourier矩阵${\mathit{F}}_n$中的元素为

$$[{\mathit{F}}_n{]}_{j,k}=\frac{\mathrm{1}}{\sqrt[]{n}}{e}^{\frac{\mathrm{2}\pi \mathrm{ }jk}{n}}\mathrm{0}\le \mathrm{ }i,j\le \mathrm{ }n-\mathrm{1}$$

${\mathit{\Lambda }}_n$是对角矩阵且对角线元素为${\mathit{C}}_n$的特征值。

定义3  ${}^{[\mathrm{16}-\mathrm{17}]}$ 对所有的维的循环矩阵${\mathit{C}}_n$，使得

$$\mathrm{m}\mathrm{i}\mathrm{n}{‖{\mathit{C}}_n-\mathit{T}‖}_{\mathrm{1}}$$

和

$$\mathrm{m}\mathrm{i}\mathrm{n}{‖{\mathit{C}}_n-\mathit{T}‖}_{\mathrm{\infty }}$$

（24）

成立的矩阵${\mathit{C}}_n$，称为Strang循环预处理子，记为${\mathit{S}}_{}$，该预处理子表示为

$$\mathit{S}=[s_{k-l}]\hspace{1em}\mathrm{0}\le \mathrm{ }k,l\le \mathrm{ }n-\mathrm{1}$$

$s_j$为Strang预处理子对角线上的元素，且有

$$s_j=\left\{\begin{array}{l}{t}_j\hspace{1em}\hspace{1em}\mathrm{0}\le \mathrm{ }j\le \mathrm{ }⌊\frac{n}{\mathrm{2}}⌋\\ t_{j-n}\hspace{1em}⌊\frac{n}{\mathrm{2}}⌋\le \mathrm{ }j\le \mathrm{ }n-\mathrm{1}\\ s_{j+n}\hspace{1em}\mathrm{0}\le \mathrm{ }-j\le \mathrm{ }n-\mathrm{1}\end{array}\right.$$

（25）

利用GMRES方法求解如下与原线性方程组（6）式等价的预处理后线性方程组

$${\mathit{P}}^{-\mathrm{1}}\mathbf{A}\mathbf{u}={\mathit{P}}^{-\mathrm{1}}\mathit{f}$$

（26）

其中，$\mathit{P}$为$\mathit{A}$的预处理子，得到如下预处理GMRES（PGMRES）方法。

算法（PGMRES方法）给定系数矩阵$\mathit{A}$，右端向量$\mathit{f}$，预处理子$\mathit{P}$，初始向量${\mathit{u}}^{(\mathrm{0})}$和允许误差.本算法计算${\mathit{u}}^{(k)}\in {\mathit{C}}^{}$，使得$\frac{{‖{\mathit{r}}^{(k)}‖}_{\mathrm{2}}}{{‖{\mathit{r}}^{(\mathrm{0})}‖}_{\mathrm{2}}}\le \mathrm{ }\epsilon$，其中${\mathit{r}}^{(k)}={\mathit{f}}^{}-{\mathit{A}}^{}{\mathit{u}}^{(k)},{\mathit{r}}^{(\mathrm{0})}={\mathit{P}}^{-\mathrm{1}}{\left({\mathit{f}}^{}-{\mathit{A}}^{}{\mathit{u}}^{(\mathrm{0})}\right)}^{},\beta ={‖{\mathit{r}}^{(\mathrm{0})}‖}_{\mathrm{2}},$${\mathit{v}}_{\mathrm{1}}=\frac{{\mathit{r}}^{(\mathrm{0})}}{\beta };$

取 $k=\mathrm{1,2},......$

对Krylov子空间${\mathit{𝒦}}_k\left({\mathit{P}}^{-\mathrm{1}}\mathit{A},{\mathit{v}}_{\mathrm{1}}\right)=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{{\mathit{v}}_{\mathrm{1}},{\mathit{P}}^{-\mathrm{1}}\mathit{A}{\mathit{v}}_{\mathrm{1}},\cdots ,{\left({\mathit{P}}^{-\mathrm{1}}\mathit{A}\right)}^{k-\mathrm{1}}{\mathit{v}}_{\mathrm{1}}\right\},$

用Arnoldi分解计算${\mathit{V}}_k=\left[{\mathit{v}}_{\mathrm{1}},{\mathit{v}}_{\mathrm{2}},\cdots ,{\mathit{v}}_k\right],$

以及阶矩阵${\tilde{\mathit{H}}}_{k+\mathrm{1},k}$，

$$\begin{array}{l}{\tilde{\mathit{H}}}_{k+\mathrm{1},k}=\left[\begin{array}{c}{\mathit{H}}_k\\ {\mathit{\beta }}_k\mathrm{ }{e}_k^T\end{array}\right]\mathrm{ }\mathrm{ }{\beta }_k=h_{k+\mathrm{1},k}\\ \\ {\mathit{H}}_k=\left[\begin{array}{ccccc}{h}_{\mathrm{11}}& h_{\mathrm{12}}& \cdots & h_{\mathrm{1},k-\mathrm{1}}& h_{\mathrm{1},k}\\ h_{\mathrm{21}}& h_{\mathrm{22}}& h_{\mathrm{23}}& \cdots & h_{\mathrm{2}k}\\ \mathrm{0}& \ddots & \ddots & \ddots & ⋮\\ ⋮& \cdots & \ddots & \ddots & \ddots \\ \mathrm{0}& ⋮& ⋮& h_{k,k-\mathrm{1}}& h_{k,k}\end{array}\right]\end{array}$$

计算${\tilde{\mathit{H}}}_{k+\mathrm{1},k}$的QR分解：${\tilde{\mathit{H}}}_{k+\mathrm{1},k}=G^{\mathrm{{\rm T}}}\left[\begin{array}{c}{\mathit{R}}_k\\ \mathrm{0}\end{array}\right]$，

$${\mathit{G}}_{}={\mathit{G}}_{\mathit{k}}{\mathit{G}}_{\mathit{k}-\mathrm{1}}...{\mathit{G}}_{\mathrm{1}}{\mathit{G}}_i=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{I}}_{i-\mathrm{1}},\left[\begin{array}{cc}{c}_i& s_i\\ -s_i& c_i\end{array}\right],{\mathit{I}}_{k-i}){c_i}^{\mathrm{2}}+{s_i}^{\mathrm{2}}=\mathrm{1}$$，

${\mathit{R}}_k$为非奇异上三角矩阵；

按式（27）计算和得

如果 ，

求解${\mathit{R}}^{(k)}{\mathit{z}}^{(k)}={\mathit{t}}^{(k)}$得到${\mathit{z}}^{(k)}$，

${\mathit{x}}^{(k)}={\mathit{x}}^{(\mathrm{0})}+{\mathit{V}}^{(k)}{\mathit{z}}^{(k)};$停算。

结束

$$\left\{\begin{array}{l}{\tau }_{\mathrm{1}}\mathrm{ }=\beta \mathrm{ }{c}_{\mathrm{1}}\\ {\tau }_i\mathrm{ }={(-\mathrm{1})}^{i-\mathrm{1}}\mathrm{ }\beta \mathrm{ }{s}_{\mathrm{1}}{s}_{\mathrm{2}}\cdots \mathrm{ }{s}_{i-\mathrm{1}}{c}_i\\ {\rho }_k\mathrm{ }={(-\mathrm{1})}^k\mathrm{ }\beta \mathrm{ }{s}_{\mathrm{1}}{s}_{\mathrm{2}}\cdots \mathrm{ }{s}_k\\ t_k=({\tau }_{\mathrm{1}},{\tau }_{\mathrm{2}},\cdots ,{\tau }_k{)}^{\mathrm{T}}\end{array}\right.$$

（27）

在以上PGMRES算法中，每一步需要求解残量方程$\mathbf{P}\mathbf{z}=\mathit{r}$，所以在构造预处理时要求预处理子结构简单且易于求解。

线性方程组式（6）中系数矩阵$\mathit{A}$为特殊结构，因此，构造刚度矩阵$\mathit{A}$的如下Toeplitz预处理子：

$${\mathit{P}}_{\mathrm{1}}=({d_{\mathrm{1}}}^{h,\delta }\overline{\alpha }+{d_{\mathrm{2}}}^{h,\delta }\overline{{\alpha }^{\text{'}}})\mathit{I}+\overline{{\alpha }^{-}}\mathrm{ }\mathit{T}+\overline{{\alpha }^{+}}\mathrm{ }{\mathit{T}}^{\mathrm{T}}+\overline{\alpha }(\mathit{T}+{\mathit{T}}^{\mathrm{T}})$$

其中，

$$\begin{array}{ll}\overline{\alpha }=\frac{\alpha (x_{\mathrm{1}})+\alpha (x_{\mathrm{2}})+\cdots +\alpha (x_{N-\mathrm{1}})}{N-\mathrm{1}}& \\ \overline{{\alpha }^{\text{'}}}=\frac{{\alpha }^{\text{'}}(x_{\mathrm{1}})+{\alpha }^{\text{'}}(x_{\mathrm{2}})+\cdots +{\alpha }^{\text{'}}(x_{N-\mathrm{1}})}{N-\mathrm{1}}& \\ \overline{{\alpha }^{-}}=\frac{\alpha (x_{\frac{\mathrm{1}}{\mathrm{2}}})+\alpha (x\frac{\mathrm{3}}{\mathrm{2}})+\cdots +\alpha (x_{N-\frac{\mathrm{3}}{\mathrm{2}}})}{N-\mathrm{1}}& \\ \overline{{\alpha }^{+}}=\frac{\alpha (x_{\frac{\mathrm{3}}{\mathrm{2}}})+\alpha (x\frac{\mathrm{5}}{\mathrm{2}})+\cdots +\alpha (x_{N-\frac{\mathrm{1}}{\mathrm{2}}})}{N-\mathrm{1}}& \end{array}$$

由于${\mathit{P}}_{\mathrm{1}}$为Toeplitz矩阵，在PGMRES算法中，求解以${\mathit{P}}_{\mathrm{1}}$为系数矩阵的残量方程${\mathit{P}}_{\mathrm{1}}\mathit{z}=\mathit{r}$的计算量太大。为减少计算量，我们进一步用${\mathit{P}}_{\mathrm{1}}$的Strang循环近似矩阵逼近${\mathit{P}}_{\mathrm{1}}$，令$\mathit{C}$为Toeplitz矩阵$\mathit{T}$的Strang循环逼近矩阵，即$\mathit{T}\approx \mathit{C}$，则得到刚度矩阵$\mathit{A}$的如下预处理子

$${\mathit{P}}_{\mathrm{2}}=({d_{\mathrm{1}}}^{h,\delta }\overline{\alpha }+{d_{\mathrm{2}}}^{h,\delta }\overline{{\alpha }^{\text{'}}})\mathit{I}+\overline{{\alpha }^{-}}\mathrm{ }\mathit{C}+\overline{{\alpha }^{+}}\mathrm{ }{\mathit{C}}^{\mathrm{T}}+\overline{\alpha }(\mathit{C}+{\mathit{C}}^{\mathrm{T}})$$

由于预处理子${\mathit{P}}_{\mathrm{2}}$为循环矩阵，则${\mathit{P}}_{\mathrm{2}}$可以通过离散的傅里叶变换进行对角化${\mathit{P}}_{\mathrm{2}}={\mathit{F}}^{\mathrm{*}}\mathit{\Lambda }{\mathit{F}}_n$，所以只需要一次FFT和一次逆FFT变换可求出，PGMRES算法中残量方程${\mathit{P}}_{\mathrm{2}}\mathit{z}=\mathit{r}$的解$\mathit{z}={\mathit{F}}^{\mathrm{*}}{\mathit{\Lambda }}^{-\mathrm{1}}{\mathit{F}}_n\mathit{r}$，计算量仅为.

同时，在以上PGMRES算法中还涉及到Toeplitz矩阵与向量相乘，由于Toeplitz矩阵与向量的乘积可以嵌入到阶循环矩阵与向量相乘中，也可以通过FFT计算量得到，计算量仍是所以PGMRES算法每次迭代的计算量为${}^{[\mathrm{17}-\mathrm{18}]}$。

3 数值算例

通过数值实验来验证本文构造的刚度矩阵$\mathit{A}$的循环预处理子${\mathit{P}}_{\mathrm{2}}$的有效性，分别利用GMRES方法、预处理GMRES方法求解线性方程组式（6），其中系数矩阵和右端项分别由式（18）和式（7）定义。所有迭代法都以零向量作为初始向量，迭代终止条件为

$$\frac{{‖{\mathit{r}}^{(k)}‖}_{\mathrm{2}}}{{‖{\mathit{r}}^{(\mathrm{0})}‖}_{\mathrm{2}}}\le \mathrm{ }{\mathrm{10}}^{-\mathrm{6}}$$

式中：${\mathit{r}}^{(k)}$为第步迭代残向量，${\mathit{r}}^{(\mathrm{0})}$为初始残向量。

在数值实验中，取方程

$$\left\{\begin{array}{ll}\begin{array}{l}{\int }_{x-\delta }^{x+\delta }(\alpha ( x)+\alpha (y))\\ \sigma (x-y)(u(x)-u(y))\mathrm{d}y=f(x)\end{array}& \\ x\in (a,b)& \\ u(x)=g(x)& \\ x\in (a-\delta ,a]\bigcup [b,b+\delta ]& \end{array}\right.$$

的参数如下：扩散区域=，变系数，，选取精确解为=，从而可以得到离散点处右端项的值。

算例1 取= .

算例2 取= ， .

从表1和表2可以得出以下结论：三种算法算出来的数值近似解是相同的；预处理GMRES算法需要的迭代步数和计算时间比GMRES算法需要的迭代步数和计算时间少；${\mathit{P}}_{\mathrm{2}}$做预处理子时需要的计算时间在这三种算法中是最短的。

表1 算例1的数值结果

Tab.1 The numerical results for Example 1

网格步长	GMRES方法		预处理GMRES方法 PGMRES
	迭代步数 Iter	计算时间 CPU/s	预处理子${\mathit{P}}_1$		预处理子${\mathit{P}}_2$
			迭代步数 Iter	计算时间 CPU/s	迭代步数 Iter	计算时间 CPU/s
$2^{-6}$	430	0.07	350	0.05	240	0.02
$2^{-7}$	480	0.39	378	0.45	273	0.68
$2^{-8}$	496	1.80	383	1.31	305	0.36
$2^{-9}$	500	2.86	406	2.14	328	0.56
$2^{-10}$	527	38.28	410	30.56	339	7.24
$2^{-11}$	530	570.87	420	450.67	380	110.26

表2 算例2的数值结果

Tab.2 The numerical results for Example 2

网格步长	GMRES方法		预处理GMRES方法 PGMRES
	迭代步数 Iter	计算时间 CPU/s	预处理子${\mathit{P}}_1$		预处理子${\mathit{P}}_2$
			迭代步数 Iter	计算时间 CPU/s	迭代步数 Iter	计算时间 CPU/s
$2^{-6}$	540	0.37	430	0.27	270	0.08
$2^{-7}$	650	0.06	515	0.05	321	0.13
$2^{-8}$	758	1.68	603	1.46	317	0.33
$2^{-9}$	896	2.96	547	2.42	674	0.48
$2^{-10}$	1011	40.56	803	32.67	800	8.96
$2^{-11}$	1100	318.63	870	256.89	890	63.96

图1和图2当中的横坐标表示矩阵特征值的实部，纵坐标表示特征值的虚部。从这两个图中，发现预处理矩阵的特征值分布比原始矩阵的特征值分布集中，说明预处理矩阵的条件数比原始矩阵的条件数小，故预处理GMRES方法的迭代步数比GMRES方法的迭代步数少。

图1 算例1中N为300时, 矩阵特征值分布图

Fig.1 When N is 300 in Example 1, the eigenvalue distribution of the matrix

图2 算例2中N为300时, 矩阵特征值分布图

Fig.2 When N is 300 in Example 2, the eigenvalue distribution of the matrix

4 总结与展望

利用快速配置法数值离散变系数非局部扩散模型所得到的一类具有Toeplitz结构的线性方程组的求解，针对系数矩阵的特殊结构，构造了一个Toeplitz预处理子和循环预处理子，并将其应用于GMRES方法中，即得到预处理GMRES方法。最后利用数值实验验证了预处理GMRES方法的有效性。

非局部扩散模型的数值方法有待进一步研究，接下来可以考虑从数值离散方法入手，提出新的离散方法，保证线性系统更易求解。

作者贡献声明

冉育红：构思者及负责人，指导实验设计、数据分析，论文的修改。

李存吉：数值算例设计者和执行者，完成数据分析，论文初稿的写作。

殷俊锋：参与数值实验设计和算例实验结果分析。