数值微分问题就是已知近似函数或者在一些离散点的值, 求它在某点的导数或高阶导数的近似值.这个问题产生于数学模型以及实际问题中, 例如, 积分方程问题[1]、图像处理中的边界识别问题[2]、数学物理方程中的一些反问题[3].然而, 数值微分面临的问题是观测数据的微小误差也可能会造成数值结果的巨大误差[4], 因此, 如何克服这种不稳定也就成为了数值微分方法研究中的一个重要的课题.
处理数值微分问题已经有了大量的方法:磨光法[5]、差分法[6-8]、基于Tikhonov正则化理论的方法[9, 10]以及积分法[11-13].其中利用积分算子法可以给出一致的误差估计, 而且当函数的光滑性加强时, 可以构造类似的积分算子使得误差精度提高[14].许多学者已经利用差分法和基于Tikhonov正则化理论的方法对二阶数值微分问题作出了相应的研究, 而积分算子法对于高阶甚至任意阶数值微分问题还鲜见探讨.本文基于正交多项式的相关知识, 将引入新的积分算子, 控制近似计算过程中的不稳定扰动, 以达到可以稳定逼近任意阶导数的目的, 同时提高了新算子在近似逼近时误差的精度.
在这一部分内容, 我们基于正交多项式的背景下, 引入了一个新的积分算子, 可以对近似已知函数有良好的逼近程度.对于不适定问题而言, 稳定解的收敛率分析是一个十分重要的问题, 我们所提出的积分算子, 一个重要的特征是提高了收敛率, 在本节内容中, 将指出, 在一定的条件下, 精度可以达到$O(\delta^{\frac{q+1}{n+q+1}})$, $q\in N$.$\delta$是近似数据的误差界.
记$I=[a, b]$, 假设$f^{\delta}=f+\tilde{\omega}$是一个定义在$I$上近似于$f$的已知函数, 其中$f\in C^{n+q+1}(I)$且$\tilde{\omega}$是一噪声水平为$\delta$的有界可积函数,
引入积分算子:
其中, $h>0$是参数,
并把$(D_{h}^{(n)}f^{\delta})(x)$作为函数$f^{n}(x)$的近似, 即$n$阶导数的近似逼近.
引理2.1 令
则$\{p_{q}(t)\}_{q=0}^{\infty}$就是在区间[0, 1]上关于权函数$(1-t)^{\frac{1}{2}}\cdot t^{n}$的正交多项式.
证 要证$\{p_{q}(t)\}_{q=0}^{\infty}$就是在区间[0, 1]上关于权函数$(1-t)^{\frac{1}{2}}\cdot t^{n}$的正交多项式, 即证
不妨设$l<m$, 我们先证明
因为根据行列式的性质
由式(2.5) 可知, $p_{m}(t)$与$1, t, \cdots, t^{m-1}$在区间[0, 1]上关于权函数$(1-t)^{\frac{1}{2}}t^{n}$的正交多项式.
不妨记
那么$a_{i}$就是行列式(2.3) 第$q$列的元的代数余子式.
根据引理2.1, 我们可以得到如下推论, 此推论在后续内容中有着重要的应用.
推论2.1
证 下面分对于$m$的奇偶性分两种进行证明:
(1) 若$m$为奇数, 结论显然成立;
(2) 若$m$为偶数, 不妨令$m = 2k$, 因为$\sum\limits_{i = 0}^q {a_{i}{t^i}}$是在区间[0, 1]上关于权函数${(1 - t)^{\frac{1}{2}}}{t^n}$的正交多项式, 此即
根据对称性, 并作变量代换$u = 1 - {t^2}$可得
由式(2.7) 可知
以下定理将说明, $f(x)$在$x$处不存在普通意义下的$n$阶导数, 但极限$\mathop {\lim }\limits_{h \to {0^ + }} {D_h}^{\left( n \right)}f\left( x \right)$仍有可能存在.在假设函数${f^{\left( n \right)}}$左、右导数都存在但不一定相等的情况下, 我们证明了$\left( {D_{h}^{\left( n \right)}f} \right)\left( x \right)$确为一种广义的导数, 即
定理2.1 假设$f \in {C^{n - 1}}\left( I \right)$, 如果${f_ - }^{\left( n \right)}\left( x \right)$和${f_ + }^{\left( n \right)}\left( x \right)$都存在, $x \in {I_h}$, 那么
其中${I_h} = [a+h, b-h]$.
证 根据带有佩亚诺余项的泰勒公式可知, 对于任给的$\varepsilon$, 存在$\delta > 0$, 有
其中${f_{n - 1}}(x + ht)$是$f\left( {x + ht} \right)$的前$n - 1$项泰勒展开式.
而又因为${\left[{\sum\limits_{i = 0}^q {a_{i}{{(1-{t^2})}^{n + i}}} } \right]^{\left( n \right)}}{t^n}$是关于$t$的奇函数, 则
直接计算可得
由已知, 对于任意的$x \in {I_h}$, 均有
又
那么对于任意的$\varepsilon ' > 0$, 只需要取
存在$\delta$, 有$0 < h < \delta $,
此即
式(2.8) 表明, 我们可以利用$\left( {{D_h}^{\left( n \right)}f} \right)\left( x \right)$来逼近导数${f^{\left( n \right)}}\left( x \right)$, 此时不需要限制${f^{\left( n \right)}}$具有通常意义下的可导性, 这就是之所以称$\left( {{D_h}^{\left( n \right)}f} \right)\left( x \right)$是广义导数的原因.在实际应用中, 往往知道的是$f\left( x \right)$的近似函数${f^\delta }\left( x \right)$或者离散点上的测量数据${f^\delta }\left( {{x_i}} \right)$.由于噪声的作用, ${f^\delta }\left( x \right)$可能不具有通常意义下的导数, 因此广义导数$\left( {{D_h}^{\left( n \right)}f} \right)\left( x \right)$提供了一种从带有噪声的测量数据近似求出未知函数导数的有效途径.
我们将在定理2.2及推论2.2中证明所引入的${D_h}^{\left( n \right)}f\left( x \right)$是${f^{\left( n \right)}}\left( x \right)$的一个任意阶稳定的近似逼近.为考虑方便, 我们暂时忽略噪声$\tilde \omega $的影响, 在后面的部分中, 我们会加入对噪声$\tilde \omega $的考虑.
定理2.2 假设$f \in {C^{n + q + 1}}\left( I \right)$且$q \in N$, 则${D_h}^{\left( n \right)}f\left( x \right) = {f^{\left( n \right)}}\left( x \right) + O({h^{q + 1}}).$
证 根据泰勒公式, 对于任意的$x \in {I_h}$, 存在$\xi \in \left[{x-h, x + h} \right]$, 有
事实上, 运用分部积分法:
将(2.6), (2.9) 式代入${D_h}^{\left( n \right)}f\left( x \right)$的表达式并化简, 有
推论2.2 假设$f \in {C^{n + q + 1}}\left( I \right)$且$q \in N$, $\left| {{f^{\left( {n + q + 1} \right)}}\left( x \right)} \right| \le {M_{n + q + 1}}$, 那么有
其中
证 由推论2.1及(2.9) 式, 有
由于存在${M_{n + q + 1}} > 0$, 使得对于任意的$x \in I$, 有$\left| {{f^{\left( {n + q + 1} \right)}}\left( x \right)} \right| \le {M_{n + q + 1}}$, 则
如果记
那么上式即为(2.10) 式.
下面, 我们将考虑加入噪声$\tilde \omega $的情况.如果我们要考虑带有噪声的函数${f^\delta }\left( x \right)$, 那么在估计${f^{\left( n \right)}}\left( x \right)$时, $f(x + h \cdot )$就要用${f^\delta }(x + h \cdot )$替代.以下我们将证明, 若适当选取$h = h(\delta )$, ${D_h}^{\left( n \right)}{f^\delta }\left( x \right)$仍然是${f^{\left( n \right)}}\left( x \right)$的一个稳定的近似逼近.
定理2.3 令${f^\delta }$是一个噪声函数, 假设$f \in {C^{n + q + 1}}\left( I \right)$且$\tilde \omega$是一噪声水平为$\delta $的有界可积函数, 那么
证 结合(2.1) 式,
根据三角不等式的性质, 结合(2.10) 式,
推论2.3 假设$f \in {C^{n + q + 1}}\left( I \right)$, 如果我们取$h = {\left[{\frac{{n{C_2}}}{{\left( {q + 1} \right){C_1}}}\delta } \right]^{\frac{1}{{n + q + 1}}}}$, 那么
证 记$\varphi \left( h \right) = {C_1}{h^{q + 1}} + {C_2}\frac{\delta }{{{h^n}}}$, 令$\varphi '\left( h \right) = 0$, 得${h^*} = {\left[{\frac{{n{C_2}}}{{\left( {q + 1} \right){C_1}}}\delta } \right]^{\frac{1}{{n + q + 1}}}}$, 且经过验证, 此点为最小值点,
即
从式(2.12) 可以看出, 对于新引入的积分算法, 在一定的条件下, 近似导数的收敛率达到了$O({\delta ^{\frac{{q + 1}}{{n + q + 1}}}})$, 精度较Lanczos广义导数法有了很大程度的提高.
为了验证上述估计的有效性与准确性, 我们给出了两个数值实验.这些实验通过MatlabR2012a测试执行.令${f^\delta }\left( {{x_i}} \right) = f\left( {{x_i}} \right) + c\tilde \omega \left( {{x_i}} \right)$是一系列噪声数据且采样周期为$T = {10^{ - 4}}$, $c\tilde \omega \left( {{x_i}} \right)$是模拟一个零均值的高斯白噪声序列.通过$3\sigma$原则, 我们假定噪声水平$c\tilde \omega \left( {{x_i}} \right) = 3c$.选取积分区间为${x_i} \in I = \left[{-2, 2} \right]$, 积分计算时采用复合梯形求积公式.
例1 ${f_1}\left( x \right) = {e^{ - {x^2}}} \cdot \cos 3\pi x$, 取噪声水平$\delta=0.15$, $\delta=0.015$时, 实验的相对误差
图 1表示噪声水平是$\delta = 0.15$时数值实验的结果, 实线表示${f_1}^{(n)}\left( x \right)$在$n = 1, 2, 3, 4$处的导数值, 点线表示导数的估计值${D_h}^{(n)}{f_1}\left( {{x_i}} \right)$.
例2 考虑分段函数
取噪声水平$\delta = 0.15$. ${C^2} \in I\left[{-2, 2} \right]$, ${f_2}^{(2)}\left( x \right) \equiv \left| x \right|$, ${f_2}^{(3)}\left( x \right)在$x = 0不存在. 表 2给出了在噪声水平$\delta = 0.15$时实验的相对误差
图 2表示噪声水平是$\delta = 0.15$, $\delta = 0.015$时数值实验的结果, 实线表示${f_2}^{(n)}\left( x \right)$在$n = 1, 2, 3$处的导数值, 点线表示导数的估计值${D_h}^{(n)}{f_2}\left( {{x_i}} \right)$
本文所引进的积分算子, 较之Lanczos法, 将精度提高到了$n$阶, 这样逼近程度大大增加了.上述2个数值实验表明, 在加入随机噪声的情况下, 对于给定区间的点的导数, 引入的积分算子逼近效果是比较理想的. 表 1, 表 2的相对误差均在可接受的范围之内, 表 2的相对误差仅仅只有0.01左右, 其中在不可导点$x=0$处, 并没有出现剧烈的震荡.然而, 由于拉格朗日展开式的局限性, 使得我们例子所假设的积分区间${x_i} \in I = \left[{-2, 2} \right]$, 其实际有效区间是${x_i} \in I = \left[{-2+h, 2-h} \right]$, 这就导致了端点处出现轻微震荡效果.
在实际应用中, 如果想求得函数在区间的近似导数, 数值实验表明, 利用本文所提方法去逼近此导数则是切实可行、有效的, 其计算结果也是稳定的.虽然此积分算子方法本身存在一定的局限性, 使得其无法得出端点附近的导数, 但是通过数值实验结果显示, 我们所得出的近似导数逼近效果仍然是十分令人满意的.对于端点附近的导数的计算问题, 我们可以考虑对近似已知函数作适当的延拓, 考虑构造新的正则化方法解决.