基于核贪婪算法的微分方程数值求解——贪婪算法的微分方程数值求解

Numerical Solution of Differential Equations Based on Kernel Greedy Algorithm—Greedy Algorithm for Numerical Solution of Differential Equations

摘要: 基于再生核Hilbert空间的相关知识，再生核Hilbert空间中的函数可以表示用n个再生点构成的核函数的线性表示进行逼近。为了实现n个再生点的最优选取，文章借助Santin等人提出f-greedy、P-greedy、f/P-greedy等贪婪类算法，实现上述方法在微分方程数值求解的理论推导及算法实验。与传统的均匀取点相比，P-greedy和预正交自适应Fourier分解(POAFD)的贪婪算法逼近误差更优，但f-greedy和f/P-greedy效果还存在改进的空间。

Abstract: Based on the knowledge of real regenerated kernel Hilbert space, the function in the regenerated kernel Hilbert space can be approximated by a linear representation of the kernel function composed of n regenerated points. In order to achieve the optimal selection of n regeneration points, Santin et al. put forward a greedy algorithm such as F-greedy, P-greedy, and f/P-greedy to achieve the theoretical derivation and algorithm experiment of the above method in numerical solution of differential equations. Compared to the traditional method of uniformly taking points, the P-greedy and Pre-Orthogonal Adaptive Fourier Decomposition (POAFD) greedy algorithms provide better approximation errors, while the f-greedy and f/P-greedy algorithms do not perform as well and have room for improvement.

1. 引言

设二阶线性微分方程：

$u^{″}\left(x\right)+a\left(x\right)u^{\prime }\left(x\right)+b\left(x\right)u\left(x\right)=f\left(x\right),x\in \left(0,1\right)$ (1)

满足边值条件：

$u\left(0\right)=u_0,u\left(1\right)=u_1$ (2)

其解存在且唯一，其中$a\left(x\right)$ ，$b\left(x\right)$ 是在区间[0, 1]上具有连续的导函数，$f\left(x\right)$ 为[0, 1]上的绝对连续函数。首先，我们可以通过线性平移的方式齐次化，将上述微分方程转化为以下的一般形式：

$\{\begin{array}{l}{u}^{″}\left(x\right)+a\left(x\right)u^{\prime }\left(x\right)+b\left(x\right)u\left(x\right)=f\left(x\right),x\in \left(0,1\right)\hfill \\ u\left(0\right)=0,u\left(1\right)=0\hfill \end{array}$ (3)

其中$f\left(x\right)$ ，$a\left(x\right)$ ，$b\left(x\right)$ 为已知函数。为了利用再生核Hilbert空间的知识，求解微分方程的数值解，我们假设$u\in W_3^2\left[0,1\right]$ ，$f\in W_2^1\left[0,1\right]$ ，关于再生核Hilbert空间$W_2^1\left[0,1\right]$ 及$W_3^2\left[0,1\right]$ 的具体介绍参见本文第2节。由于大多数偏微分方程的解析解较难得到，研究者们已经提供了不同的数值方法用于求解偏微分方程的近似解。

在1999年，L. Lapidus及G. F. Pinder提出加权残差法用于求解科学与工程中的偏微分方程数值解[1]；在2009年，C. Johnson等人提出了有限元方法求偏微分方程数值解[2]；在2017年，D. Watson等人提出径向基函数微分求积法求解微分方程[3]；同年，M. N. Özisik等学者提出有限差分方法(FDMs)求解热传导偏微分方程的数值解[4]；在2020年，D. Sharma、K. Goyal和R. K. Singla提出小波法求解偏微分方程数值解[5]。

利用再生核方法求解微分方程的历史也很悠久了。在2009年，崔明根等学者提出再生核方法求解偏微分方程[6]。后来有很多学者对此方法进行不断改进及推广应用，比如：在2017年，S. Abbasbandy，R. A. Van Gorder和P. Bakhtiari使用再生核方法求解Brinkman-Forchheimer动量方程[7]；在2018年，P. Bakhtiari、S. Abbasbandy和R. A. Van Gorder使用再生核方法求解一维Swift-Hohenberg方程[8]；在2020年，牛军等学者使用再生核方法求解自热传导方程[9]。

近年来，随着稀疏表示理论的兴起，钱涛教授等提出了自适应傅里叶分解(AFD)方法[10]，预正交自适应傅里叶分解方法(简称POAFD) [11]和及其他几种变体[12]。此方法利用能量下降最大原理，可以实现函数的快速逼近。基于核贪婪选点的思想，它被进一步应用于求解偏微分方程数值解等领域[13]-[15]。

基于以往学者关于再生核求解偏微分方程的工作，本文希望进一步借助Santin等人提出的最新P-greedy、f-greedy、f/P-greedy及预正交自适应Fourier分解等贪婪类算法策略，实现上述方法在微分方程数值求解的理论推导及算法实验。传统的均匀取点方法效率较低，而POAFD方法虽然有效，但计算复杂度较高。核贪婪算法为中心点选取提供了一种新的思路，具有一定的创新性。

2. 再生核的定义和性质

设$\text{Ω}$ 为一个非空集合，${\mathscr{H}}_k\left(\Omega \right)$ 是一个Hilbert函数空间，内积定义为$\langle \cdot ,\cdot \rangle$ ，函数$K:\Omega \times \Omega \to R$ 是空间${\mathscr{H}}_k\left(\Omega \right)$ 的再生核，当且仅当：

a) $K\left(\cdot ,x\right)\in {\mathscr{H}}_k\left(\Omega \right),\forall x\in \Omega$ ；

b) $f\left(x\right)={\langle f,K\left(\cdot ,x\right)\rangle }_{{\mathscr{H}}_k\left(\Omega \right)},\forall x\in \Omega ,\forall f\in {\mathscr{H}}_k\left(\Omega \right)$ 。

最后一个条件称为再生性：函数f在x点的值可以通过f和再生核$K\left(\cdot ,x\right)$ 的内积再生，具有再生核的Hilbert空间就称为再生核Hilbert空间，简称再生核空间[6]。

下面介绍再生核Hilbert空间的部分性质：

唯一性。如果Hilbert函数空间${\mathscr{H}}_k\left(\Omega \right)$ 有再生核，那么此再生核是唯一的。

设$K\left(x,y\right)$ 是Hilbert函数空间${\mathscr{H}}_k\left(\Omega \right)$ 的再生核，则$K\left(x,x\right)\ge 0$ 且$K\left(x,x\right)=0$ ，当且仅当${\mathscr{H}}_k\left(\Omega \right)=0$ 。

再生核是正定的，也就是对任意的$n\ge 1$ ，任意的$\left(a_1,\cdots ,a_n\right)\in R$ ，任意的$\left(x_1,\cdots ,x_n\right)\in \Omega$ ，总有${\displaystyle \sum _{i,j=1}^n{a}_i}{a}_{j}K\left(x_i,x_j\right)\ge 0$ 。

设$K\left(x,y\right)$ 是Hilbert函数空间${\mathscr{H}}_k\left(\Omega \right)$ 的再生核，V是${\mathscr{H}}_k\left(\Omega \right)$ 的闭子空间，则V是一个再生核空间，且其再生核为$K_V\left(x,y\right)=\left[P_V\left(K\left(\cdot ,y\right)\right)\right]\left(x\right)$ ，这里$P_V$ 表示${\mathscr{H}}_k\left(\Omega \right)$ 到V上的正交投影。

2.1. 再生核空间$W_2^1\left[0,1\right]$

$W_2^1\left[0,1\right]=\left\{u\left(x\right)|u\left(x\right)是�����值函�,u^{\prime }\in L^2\left[0,1\right],u\left(0\right)=0,u\left(1\right)=0\right\}$ 。对任意的

$u\left(x\right),v\left(x\right)\in W_2^1\left[0,1\right]$ ，其上定义的内积和范数为：

${\langle u\left(x\right),v\left(x\right)\rangle }_{W_1}=u\left(0\right)v\left(0\right)+{\displaystyle {\int }_0^1{u}^{\prime }}\left(x\right)v^{\prime }\left(x\right)\text{d}x$ (4)

${\Vert u\Vert }_{W_2^1}=\sqrt{{\langle u,u\rangle }_{W_2^1}}$

文献[16]证明了其是一个再生核Hilbert空间并且再生核为：

$K_x^{\left\{1\right\}}\left(y\right)=\{\begin{array}{l}1+y,y\le x\hfill \\ 1+x,y>x\hfill \end{array}$ (5)

2.2. 再生核空间$W_3^2\left[0,1\right]$

$W_3^2\left[0,1\right]=\left\{u\left(x\right)|u\left(x\right),u^{\prime }\left(x\right),u^{″}\left(x\right)是�����值函�,u^{(3)}\in L^2\left[0,1\right],u\left(0\right)=0,u\left(1\right)=0\right\}$ 对任意的$u\left(x\right),v\left(x\right)\in W_2^3\left[0,1\right]$ 其上定义的内积和范数为：

${\langle u\left(x\right),v\left(x\right)\rangle }_{W_3^2}=u\left(0\right)v\left(0\right)+u^{\prime }\left(0\right)v^{\prime }\left(0\right)+u^{″}\left(0\right)v^{″}\left(0\right)+{\displaystyle {\int }_0^1{u}^{\left(3\right)}}\left(x\right)v^{\left(3\right)}\left(x\right)\text{d}x$ (6)

${\Vert u\Vert }_{W_2^3}=\sqrt{{\left(u,u\right)}_{W_2^3}}$

文献[17]证明了$W_3^2\left[0,1\right]$ 是一个再生核Hilbert空间并且其再生核为：

$K_x^{\left\{3\right\}}\left(y\right)=\{\begin{array}{l}{\displaystyle \sum _{i=0}^5{c}_i}{y}^i,y\le x\hfill \\ {\displaystyle \sum _{i=0}^5{d}_i}{y}^i,y>x\hfill \end{array}$ (7)

其中${\left\{c_i\right\}}_{i=0}^5$ 和${\left\{d_i\right\}}_{i=0}^5$ 的具体表达式在文献[16]中找到。

3. 贪婪策略与稀疏表示

二阶线性微分方程(3)可以抽象表示为：

$\left(Lu\right)=u^{″}\left(x\right)+a\left(x\right)u^{\prime }\left(x\right)+b\left(x\right)u\left(x\right),u\left(x\right)\in W_3^2\left[0,1\right]$ (8)

定义一个算子$L:W_3^2\left[0,1\right]\to W_2^1\left[0,1\right]$ ，则在再生核空间中，问题(1.2)等价于求解算子方程：

$Lu=f$ (9)

因为方程(3)的解存在且唯一，则线性算子$L:W_3^2\left[0,1\right]\to W_2^1\left[0,1\right]$ 是有界且可逆的。

下面介绍有关线性算子L的性质。

引理1. 设$L^{*}$ 是L的共轭算子，则

$\left(L_x^{*}{K}_x^{\left\{1\right\}}\right)\left(y\right)=\left(L_x{K}_x^{\left\{3\right\}}\right)\left(y\right)$ (10)

证明：

$\begin{array}{c}{\psi }_i\left(x_i\right)=L_x^{*}{K}_x^{\left\{1\right\}}\left(x_i\right)\\ ={\langle L_x^{*}{K}_x^{\left\{1\right\}}\left(s\right),K_{x_i}^{\left\{3\right\}}\left(s\right)\rangle }_{W_2^3}\\ ={\langle K_x^{\left\{1\right\}}\left(s\right),L_x{K}_{x_i}^{\left\{3\right\}}\left(s\right)\rangle }_{W_2^1}\\ =L_x{K}_x^{\left\{3\right\}}\left(x_i\right)\end{array}$

引理2. 设${\left\{x_i\right\}}_{i=1}^{\infty }$ 是区间[0, 1]的一组稠密的互异点集，${\psi }_i\left(x\right)=\left(L_x{K}_x^{\left\{3\right\}}\right)\left(x_i\right),i=1,2,\cdots$ 。则${\left\{{\psi }_i\left(x\right)\right\}}_{i=1}^{\infty }$ 是线性无关的，并且是$W_3^2\left[0,1\right]$ 中一组完全系。

关于上述引理文献[17]已给出详细证明。

得到方程(3)的近似解，首先，我们对${\left\{{\psi }_i\left(x\right)\right\}}_{i=1}^{\infty }$ 进行施密特正交化得到再生核$W_3^2\left[0,1\right]$ 空间的一组标准正交基${\left\{{\tilde{\psi }}_i\left(x\right)\right\}}_{i=1}^{\infty }$ 。设$X_{\text{N}}=\left\{x_1,\cdots ,x_{\text{N}}\right\}\subset \left[0,1\right]$ 为一组数目为N的离散点集。设${\psi }_i\left(x\right)=\left(L_x{K}_x^{\left\{3\right\}}\right)\left(x_i\right)$ ，其中$i=1,\cdots ,\text{N}$ ，则其对应的正交基空间${\left\{{\tilde{\psi }}_i\left(x\right)\right\}}_{i=1}^{\text{N}}$可通过如下Gram-Schmidt过程获取，具体正交化的流程如下：

$\begin{array}{l}{\tilde{\psi }}_1\left(x\right):=\frac{{\psi }_1\left(x\right)}{\sqrt{{\langle {\psi }_1\left(x\right),{\psi }_1\left(x\right)\rangle }_{W_3^2}}}=\frac{{\psi }_1\left(x\right)}{\sqrt{{L_x{L}_y{K}_x^{\left\{3\right\}}\left(y\right)|}_{x=x_1}^{y=x_1}}},\\ {\widehat{\psi }}_k\left(x\right):={\psi }_k\left(x\right)-{\displaystyle \sum _{j=1}^{k-1}\langle {\psi }_j\left(x\right),{\tilde{\psi }}_j\left(x\right)\rangle {\tilde{\psi }}_j}\left(x\right),\end{array}$

${\tilde{\psi }}_k\left(x\right):=\frac{{\widehat{\psi }}_k\left(x\right)}{\sqrt{{\langle {\widehat{\psi }}_k\left(x\right),{\widehat{\psi }}_k\left(x\right)\rangle }_{W_3^2}}},1<k\le \text{N}.$ (11)

事实上，上述内积系数规范整理一下，可进一步表示为：

${\tilde{\psi }}_i\left(x\right)={\displaystyle \sum _{k=1}^i{\beta }_{ik}}{\psi }_k\left(x\right)$ (12)

其中${\beta }_{ik}$ 为正交化系数及其系数构成的矩阵$B={\left({\beta }_{ij}\right)}_{i,j=1}^N$ 满足：

$K_{\text{New}}=B^{-1}{\left(B^{\text{T}}\right)}^{-1}={\left(B^{\text{T}}B\right)}^{-1}$ (13)

而$K_{\text{new}}={\left({L_x{L}_y{K}_x^3\left(y\right)|}_{x=x_i}^{y=y_j}\right)}_{i,j=1}^N$ 为对再生核进行算子运算形成的矩阵。

引理3. 设$u\left((x)\right)$ 是方程(3)的解及${\left\{x_i\right\}}_{i=1}^{\infty }$ 是[0, 1]的一组稠密的互异点集，则：

$u\left(x\right)={\displaystyle \sum _{i=1}^{\infty }{\displaystyle \sum _{k=1}^i{\beta }_{ik}}}f\left(x_k\right){\tilde{\psi }}_i\left(x\right)$ (14)

证明：根据引理2，知${\left\{{\tilde{\psi }}_i\left(x\right)\right\}}_{i=1}^{\infty }$ 构成$W_3^2\left[0,1\right]$ 空间中的一组正交基。因此对于$u\left(x\right)\in W_3^2\left[0,1\right]$ 可以展成如下正交级数：

$\begin{array}{c}u\left(x\right)={\displaystyle \sum _{i=1}^{\infty }\langle u,{\tilde{\psi }}_i\rangle {\tilde{\psi }}_i\left(x\right)}\\ ={\displaystyle \sum _{i=1}^{\infty }\langle u,{\displaystyle \sum _{k=1}^i{\beta }_{ik}}{\psi }_k\left(x\right)\rangle {\tilde{\psi }}_i\left(x\right)}\\ ={\displaystyle \sum _{i=1}^{\infty }{\displaystyle \sum _{k=1}^i{\beta }_{ik}}}\langle u,{\psi }_k\rangle {\tilde{\psi }}_i\left(x\right)\\ ={\displaystyle \sum _{i=1}^{\infty }{\displaystyle \sum _{k=1}^i{\beta }_{ik}}}\langle u,L^{*}{K}_{x_k}^{\left\{1\right\}}\rangle {\tilde{\psi }}_i\left(x\right)\\ ={\displaystyle \sum _{i=1}^{\infty }{\displaystyle \sum _{k=1}^i{\beta }_{ik}}}\langle Lu,K_{x_k}^{\left\{1\right\}}\rangle {\tilde{\psi }}_i\left(x\right)\\ ={\displaystyle \sum _{i=1}^{\infty }{\displaystyle \sum _{k=1}^i{\beta }_{ik}}}f\left(x_k\right){\tilde{\psi }}_i\left(x\right).\end{array}$

对上面的正交级数作n项截断，便得到方程(3)的n项近似解：

$u_n\left(x\right)={\displaystyle \sum _{i=1}^n{\displaystyle \sum _{k=1}^i{\beta }_{ik}}}f\left(x_k\right){\tilde{\psi }}_i\left(x\right)$ (15)

优化选点策略

对于如何优化选取参数点$x_k$ 的问题，我们可以直接从左边利用文献[18]等人提出的预正交自适应傅里叶分解(简称POAFD)策略，在区间[0, 1]通过离散化得到N个点$X_{\text{N}}$ ，从N个点集中$X_{\text{N}}=\left(x_1,x_2,\cdots ,x_{\text{N}}\right)$

逐步优化选取一组最优子集$X_{I_{n-1}}=\left\{x_{i_1},x_{i_2},\cdots ,x_{i_{n-1}}\right\}\subset X_{\text{N}}$ ：

$\text{POAFD}:x_{i_n}=\underset{x\in X_{\text{N}}\backslash X_{I_{n-1}}}{\arg \max }\left|\langle u,{\tilde{\psi }}_i\rangle \right|$(16)

由于$L\left(u_n-u\right)=f_n-f$ ，因此希望从右边先寻找到最优再生点构成的函数快速逼近f。因为$f\in W_2^1\left[0,1\right]$ ，而$W_2^1\left[0,1\right]$ 空间也为再生核空间，其再生核$K_x^{\left\{1\right\}}\left(y\right)$ 在离散点取值构成的线性表示可以快速逼近函数值f，因此本文希望借助于Santin等人提出f-greedy、P-greedy、f/P-greedy贪婪策略从右端选点${\left\{x_{i_k}\right\}}_{k=1}^n$ ，再将所选出来的点代入左端基函数$\left(L_x{K}_x^{\left\{3\right\}}\right)\left(x_{i_k}\right)$ 进行Gram-Schmidt得到$W_3^2\left[0,1\right]$ 空间的标准正交基${\left\{{\tilde{\psi }}_{i_k}\left(x\right)\right\}}_{k=1}^n$ ，再根据引理3求得近似解$u_n(x)$ 。

接下来我们介绍f-greedy、P-greedy、f/P-greedy算法对函数f的选点策略。根据文献[18]和文献[19]，对于$f\in W_2^1\left[0,1\right]$ ，由于$K_x^{\left\{1\right\}}\left(y\right)$ 是再生核空间$W_2^1\left[0,1\right]$ 的再生核，将其离散化，得到一组数目为N的离散系${\left\{K_x^{\left\{1\right\}}\left(x_i\right)\right\}}_{i=1}^{\text{N}}$ 。类似于(11)，我们可以得到$W_2^1\left[0,1\right]$ 空间中的一组新的正交基：

$\begin{array}{l}{v}_1\left(x\right):=\frac{K_x^{\left\{1\right\}}\left(x_1\right)}{\sqrt{{\langle K_x^{\left\{1\right\}}\left(x_1\right),K_x^{\left\{1\right\}}\left(x_1\right)\rangle }_{W_2^1}}}=\frac{K_x^{\left\{1\right\}}\left(x_1\right)}{\sqrt{K_{x_1}^{\left\{1\right\}}\left(x_1\right)}},\\ {\widehat{v}}_k\left(x\right):=K_x^{\left\{1\right\}}\left(x_k\right)-{\displaystyle \sum _{j=1}^{k-1}{v}_j}\left(x_k\right)v_j\left(x\right),\end{array}$

$v_k\left(x\right):=\frac{{\widehat{v}}_k\left(x\right)}{\sqrt{K_x^{\left\{1\right\}}\left(x_k\right)-{\displaystyle \sum _{j=1}^{k-1}{v}_j}\left(x_k\right)v_j\left(x\right)}},1<k\le \text{N}.$(17)

因此，对于$f\in W_2^1\left[0,1\right]$ ，我们可以直接用N项正交展开函数$\left(s_{\text{N}}f\right)\left(x\right)={\displaystyle \sum _{k=1}^{\text{N}}\langle f,v_k\rangle v_k\left(x\right)}$ 逼近$f\left(x\right)$ 。为了实现对离散点集$\left(X_{\text{N}},Y_{\text{N}}\right)$ 的优化选取，根据文献[20]，N项逼近函数满足如下误差界：

$\left|f\left(x\right)-\left(s_{\text{N}}f\right)\left(x\right)\right|\le P_{\text{N}}\left(x\right){\Vert f\Vert }_{{\mathscr{H}}_k\left(\Omega \right)}$ (18)

其中$X_{\text{N}}=\left(x_1,x_2,\cdots ,x_{\text{N}}\right)$ ，$Y_{\text{N}}=\left(f\left(x_1\right),\cdots ,f\left(x_{\text{N}}\right)\right)$ ，而能量幂函数$P_{\text{N}}\left(x\right)$ 表示为：

$P_{\text{N}}{\left(x\right)}^2:=K_x^{\left\{1\right\}}\left(x\right)-{\displaystyle \sum _{j=1}^{\text{N}}{v}_j}{\left(x\right)}^2$ (19)

为了从N个点集$X_{\text{N}}=\left(x_1,x_2,\cdots ,x_{\text{N}}\right)$ 选取一组最优子集$X_{I_n}=\left\{x_{i_1},x_{i_2},\cdots .x_{i_n}\right\}\subset X_{\text{N}}$ ，实现核函数的贪婪分解策略，我们根据文献[21]提出的如下贪婪策略进行逐步选点：

$\begin{array}{l}\text{P-greedy}:\underset{\text{ x}\in X_{\text{N}}\backslash X_{I_{n-1}}}{x_{i_n}=\arg \max }{P}_{n-1}\left(x\right)\\ \text{f-greedy}:x_{i_n}=\underset{\text{x}\in X_{\text{N}}\backslash X_{I_{n-1}}}{\arg \max }\left|y_i-\left(s_{\text{N}}f\right)\left(x\right)\right|\\ \text{fP-greedy}:x_{i_n}=\underset{\text{x}\in X_{\text{N}}\backslash X_{I_{n-1}}}{\arg \max }\frac{\left|y_i-\left(s_{\text{N}}f\right)\left(x\right)\right|}{P_n\left(x\right)}\end{array}$ (20)

其中$X_{I_{n-1}}=\left\{x_{i_1},x_{i_2},\cdots ,x_{i_{n-1}}\right\}\subset X_{\text{N}}$ 表示已经根据上述条件选择$n-1$ 个最优点形成的点集。

根据上述贪婪策略，$\left(s_{\text{N}}f\right)\left(x\right)$ 的迭代更新可以通过下面的定理1实现，而$P_{n-1}\left(x\right)$ 的迭代更新通过式(18)。事实上，每一种贪婪策略的目的都不一样：P-greedy的目的是使函数最小，根据式子(17)进而减小误差上界。f-greedy的目的是找到残差最大的点，达到非插值点对应残差最小，进而减少误差。f/P-greedy

的目的是找到${\Vert f\Vert }_{{\mathscr{H}}_k\left(\Omega \right)}$ 范数最小，然后缩小误差。

设$X_n=\left\{x_{i_1},x_{i_2},\cdots ,x_{i_n}\right\}$ 是由P-greedy、f-greedy、f/P-greedy策略选取的n个最优子集，${\left\{v_{i_k}\left(x\right)\right\}}_{k=1}^n$ 是由Gram-Schmidt正交化$\left\{K_{x_{i_1}}^{\left\{1\right\}}\left(x\right),K_{x_{i_2}}^{\left\{1\right\}}\left(x\right),\cdots ,K_{x_{i_n}}^{\left\{1\right\}}\left(x\right)\right\}$ 得到的一组正交基。

令$\left(s_{n}f\right)\left(x\right)={\displaystyle \sum _{k=1}^n\langle f\left(x\right),v_{i_k}\left(x\right)\rangle v_{i_k}\left(x\right)}$ ，则我们有如下结论：

定理1. 设$X_{\text{N}}$ 为$\left[0,1\right]$ 区间上N个不同的离散点集。$X_n=\left\{x_{i_1},x_{i_2},\cdots ,x_{i_n}\right\}\subset X_{\text{N}}$ 是由P-Greedy、f-greedy、f/P-greedy策略选取的n个最优子集。则对于$\left(s_{n}f\right)\left(x\right)={\displaystyle \sum _{k=1}^n\langle f,v_{i_k}\left(x\right)\rangle v_{i_k}\left(x\right)}$ 也满足以下条件：

$\left(s_{n}f\right)\left(x\right)={\displaystyle \sum _{j=1}^n{\alpha }_j^{\ast }}{K}_x^{\left\{1\right\}}\left(x_{i_j}\right)={\displaystyle \sum _{k=1}^n{c}_k}{v}_k\left(x\right)=c_n{v}_n\left(x\right)+\left(s_{n-1}f\right)\left(x\right)$ ，

其中$c_n=\frac{f\left(x_{i_n}\right)-s_{n-1}\left(x_{i_n}\right)}{v_n\left(x_{i_n}\right)}$ ，${\alpha }^{*}:=K_{n\times n}^{-1}Y$ ，$Y={\left(f\left(x_{i_1}\right),\cdots ,f\left(x_{i_n}\right)\right)}^{\text{T}}$ ，$K_{n\times n}={\left(K_{x_{i_j}}^{\left\{1\right\}}\left(x_{i_k}\right)\right)}_{j,k=1}^n$ 为核矩阵，${\alpha }^{*}=\left({\alpha }_1,\cdots ,{\alpha }_n\right)$ 。

证明：根据最佳平方逼近性质：

${\Vert f\left(x\right)-{\displaystyle \sum _{k=1}^n\langle f\left(x\right),v_k\left(x\right)\rangle v_k\left(x\right)}\Vert }^2=\underset{{\alpha }_k^{*}}{\min }{\Vert f\left(x\right)-{\displaystyle \sum _{k=1}^n{\alpha }_k^{*}}{K}_x^{\left\{1\right\}}\left(x_{i_k}\right)\Vert }^2$ ，

而满足(3.11)最优值${\left\{{\alpha }_k^{*}\right\}}_{k=1}^n$ 使得$\left(s_{n}f\right)\left(x\right)={\displaystyle \sum _{k=1}^n{\alpha }_k^{\text{*}}}{K}_x^{\left\{1\right\}}\left(x_{i_k}\right)$ 可以通过求解以下线性方程组得到：

$K_{n\times n}{\alpha }_k^{*}=Y$ ，

其中根据(3.12)，知其也满足插值条件：$\left(s_{n}f\right)\left(x_{i_k}\right)=f\left(x_{i_k}\right)$ ，$\forall k=1,\cdots ,n$ ，${\alpha }^{*}:=K_{n\times n}^{-1}Y$ ， $Y={\left(f\left(x_{i_1}\right),\cdots ,f\left(x_{i_n}\right)\right)}^{\text{T}}$ ，$K_{n\times n}={\left(K_{x_{i_j}}^{\left\{1\right\}}\left(x_{i_k}\right)\right)}_{j,k=1}^n$ 为核矩阵，${\alpha }^{*}=\left({\alpha }_1,\cdots ,{\alpha }_n\right)$ 。

当对核矩阵$K_{n\times n}$ 进行cholesky分解$K_{n\times n}=V{V}^{\top }$ ，则$\left(s_{n}f\right)\left(x\right)$ 进一步有下面的矩阵变形式：

$\begin{array}{c}\left(s_{n}f\right)\left(x\right)={\displaystyle \sum _{k=1}^n\langle f\left(x\right),v_k\left(x\right)\rangle v_k\left(x\right)}\\ ={\displaystyle \sum _{k=1}^n{\alpha }_k^{*}}{K}_x^{\left\{1\right\}}\left(x_{i_k}\right)\\ =\left(K_x^{\left\{1\right\}}\left(x_{i_1}\right),\cdots ,K_x^{\left\{1\right\}}\left(x_{i_1}\right)\right)\left[\begin{array}{c}{\alpha }_1\\ ⋮\\ {\alpha }_n\end{array}\right]\\ =\left(K_x^{\left\{1\right\}}\left(x_{i_1}\right),\cdots ,K_x^{\left\{1\right\}}\left(x_{i_1}\right)\right)K_{n\times n}^{-1}\left[\begin{array}{c}f\left(x_{i_1}\right)\\ ⋮\\ f\left(x_{i_n}\right)\end{array}\right]\\ =\left(K_x^{\left\{1\right\}}\left(x_{i_1}\right),\cdots ,K_x^{\left\{1\right\}}\left(x_{i_1}\right)\right){\left(V^{\top }\right)}^{-1}{V}^{-1}\left[\begin{array}{c}f\left(x_{i_1}\right)\\ ⋮\\ f\left(x_{i_n}\right)\end{array}\right]\\ =\left(K_x^{\left\{1\right\}}\left(x_{i_1}\right),\cdots ,K_x^{\left\{1\right\}}\left(x_{i_1}\right)\right)B^{\top }B\left[\begin{array}{c}f\left(x_{i_1}\right)\\ ⋮\\ f\left(x_{i_n}\right)\end{array}\right]\\ =\left(v_1,\cdots ,v_n\right)\left[\begin{array}{c}{c}_1\\ ⋮\\ c_n\end{array}\right]\end{array}$

其中$V^{-1}=B$ 。因此，当令$\left(s_{n}f\right)\left(x\right)=\left(s_{n-1}f\right)\left(x\right)+c_n{v}_n$ 时，根据插值的性质，即：

$c_n=\frac{y\left(x_{i_n}\right)-s_{n-1}\left(x_{i_n}\right)}{v_n\left(x_{i_n}\right)}$ 。

4. 数值实验

在本节中，我们使用预正交自适应Fourier分解算法(POAFD)、f-greedy，P-greedy，f/P-greedy四种不同的贪婪算法，结合再生核方法求解偏微分方程问题。为此，我们首先考虑如下方程：

$\{\begin{array}{l}{u}^{″}\left(x\right)-2u^{\prime }\left(x\right)+u\left(x\right)=3x-1,x\in \left(0,1\right)\hfill \\ u\left(0\right)=0,u\left(1\right)=0\hfill \end{array}$ $$ (21)

方程的真解是，其图像如图1所示：

Figure 1. The true solution of the equation

图1. 方程的真实解图像

为了比较不同贪婪选点策略在偏微分方程数值逼近上的优势，将其与等间距取点方式进行比较，图2~5分别展示了n = 5时，均匀取点与不同贪婪选点策略n = 5上稀疏点的分布与近似解的逼近情况。

Figure 2. n = 5 uniform selection

图2. n = 5均匀取点

Figure 3. n = 5 P-greedy selection

图3. n = 5 P-greedy

Figure 4. n = 5 POAFD selection

图4. n = 5 POAFD

Figure 5. n = 5 f-greedy selection

图5. n = 5 f-greedy

Table 1. The approximation error of greedy point selection and uniform point selection in the kernel

表1. 贪婪选点核均匀选点的近似误差

n	POAFD	P-greedy	f-greedy	f\ P-greedy	均匀取点
1	$1.00\times {10}^{-3}$	$7.00\times {10}^{-3}$	$9.00\times {10}^{-3}$	$1.00\times {10}^{-3}$	$8.00\times {10}^{-3}$
4	$3.18\times {10}^{-5}$	$1.45\times {10}^{-5}$	$6.09\times {10}^{-5}$	$2.00\times {10}^{-4}$	$4.33\times {10}^{-5}$
6	$2.66\times {10}^{-6}$	$1.48\times {10}^{-6}$	$5.31\times {10}^{-5}$	$2.12\times {10}^{-5}$	$6.85\times {10}^{-6}$
8	$7.03\times {10}^{-7}$	$3.75\times {10}^{-7}$	$3.02\times {10}^{-6}$	$1.81\times {10}^{-5}$	$1.82\times {10}^{-6}$
10	$1.65\times {10}^{-7}$	$1.78\times {10}^{-7}$	$3.00\times {10}^{-6}$	$1.6\times {10}^{-5}$	$5.123\times {10}^{-7}$
12	$9.00\times {10}^{-8}$	$8.48\times {10}^{-8}$	$2.58\times {10}^{-6}$	$1.4\times {10}^{-5}$	$2.84\times {10}^{-7}$
14	$6.17\times {10}^{-8}$	$4.94\times {10}^{-8}$	$2.37\times {10}^{-6}$	$1.2\times {10}^{-5}$	$1.51\times {10}^{-7}$

实验结论

如表1所示，从实验结果分析得出，POAFD和P-greedy选取稀疏点所构造的截断误差比均匀取点的误差更接近方程的真实解，P-greedy的选点策略比POAFD效果更好，截断误差更小，收敛速度更快，然而f-greedy和f/P-greedy选取稀疏点所构造的截断误差比均匀取点的误差大，还存在改进的空间。