物理科学中的许多问题都可以归结为常微分方程的问题。常微分方程对人类生活、学习和工作的各个方面都有着很大的影响。实际上，常微分方程常用于模拟数学物理中的问题。因此，常微分方程的求解一直受到学者们高度关注。

然而，只有少数简单的常微分方程及其逆问题可以解析求解，在科学研究和社会生产中遇到的大多数常微分方程都是非常复杂的，很难得到它们的解析解。因此，考虑常微分方程及其逆问题的数值方法具有十分重要意义。逆问题由于不定性更是极其不易求解，这不仅是在数学理论上，而且在数值计算中都会有许多困难。

传统的求解微分方程的数值方法包括欧拉方法、龙格库塔法、有限差分等，这些方法可以得到满足某一精度的解。然而，它们需要对整个区域进行插值，这会增加计算复杂度，耗费较多的计算时间。

近年来，为了克服传统数值方法的缺陷，研究者们试图寻找更方便的方法来获得近似解析解，从而避免插值过程。利用神经网络求解常微分方程成为一个新的方向^[1]，此方法虽然可以得到具有解析表达式的近似解，但基于神经网络的常微分方程求解方法依赖于神经网络的函数逼近能力^[2]，将受限于它们的理论弱点，比如BP神经网络算法可能收敛不到局部最优解。

支持向量机(SVMs)由Vapnik和他的合作者在20世纪90年代提出，是建立在小样本统计学习理论和结构风险最小化原理基础上的通用机器学习算法^[3]。这个算法是在有限数据的情况下，寻找由置信区间和经验风险组成的结构风险最小化。另外，这个算法采用核函数可以非常容易地处理高维数问题。其基本思想是将原始问题通过非线性变换映射到高维特征空间，在高维空间中构造决策函数，从而将原问题转化为一个凸二次规划问题，进而得到唯一的全局最优解。因此，该算法能够有效地处理神经网络等经典算法中的维数问题、局部极小问题和过学习问题。支持向量机用于处理函数估计和分类问题，在模式识别等方面已得到广泛的应用^[4-10]。

Suykens等^[11]提出了最小二乘支持向量机(LS-SVMs)，用等式约束替代了支持向量机的不等式约束，避免了求解二次规划问题，运算速度快。文献[12-13]将最小二乘支持向量机用于求解微分方程，通过将解常微分问题转化为函数估计问题进行求解。然而，相比于支持向量机，最小二乘支持向量机虽然求解速度快，但其解不具有稀疏性，对噪声更敏感。

在许多物理现象中，观测数据是在离散的有限个点上获得的。因此对于离散微分方程及其逆问题的解的探究具有十分重要的意义。而上述所提到无论是传统的数值计算方法，还是基于神经网络及最小二乘支持向量机方法，在求解离散微分方程及其逆问题时都将失效。

在文献[14]中，Saitoh等通过利用Tikhonov正则化、迭代以及再生核理论^[15-17]给出了一种求解离散线性微分方程及其逆问题近似解的方法。该方法基于再生核理论，选定合适的再生核Hilbert空间，相应的算子就可以进行定义。该方法适用于不同的离散微分方程，且对于原问题以及逆问题都可以得到相同结构的解。然而，该方法需要通过迭代计算 $N$ 个再生核( $N$ 为给定离散点的个数)，当 $N$ 增大时，计算量也相应地增加，计算过程变得更加繁重且耗时，其算法的有效性严重依赖于所使用计算机的效率。

为克服文献[14]的缺陷，在Saitoh等工作的基础上^[14]，本文考虑结合支持向量机、Tikhonov正则化和再生核理论来求解离散线性微分方程及其逆问题。该方法基于再生核理论和Tikhonov正则化方法，将离散微分方程及其逆问题转化为支持向量机的函数逼近问题，最终转化为求解一个二次规划问题，进而求得具有解析表达式的稀疏解，具有较强的泛化能力，便于后续应用。对比于文献[14]中的方法，该方法具有文献[14]方法的优点，适用于不同的离散微分方程，原问题及其逆问题均具有相同结构的解，而且克服了需要迭代计算 $N$ 个再生核的缺点。在所提方法中，只需要所选择再生核Hilbert空间的再生核即可进行求解，计算复杂度大大减轻，节约大量计算时间。

1 初步内容

本文将提出一种新的方法来构造离散常微分方程及其相应的逆问题的近似解。所提方法适用于一般的离散微分方程，为便于说明，考虑文献[14]中所求解的微分方程：

$ LF = g $

(1)

其中， $L$ 是微分算子

$LF = F'' + \alpha F' + \beta F$

$F$ 属于某个再生核空间， $\alpha $ 和 $\beta $ 是给定的常数。

在许多物理现象中，观测数据是在离散的有限个点上获得的。因此，在本文中，笔者将在离散的有限个点上构造微分方程 $(1)$ 的近似解。离散微分方程问题：给定函数 $g$ 在某些离散点上的值，求解函数 $F$ 。离散微分方程的逆问题：给定函数 $F$ 在某些离散点上的值，求解函数 $g$ 。

对于 ${\mathbb R}$ 中的不同点 ${x_j}$ ，定义微分算子

$\begin{array}{l} {{L}}\left( F \right)\left( {{x_1},{x_2}, \ldots ,{x_N}} \right): = \\ \left( {\left( {LF} \right)\left( {{x_1}} \right),\left( {LF} \right)\left( {{x_2}} \right), \ldots ,\left( {LF} \right)\left( {{x_N}} \right)} \right) \in {{\mathbb R}^N}\\ \end{array} $

则 ${{L}}$ 是从一些再生核Hilbert空间到 ${{\mathbb R}^N}$ 空间的有界线性算子。

在本文中将考虑在两个典型的再生核Hilbert空间中求解微分方程及逆问题，这两个空间分别是Sobolev空间和Paley-Wiener空间。

对于实数 $s$ ， $\displaystyle s > \frac{1}{2}$ 的Sobolev空间 ${H_{{K_s}}}$ ，其再生核定义为

${K_s}\left( {x,y} \right) = \frac{1}{{2{\rm{\text{π}}} }}\int_{\mathbb R} {{{\left( {1 + {{\left| \zeta \right|}^2}} \right)}^{ - s}}{{\rm{e}}^{{\rm{i}}\left( {x - y} \right) \cdot \zeta }}{\rm{d}}\zeta } $

其中 $x,y \in {\mathbb R}$ 。

${H_{{K_s}}}$ 中的函数 $f$ 的范数定义为

${\left\| f \right\|_{{H_{{K_s}}}}} = {\left( {\int_{\mathbb R} {{{\left( {1 + {{\left| \zeta \right|}^2}} \right)}^s}{{\left| {\left( {{\cal F}f} \right)\left( \zeta \right)} \right|}^2}{\rm{d}}\zeta } } \right)^{\tfrac{1}{2}}}$

其中 ${\cal F}$ 是傅里叶变换^[16]。

笔者还考虑求解Paley-Wiener空间 ${H_{{K_h}}}$ 中的问题. Paley-Wiener空间 ${H_{{K_h}}}$ 中的函数 $f$ 是由以下积分变换得到

$f\left( z \right): = \frac{1}{{2{\text{π}} }}\int_{\mathbb R} {{\chi _h}\left( t \right)g\left( t \right){{\rm{e}}^{ - {\rm{i}}z \cdot t}}{\rm{d}}t}$

其中 ${\;\chi _h}\left( t \right)$ 是 ${\left( { - \text{π}/h, + \text{π}/h} \right)}$ 的特征函数， $g \in {L_2}{( { - \text{π}/h, }}$ $+ \text{π}/h)$ ， $h > 0$ 。

当 $z$ 趋向于 $\infty $ 时，在 ${H_{{K_h}}}$ 中的函数 $f$ 满足

$\left| {f\left( z \right)} \right| \leqslant C\;\exp\left( {\frac{{\text{π} \left| z \right|}}{h}} \right)$

和

${\int_{\mathbb R} {\left| {f\left( x \right)} \right|} ^2}{\rm{d}}x < \infty $

其中 $C$ 为某一常数。

${H_{{K_h}}}$ 空间的再生核为

$\begin{array}{l} \displaystyle {K_h}\left( {z,\bar u} \right) = \frac{1}{{2 \text{π} }}\int_{\mathbb R} {{\chi _h}\left( t \right){{\rm{e}}^{ - {\rm{i}}z \cdot t}}\overline {{{\rm{e}}^{ - {\rm{i}}u \cdot t}}} {\rm{d}}t} = \\ \displaystyle \frac{1}{{\text{π} \left( {z - \bar u} \right)}}\sin \frac{\text{π} }{h}\left( {z - \bar u} \right)\\ \end{array} $

其中 $z,u \in {\mathbb C}$ (参见文献[14])。

2 基于SVMs的离散微分方程求解方法

为了求解离散线性微分方程，在文献[14]中，作者们通过应用再生核理论和Tikhonov正则化，将问题表述如下：

对于任意 $N$ 个给定值 $\left\{ {g\left( {{x_j}} \right)} \right\}_{j = 1}^N$ ，求解 ${H_{{K_h}}}$ 中的函数 $F$ 等价于找到如下问题的解：

$ \mathop {\inf }\limits_{F \in {H_{{K_h}}}} \left\{ {\lambda \left\| F \right\|_{{H_{{K_h}}}}^2 + \sum\limits_{j = 1}^N {{{\left| {\left( {LF} \right)\left( {{x_j}} \right) - g\left( {{x_j}} \right)} \right|}^2}} } \right\} $

(2)

其中 $\lambda > 0$ 。

为了求解问题 $(2)$ ，文献[14]的算法需要从一个离散点开始，迭代计算 $N$ 个再生核。该算法的缺点是计算负担相对比较重，算法的有效性在很大程度上受正在使用的计算机效率的影响。

事实上，问题 $(2)$ 可以适应于学习理论中的更一般的框架，如下：

对于任意常数 $C > 0$ ，求解

$\mathop {\min}\limits_{F \in {H_{{K_h}}}} \left\{ {\left\| F \right\|_{{H_{{K_h}}}}^2 + \frac{C}{N}\sum\limits_{j = 1}^N {c\left( {\left( {LF} \right)\left( {{x_j}} \right),g\left( {{x_j}} \right)} \right)} } \right\}$

其中 $c\left( { \cdot , \cdot } \right)$ 是一个损失函数，用来度量模型的预测值与真实值之间的差异。问题 $(2)$ 可以看作是选取 $c\left( {\left( {LF} \right)\left( {{x_j}} \right),g\left( {{x_j}} \right)} \right) = {\left| {\left( {LF} \right)\left( {{x_j}} \right) - g\left( {{x_j}} \right)} \right|^2}$ 。不同的损失函数可以导出不同类型的解。

在学习理论的指导下，特别地，笔者选取

$ \begin{array}{l} c\left( {\left( {LF} \right)\left( {{x_j}} \right),g\left( {{x_j}} \right)} \right) = {\left| {\left( {LF} \right)\left( {{x_j}} \right) - g\left( {{x_j}} \right)} \right|_\varepsilon } = \\ \max\left\{ {0,\left| {\left( {LF} \right)\left( {{x_j}} \right) - g\left( {{x_j}} \right)} \right| - \varepsilon } \right\} \\ \end{array} $

其中 $\varepsilon $ 是一个选定的正数。这个损失函数被称为 $\varepsilon - $ 不敏感损失函数^[18]，它能够处理噪声破坏的数据并提供稀疏解，在应用阶段可以极大地简化解的计算。这种损失函数常被用于支持向量机。

接下来，笔者将在支持向量机的框架中，结合再生核理论和Tikhonov正则化方法进行求解。我们将重点讨论在Paley-Wiener空间 ${H_{{K_h}}}$ 中的求解算法，然后将该算法直接扩展到Sobolev空间 ${H_{{K_s}}}$ 中。

2.1 Paley-Wiener空间 ${H_{{K_h}}}$ 中的解法

对于任意 $N$ 个给定值 $\left\{ {g\left( {{x_j}} \right)} \right\}_{j = 1}^N$ ，参数 $C > 0$ ，求解

$ \mathop {\min}\limits_{F \in {H_{{K_h}}}} \left\{ {\left\| F \right\|_{{H_{{K_h}}}}^2 + \frac{C}{N}\sum\limits_{j = 1}^N {{{\left| {\left( {LF} \right)\left( {{x_j}} \right) - g\left( {{x_j}} \right)} \right|}_\varepsilon }} } \right\} $

(3)

对于不同点 $\left\{ {{x_j}} \right\}_{j = 1}^N$ ，设 ${{L}}$ 是 ${H_{{K_h}}}$ 到 ${{\mathbb R}^N}$ 空间的有界线性算子。即

$\begin{array}{l} {{L}}:F \in {H_{{K_h}}} \to \\ \left( {\left( {LF} \right)\left( {{x_1}} \right),\left( {LF} \right)\left( {{x_2}} \right), \cdots ,\left( {LF} \right)\left( {{x_N}} \right)} \right) \in {{\mathbb R}^N} \\ \end{array} $

在 ${{\mathbb R}^N}$ 空间中取一个标准正交系 $\left\{ {{e_j}} \right\}_{j = 1}^N$ ，有

$\left( {LF} \right)\left( {{x_j}} \right) = {\left\langle {{{L}}F,{e_j}} \right\rangle _{{{\mathbb R}^N}}} = {\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle _{{H_{{K_h}}}}}$

则线性映射 ${{L}}$ 的像空间是一个再生核Hilbert空间，记为 ${H_A}$ ，其再生核矩阵为^[14]

$ \begin{array}{l} A = {\left\{ {{{\left\langle {{{{L}}^ * }{e_j},{{{L}}^ * }{e_{j'}}} \right\rangle }_{{H_{{K_h}}}}}} \right\}_{j,j' = 1,2, \ldots ,N}} = \\ {\left\{ {{a_{jj'}}} \right\}_{j,j' = 1,2, \ldots ,N}} \end{array} $

(4)

其中

$\begin{array}{l} \left( {{{{L}}^ * }{e_j}} \right)\left( x \right) = {\left\langle {\left( {{{{L}}^ * }{e_j}} \right)\left( \cdot \right),{K_h}\left( { \cdot ,x} \right)} \right\rangle _{{H_{{K_h}}}}} = \\ {\left\langle {{e_j},{{L}}{K_h}\left( { \cdot ,x} \right)} \right\rangle _{{{\mathbb R}^N}}} = \\ \frac{1}{{2 \text{π} }}\int_{ - \text{π} /h}^{\text{π} /h} {\left( { - {\eta ^2} - {\rm{i}}\alpha \eta + \beta } \right)} {{\rm{e}}^{{\rm{i}}\eta \left( {x - {x_j}} \right)}}{\rm{d}}\eta \\ \end{array} $

定理1^[14] 　式(4)中引入的元素 ${a_{jj'}}$ 为

$ \begin{array}{l} \displaystyle {a_{jj'}} = \frac{{4X\left( { - 6 + {b^2}{X^2}} \right)\cos \left( {bX} \right)}}{{\text{π} {b^4}}} + \\ \displaystyle \frac{{\left( {24 - 12{b^2}{X^2} + {b^4}{X^4}} \right)\sin \left( {bX} \right)}}{{\text{π} {b^5}}} + \\ \displaystyle \frac{{\left( {{\alpha ^2} - 2\beta } \right)\left( {2 - {b^2}{X^2}} \right)\sin \left( {bX} \right)}}{{\text{π} {b^3}}} + \\ \displaystyle \frac{{ - 2\left( {{\alpha ^2} - 2\beta } \right)X\cos \left( {bX} \right)}}{{\text{π} {b^2}}} + \frac{{{\beta ^2}}}{{\text{π} b}}\sin \left( {bX} \right) \\ \end{array} $

其中 $X = \text{π}/h$ ， $b = \left| {{x_j} - {x_{j'}}} \right|$ ，矩阵 ${\left\{ {{a_{jj'}}} \right\}_{j,j' = 1,2, \cdots ,N}}$ 是正定矩阵。

基于上述结果，求解问题 $(3)$ 等价于求解

$ \mathop {\min}\limits_{F \in {H_{{K_h}}}} \left\{ {\frac{1}{2}\left\| F \right\|_{{H_{{K_h}}}}^2 + \frac{C}{N}\sum\limits_{j = 1}^N {{{\left| {{{\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle }_{{H_{{K_h}}}}} - g\left( {{x_j}} \right)} \right|}_\varepsilon }} } \right\} $

(5)

其中，添加因子 $\frac{1}{2}$ 是为了简化计算，不会影响求解。

为了处理 $\varepsilon - $ 不敏感损失函数，引入松弛变量 ${\xi ^{\left( * \right)}} = \left[ {\xi _1^{\left( * \right)},\xi _2^{\left( * \right)}, \cdots ,\xi _N^{\left( * \right)}} \right]$ ，其中 $\xi _j^{\left( * \right)}$ 代表 ${\xi _j}$ 和 $\xi _j^ * $ 。笔者将问题 $(5)$ 转化为以下可以具有相同解的等价问题：

$ \mathop {\min}\limits_{F \in {H_{{K_h}}},\xi ,\xi *} \left\{ {\frac{1}{2}\left\| F \right\|_{{H_{{K_h}}}}^2 + \frac{C}{N}\sum\limits_{j = 1}^N {\left( {{\xi _j} + \xi _j^ * } \right)} } \right\} $

(6)

其中

${\xi _j} = \max\left\{ {0,{{\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle }_{{H_{{K_h}}}}} - g\left( {{x_j}} \right) - \varepsilon } \right\}$

和

$\xi _j^ * = \max\left\{ {0,g\left( {{x_j}} \right) - {{\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle }_{{H_{{K_h}}}}} - \varepsilon } \right\}$

可以容易地证明

$\begin{array}{l} {\xi _j} + \xi _j^ * = \max\left\{ {0,\left| {{{\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle }_{{H_{{K_h}}}}} - g\left( {{x_j}} \right)} \right| - \varepsilon } \right\} = \\ {\left| {\left( {LF} \right)\left( {{x_j}} \right) - g\left( {{x_j}} \right)} \right|_\varepsilon } \\ \end{array} $

以及

$\begin{array}{c} {\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle _{{H_{{K_h}}}}} - g\left( {{x_j}} \right) \leqslant \varepsilon + {\xi _j},j = 1, \cdots ,N, \\ g\left( {{x_j}} \right) - {\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle _{{H_{{K_h}}}}} \leqslant \varepsilon + \xi _j^ * ,j = 1, \cdots ,N, \\ \xi _j^{\left( * \right)} \geqslant 0,j = 1, \cdots ,N \\ \end{array} $

问题 $(6)$ 可以看作是一个二次规划问题，将其称为原始问题。原始问题可以利用拉格朗日乘子通过对偶方法求解^[18]。为此，引入拉格朗日乘子来构造拉格朗日函数。

$\begin{array}{l} {\cal L}\left( {F,{\xi ^{\left( * \right)}},{\mu ^{\left( * \right)}},{\nu ^{\left( * \right)}}} \right) = \\ \displaystyle \frac{1}{2}{\left\| F \right\|^2} + \frac{C}{N}\sum\limits_{j = 1}^N {\left( {{\xi _j} + \xi _j^ * } \right)} \; - \\ \sum\limits_{j = 1}^N {{\mu _j}} \left[ {\varepsilon + {\xi _j} - {{\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle }_{{H_{{K_h}}}}} + g\left( {{x_j}} \right)} \right] - \\ \sum\limits_{j = 1}^N {\mu _j^ * } \left[ {\varepsilon + \xi _j^ * + {{\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle }_{{H_{{K_h}}}}} - g\left( {{x_j}} \right)} \right]\; - \\ \sum\limits_{j = 1}^N {\left( {{\nu _j}{\xi _j} + \nu _j^ * \xi _j^ * } \right)} \end{array} $

其中 ${\mu ^{\left( * \right)}} = \left[ {\mu _1^{\left( * \right)}, \ldots ,\mu _N^{\left( * \right)}} \right]$ 和 ${\nu ^{\left( * \right)}} = \left[ {\nu _1^{\left( * \right)}, \ldots ,\nu _N^{\left( * \right)}} \right]$ 是拉格朗日乘子，它们对应于原始函数的每个约束，而且必须满足正约束，如 $\mu _j^{\left( * \right)} \geqslant 0$ 和 $\nu _j^{\left( * \right)} \geqslant 0$ 。可以看出，原始问题 $(6)$ 等价于 ${\min _{F,{\xi ^{\left( * \right)}}}}{\max _{{\mu ^{\left( * \right)}},{\nu ^{\left( * \right)}}}}{\cal L}$ ，其对偶问题是 ${\max _{{\mu ^{\left( * \right)}}}}, $ $_{{{\nu ^{\left( * \right)}}}}{\min _{F,{\xi ^{\left( * \right)}}}}{\cal L}$ 。

如果原始问题是凸的，则一个解既是原始问题的最优解也是对偶问题的最优解当且仅当这个解满足Karush-Kuhn-Tucker (KKT)条件^[19]。显然，原始问题 $(6)$ 是凸的，相应的KKT条件为

$\begin{array}{c} {\mu _j}\left[ {\varepsilon + {\xi _j} - {{\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle }_{{H_{{K_h}}}}} + g\left( {{x_j}} \right)} \right] = 0 \\ \mu _j^ * \left[ {\varepsilon + \xi _j^ * + {{\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle }_{{H_{{K_h}}}}} - g\left( {{x_j}} \right)} \right] = 0 \\ {\nu _j}{\xi _j} = 0 \\ \nu _j^ * \xi _j^ * = 0 \end{array} $

和

$\begin{array}{c} \displaystyle \frac{{\partial {\cal L}}}{{\partial F}} = F - \sum\limits_{j = 1}^N {\left( {\mu _j^ * - {\mu _j}} \right)\left( {{{{L}}^ * }{e_j}} \right)} = 0 \\ \displaystyle \frac{{\partial {\cal L}}}{{\partial {\xi _j}}} = \frac{C}{N} - {\mu _j} - {\nu _j} = 0 \\ \displaystyle \frac{{\partial {\cal L}}}{{\partial \xi _j^ * }} = \frac{C}{N} - \mu _j^ * - \nu _j^ * = 0 \\ \end{array} $

对于 $j = 1, \cdots ,N$ 。

在这些条件下，原始问题 $(6)$ 可以通过其对偶问题进行求解，其对偶问题为

$\begin{array}{l} \mathop {\max }\limits_{{\mu ^{\left( * \right)}}} - \displaystyle \frac{1}{2}\sum\limits_{j,j' = 1}^N {\left( {\mu _j^ * - {\mu _j}} \right)} \left( {\mu _{j'}^ * - {\mu _{j'}}} \right)\left\langle {{{{L}}^ * }{e_j},{{{L}}^ * }{e_{j'}}} \right\rangle - \\ \displaystyle \varepsilon \sum\limits_{j = 1}^N {\left( {\mu _j^ * + {\mu _j}} \right)} + \sum\limits_{j = 1}^N {g\left( {{x_j}} \right)\left( {\mu _j^ * - {\mu _j}} \right)}\\ \end{array} $

对偶问题的约束条件为

$0 \leqslant \mu _j^{\left( * \right)} \leqslant \frac{C}{N},j = 1, \cdots ,N$

从而可得到解

${\overset{\frown} F} \left( x \right) = \sum\limits_{j = 1}^N {\left( {\mu _j^ * - {\mu _j}} \right)} \left( {{{{L}}^ * }{e_j}} \right)\left( x \right)$

实际上，通过设置适当的 $\varepsilon > 0$ ，只有拉格朗日乘子的一个子集将是非零的，这会导致产生一个稀疏解。因此，解可以写成

${\overset{\frown} F} \left( x \right) = \sum\limits_{j \in J} {\left( {\mu _j^ * - {\mu _j}} \right)} \left( {{{{L}}^ * }{e_j}} \right)\left( x \right)$

其中 $J$ 是关于 $\mu _j^ * $ 和 ${\mu _j}$ 不为零的输入数据点 ${x_j}$ 的指标集。这些 ${x_j}$ $\left( {j \in J} \right)$ 称为支持向量。(注释：参数 $C$ 是可以调整的，为了得到更优的解，可以调整参数 $C$ 的值。)

2.2 Sobolev空间 ${H_{{K_s}}}$ 中的解法

在光滑 $s\left( {s \geqslant 3} \right)$ 的Sobolev再生核Hilbert空间 ${H_{{K_s}}}$ 中，有

$\begin{array}{l} \left( {{{{L}}^ * }{e_j}} \right)\left( x \right) = {\left\langle {\left( {{{{L}}^ * }{e_j}} \right)\left( \cdot \right),{K_s}\left( { \cdot ,x} \right)} \right\rangle _{{H_{{K_s}}}}} = \\ {\left\langle {{e_j},{{L}}{K_s}\left( { \cdot ,x} \right)} \right\rangle _{{{\mathbb R}^N}}} = \\ \displaystyle \frac{1}{{2 \text{π} }}\int_{\mathbb R} {\frac{{\left( { - {\eta ^2} - {\rm{i}}\alpha \eta + \beta } \right){{\rm{e}}^{ - {\rm{i}}\eta \left( {x - {x_j}} \right)}}}}{{{{\left( {1 + {\eta ^2}} \right)}^s}}}} {\rm{d}}\eta \\ \end{array} $

则线性映射 ${{L}}$ 的像空间是一个再生核Hilbert空间，其再生核为^[14]

${a_{jj'}} = \frac{1}{{2 \text{π} }}\int_{\mathbb R} {\frac{{{{\left| { - {\eta ^2} - {\rm{i}}\alpha \eta + \beta } \right|}^2}{{\rm{e}}^{{\rm{i}}\eta \left( {{x_j} - {x_{j'}}} \right)}}}}{{{{\left( {1 + {\eta ^2}} \right)}^s}}}} {\rm{d}}\eta$

对于任意 $N$ 个给定值 $\left\{ {g\left( {{x_j}} \right)} \right\}_{j = 1}^N$ ，参数 $C > 0$ ，想要求解

$ \mathop {\min}\limits_{F \in {H_{{K_s}}}} \left\{ {\frac{1}{2}\left\| F \right\|_{{H_{{K_s}}}}^2 + \frac{C}{N}\sum\limits_{j = 1}^N {{{\left| {{{\left\langle {F,{{{L}}^ * }{e_j}} \right\rangle }_{{H_{{K_s}}}}} - g\left( {{x_j}} \right)} \right|}_\varepsilon }} } \right\} $

(7)

求解过程类似于Paley-Wiener空间 ${H_{{K_h}}}$ 的情况，可以得到问题 $(7)$ 的解

$\tilde F\left( x \right) = \sum\limits_{j \in J} {\left( {\upsilon _j^ * - {\upsilon _j}} \right)} \left( {{{{L}}^ * }{e_j}} \right)\left( x \right)$

其中 $\upsilon _j^ * \geqslant 0$ 和 ${\upsilon _j} \geqslant 0$ 是拉格朗日乘子。

3 基于SVMs的逆问题求解方法

在这一节中，将采用基于支持向量机的方法来考虑微分方程 $(1)$ 的逆问题。

那就是，对于微分方程

$LF = g$

如果给定函数 $F$ 在 $N$ 个点上的值 $\left\{ {F\left( {{x_j}} \right)} \right\}_{j = 1}^N$ ，笔者想要在某些再生核Hilbert空间中求解函数 $g$ 。

3.1 Paley-Wiener空间 ${H_{{K_h}}}$ 中的解法

首先，考虑在Paley-Wiener空间 ${H_{{K_h}}}$ 中的逆问题. 对于不同点 $\left\{ {{x_j}} \right\}_{j = 1}^N$ ，设 ${{M}}$ 是 ${H_{{K_h}}}$ 到 ${{\mathbb R}^N}$ 空间的有界线性算子。即

${{M}}:g \in {H_{{K_h}}} \to \left( {F\left( {{x_1}} \right),F\left( {{x_2}} \right), \ldots ,F\left( {{x_N}} \right)} \right) \in {{\mathbb R}^N}$

然后有

$F\left( {{x_j}} \right) = {\left\langle {{{M}}g,{e_j}} \right\rangle _{{{\mathbb R}^N}}} = {\left\langle {g,{{{M}}^ * }{e_j}} \right\rangle _{{H_{{K_h}}}}}$

则线性映射 ${{M}}$ 的像空间是一个再生核Hilbert空间，记为 ${H_B}$ ，其再生核矩阵为^[14]

$B = {\left\langle {{{{M}}^ * }{e_j},{{{M}}^ * }{e_{j'}}} \right\rangle _{{H_{{K_h}}}}} = {\left\{ {{b_{jj'}}} \right\}_{j,j' = 1,2, \ldots ,N}}$

其中

$ \begin{array}{l} \left( {{{{M}}^ * }{e_j}} \right)\left( x \right) = {\left\langle {\left( {{{{M}}^ * }{e_j}} \right)\left( \cdot \right),{K_h}\left( { \cdot ,x} \right)} \right\rangle _{{H_{{K_h}}}}} = \\ {\left\langle {{e_j},{{M}}{K_h}\left( { \cdot ,x} \right)} \right\rangle _{{{\mathbb R}^N}}} = \\ \displaystyle \frac{1}{{2 \text{π} }}\int_{-\text{π}/h}^{\text{π}/h} {\frac{{\left( { - {\eta ^2} - {\rm{i}}\alpha \eta + \beta } \right){{\rm{e}}^{ - {\rm{i}}\eta \left( {x - {x_j}} \right)}}}}{{\lambda {{\left( {1 + {\eta ^2}} \right)}^2} + {{\left| { - {\eta ^2} + {\rm{i}}\alpha \eta + \beta } \right|}^2}}}} {\rm{d}}\eta \end{array} $

元素 ${b_{jj'}}$ 的表达式为

$ {b_{jj'}} = \frac{1}{{2 \text{π} }}\int_{-\text{π}/h}^{\text{π}/h} {\frac{{{{\left| { - {\eta ^2} - {\rm{i}}\alpha \eta + \beta } \right|}^2}{{\rm{e}}^{ - {\rm{i}}\eta \left( {{x_j} - {x_{j'}}} \right)}}}}{{{{\left( {\lambda {{\left( {1 + {\eta ^2}} \right)}^2} + {{\left| { - {\eta ^2} + {\rm{i}}\alpha \eta + \beta } \right|}^2}} \right)}^2}}}} {\rm{d}}\eta $

笔者所考虑的逆问题表述为对于任意 $N$ 个给定值 $\left\{ {F\left( {{x_j}} \right)} \right\}_{j = 1}^N$ ，参数 $C > 0$ ，想要求解

$ \mathop {\min}\limits_{g \in {H_{{K_h}}}} \left\{ {\frac{1}{2}\left\| g \right\|_{{H_{{K_h}}}}^2 + \frac{C}{N}\sum\limits_{j = 1}^N {{{\left| {{{\left\langle {g,{{{M}}^ * }{e_j}} \right\rangle }_{{H_{{K_h}}}}} - F\left( {{x_j}} \right)} \right|}_\varepsilon }} } \right\} $

(8)

问题 $(8)$ 是一个二次规划问题，类似于第2节的解法，通过求解其对偶问题，可以得到问题 $(8)$ 的解

${\overset{\frown} g} \left( x \right) = \sum\limits_{j \in J} {\left( {\theta _j^ * - {\theta _j}} \right)} \left( {{{{M}}^ * }{e_j}} \right)\left( x \right)$

其中 $\theta _j^ * \geqslant 0$ 和 ${\theta _j} \geqslant 0$ 是拉格朗日乘子。

3.2 Sobolev空间 ${H_{{K_s}}}$ 中的解法

在光滑 $s\left( {s \geqslant 1} \right)$ 的Sobolev空间 ${H_{{K_s}}}$ 中，有

$\begin{array}{l} \left( {{{{M}}^ * }{e_j}} \right)\left( x \right) = {\left\langle {\left( {{{{M}}^ * }{e_j}} \right)\left( \cdot \right),{K_s}\left( { \cdot ,x} \right)} \right\rangle _{{H_{{K_s}}}}} = \\ {\left\langle {{e_j},{{M}}{K_s}\left( { \cdot ,x} \right)} \right\rangle _{{{\mathbb R}^N}}} = \\ \displaystyle \frac{1}{{2 \text{π} }}\int_{\mathbb R} {\frac{{\left( { - {\eta ^2} - {\rm{i}}\alpha \eta + \beta } \right){{\rm{e}}^{ - {\rm{i}}\eta \left( {x - {x_j}} \right)}}}}{{{{\left( {1 + {\eta ^2}} \right)}^s}\left( {\lambda {{\left( {1 + {\eta ^2}} \right)}^2} + {{\left| { - {\eta ^2} + {\rm{i}}\alpha \eta + \beta } \right|}^2}} \right)}}} {\rm{d}}\eta \\ \end{array} $

则线性映射 ${{M}}$ 的像空间是一个再生核Hilbert空间，其再生核为^[14]

$\begin{array}{l} {b_{jj'}} = \\ \displaystyle \frac{1}{{2 \text{π} }}\int_{\mathbb R} {\frac{{{{\left| { - {\eta ^2} - {\rm{i}}\alpha \eta + \beta } \right|}^2}{{\rm{e}}^{{\rm{i}}\eta \left( {{x_j} - {x_{j'}}} \right)}}}}{{{{\left( {1 + {\eta ^2}} \right)}^s}{{\left( {\lambda {{\left( {1 + {\eta ^2}} \right)}^2} + {{\left| { - {\eta ^2} + {\rm{i}}\alpha \eta + \beta } \right|}^2}} \right)}^2}}}} {\rm{d}}\eta \end{array} $

对于任意 $N$ 个给定值 $\left\{ {F\left( {{x_j}} \right)} \right\}_{j = 1}^N$ ，参数 $C > 0$ ，笔者想要求解

$\mathop {\min}\limits_{g \in {H_{{K_s}}}} \left\{ {\frac{1}{2}\left\| g \right\|_{{H_{{K_s}}}}^2 + \frac{C}{N}\sum\limits_{j = 1}^N {{{\left| {{{\left\langle {g,{{{M}}^ * }{e_j}} \right\rangle }_{{H_{{K_s}}}}} - F\left( {{x_j}} \right)} \right|}_\varepsilon }} } \right\}$

类似于3.1节，可以得到上述问题的解。

4 误差分析

在这一节中，笔者将给出逼近问题的误差界。为此，首先需要找到解的稳定性和一致性估计，然后得到有界逼近误差。

假设 $\varOmega $ 是一个具有满足内锥条件的Lipschitz边界的有界域，且Hilbert空间 ${H_{{K_\tau }}}\left( \varOmega \right) \subset {H_{{K_s}}}$ 。对于一些离散点 $\left\{ {{x_j}} \right\}_{j = 1}^N \subset \varOmega $ 和任意给定的 $N$ 个值 $\left\{ {g\left( {{x_j}} \right)} \right\}_{j = 1}^N$ ， $\tilde F$ 是问题 $(7)$ 的解。基于文献[20]中的工作，对于解 $\tilde F$ ，笔者有以下稳定性和一致性估计。

定理2 　根据 $LF\left( {{x_j}} \right) = g\left( {{x_j}} \right)$ ，笔者发现，对于每个 $\left\{ {{x_j}} \right\}_{j = 1}^N$ 和每个固定的 $\varepsilon \in {{\mathbb R}^ + }$ ，问题 $(7)$ 的解 $\tilde F$ 满足

${\left\| {\tilde F} \right\|_{{H_{{K_\tau }}}\left( \Omega \right)}} \leqslant {\left\| F \right\|_{{H_{{K_\tau }}}\left( \Omega \right)}}$

和

$\mathop {\max }\limits_{1 \leqslant j \leqslant N} \left| {L\tilde F\left( {{x_j}} \right) - g\left( {{x_j}} \right)} \right| \leqslant \frac{N}{{2C}}\left\| F \right\|_{{H_{{K_\tau }}}\left( \Omega \right)}^2 + \varepsilon$

定理3 　假设 $\varOmega \subset {{\mathbb R}^N}$ 是一个具有满足内锥条件的Lipschitz边界的有界域。令 $\tau $ 为一个正实数， $\displaystyle \left\lfloor {\tau - \frac{1}{2}} \right\rfloor > $ $\displaystyle \frac{N}{2}$ ， $1 \leqslant q \leqslant \infty $ 。假设 $F \in {H_{{K_\tau }}}\left( \varOmega \right)$ 具有 $LF\left( {{x_j}} \right) = g\left( {{x_j}} \right)$ 。设 $\tilde F$ 是问题 $(7)$ 的解，则存在一个常数 $\tilde C > 0$ ，它的取值与τ、 $N$ 和 $\varOmega $ 有关，但与 $\varepsilon $ 、 $F$ 或 $\left\{ {{x_j}} \right\}_{j = 1}^N$ 无关，从而得到逼近误差界为

$\begin{array}{l} {\left\| {F - \tilde F} \right\|_{{L_q}\left( \varOmega \right)}} \leqslant \\ \displaystyle \tilde C\left( {2{h^{\tau - N\left( {\frac{1}{2} - \frac{1}{q}} \right)}} + {{\left\| F \right\|}_{{H_{{K_\tau }}}\left( \varOmega \right)}} + \frac{N}{{2C}}\left\| F \right\|_{{H_{{K_\tau }}}\left( \varOmega \right)}^2 + \varepsilon } \right) \end{array} $

对于所有离散集 $\left\{ {{x_j}} \right\}_{j = 1}^N \subset \varOmega $ 具有填充距离 $h= $ ${h_{\left\{ {{x_j}} \right\}_{j = 1}^N,\varOmega }} \leqslant {C_\varOmega }{\tau ^{ - 2}}$ ，其中

$h = {h_{\left\{ {{x_j}} \right\}_{j = 1}^N,\varOmega }} = \mathop {\sup }\limits_{x \in \varOmega } \mathop {\min }\limits_{{x_j} \in \left\{ {{x_j}} \right\}_{j = 1}^N} {\left\| {x - {x_j}} \right\|_2}$

5 数值实验

为了验证所提出算法的有效性，笔者给出了两个实验。这两个实验都可以在几秒钟内完成。

在实验1中，考虑离散微分方程

$F''\left( x \right) + F'(x) + F(x) = g(x)$

函数 $F$ 和 $g$ 未知，给定区间 $[0,2\text{π}]$ 上函数 $g$ 的300个离散值，要求解出函数 $F$ 。分别在数据没有噪声，数据被不同信噪比的高斯白噪声干扰(SNR=20，SNR=30)3种情况下求解 $F$ 。其真实解为 $F\left( x \right) = \sin {x^2}$ 。近似结果如图1所示，用蓝色实线表示函数 $F$ 的精确解，用红色虚线表示数值结果。相对误差分别是 $1.136\;7 \times {10^{ - 6}}$ ， $4.591\;2 \times {10^{ - 5}}$ ， $9.270\;3 \times {10^{ - 6}}$ 。支持向量的数目分别为15，23，21。可以看出，在噪声干扰的情况下，笔者的方法仍然是有效的。

在实验2中，笔者考虑逆问题。对于离散微分方程

$F''(x) + F'(x) + 4F(x) = g(x)$

函数 $F$ 和 $g$ 未知，在区间 $[0,2{\text{π}} ]$ 上给定 $F\left( x \right)$ 的200个点的值，要求解函数 $g$ 。分别在数据没有噪声，数据被高斯白噪声干扰(SNR=20)两种情况下求解函数 $g$ 。其真实解为 $g\left( x \right) = 3\sin x + \cos x - 2\sin 2x$ 。结果如图2所示，用蓝色实线表示真实解，用红色虚线表示数值结果。相对误差分别是 $4.839\;2 \times {10^{ - 5}}$ ， $9.362\;1 \times $ ${10^{ - 5}}$ 。支持向量的数目分别为9和11。实验结果表明，笔者的方法是有效且对噪声是鲁棒的。

6 结论

本文通过将支持向量机、Tikhonov正则化和再生核理论相结合，提出了一种求解离散线性微分方程及其逆问题的方法。通过将离散微分方程问题及其逆问题转换为对应的目标优化问题，然后在支持向量机的框架下进行求解，得到具有解析表达式的近似解，所得解具有稀疏性，对噪声鲁棒，结构固定的特点，便于后续应用。离散微分方程及其逆问题均有相同结构的解。实验结果证明了该方法的有效性。