Andy123t
0%

SMM01 微分矩阵

发表于 更新于 分类于

本文是书籍 Spectral Method in MATLAB 第一章 Differentiation matrices 微分矩阵 笔记。

从一个基本问题出发。 给定一组网格节点 $x_i$ 和相应的函数值 $u(x_i)$，如何使这些数据逼近 $u$ 的导数？可能立即想到的方法是某种有限差分公式。 我们将从有限差分方法开始去学习谱方法。

二阶差分

考虑等距剖分 $\{ x_{1},\ldots,x_{N}\}$, 对于每一个 $j$ 有 $x_{j+1}-x_{j}=h$，这组点对应的函数值 $\{u_{1}, \ldots, u_{N}\}$：

SMM0101.png

令 $w_i$ 表示 $u’(x_i)$ 的逼近，$u’(x_i)$ 表示 $u$ 的导数在 $x_j$ 处的值。 标准二阶有限差分格式

可以通过 $u(x_{j+1})$ 与 $u(x_{j+1})$ 的泰勒展开得到。
为了简单起见，假设问题是周期性的，取 $u_0=u_N$ 和 $u_1=u_{N+1}$。然后我们可以把离散微分过程表示为矩阵向量乘积

稀疏矩阵中的省略项为零。观察这个矩阵可知是 Toeplitz 矩阵，沿着对角线有相同的数，比如元素 $a_{ij}$ 的值只与 $i-j$ 有关。 实际上，它也是循环的，这意味着 $a_{ij}$ 仅依赖于 $(i-j)(\mathrm{mod} N)$。 这个矩阵对角线 “环绕” 的。

另一种推导以上两个方程的方法是通过以下的局部插值和微分过程:

对于 $j=1,2, \ldots, N$，

  • 令 $p_j$ 是不大于 2 次的多项式，$p_{j}(x_{j-1})=u_{j-1}$, $p_{j}(x_{j})=u_{j}$, 与 $p_{j}(x_{j+1})=u_{j+1}$.
  • 设 $w_{j}=p_{j}^{\prime}(x_{j})$.

可以很容易地看出，对于固定的插值函数 $p_j$

其中

  • $a_{-1}(x)=(x-x_{j})(x-x_{j+1}) / 2 h^{2}$，
  • $a_{0}(x)=-(x-x_{j-1})(x-x_{j+1}) / h^{2}$，
  • $a_{1}(x)=(x-x_{j-1})(x-x_{j}) / 2 h^{2}$.

在 $x=x_j$ 时求微分并求值可以得到上式 $ w_{j}=\frac{u_{j+1}-u_{j-1}}{2 h}$。

四阶差分

这种局部插值的推导表明能够推广到更高阶。下面是四阶差分格式推导：

对于 $j=1,2, \ldots, N$，

  • 令 $p_j$ 是不大于 4 次的多项式，$p_{j}(x_{j \pm 2})=u_{j \pm 2}$, $p_{j}(x_{j \pm 1})=u_{j \pm 1}$, 与 $p_{j}(x_{j})=u_{j}$.
  • 设 $w_{j}=p_{j}^{\prime}(x_{j})$.

假设数据是周期的，可知方程组等价于矩阵向量乘积

这是一个五对角而不是三对角的循环矩阵。

以上两个矩阵的矩阵是微分矩阵的例子，它们的精度分别为2阶和4阶。对一个足够光滑的函数 $u$ 进行采样获得的 $u_j$，相应的离散逼近的收敛速率分别为 $O(h^{2})$ 与 $O(h^{4})$，通过泰勒展开可以验证这一点。

Program 1

第一个 MATLAB 程序，Program 1，演示了四阶差分方法的例子，选取 $u(x)=e^{\sin (x)}$ 在 $[-\pi, \pi]$ 上的周期数据。

SMM0102.png

该程序对不同 $N$ 时差分逼近 $ w_j $ 与精确导数 $ e^{\sin (x_{j})} \cos (x_{j}) $ 进行了比较。 因为使用了 MATLAB 稀疏矩阵 sparse ，这个代码在电脑上运行几分之一秒就运行出来了，即使它矩阵维数大到 4096。 四阶精度是明显的。 这是我们第一个也是最后一个没有使用谱方法的例子。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
% p1.m - convergence of fourth-order finite differences

% For various N, set up grid in [-pi,pi] and function u(x):
  Nvec = 2.^(3:12);
  clf, subplot('position',[.1 .4 .8 .5])
  for N = Nvec
    h = 2*pi/N; x = -pi + (1:N)'*h;
    u = exp(sin(x)); uprime = cos(x).*u;

    % Construct sparse fourth-order differentiation matrix:
    e = ones(N,1);
    D =   sparse(1:N,[2:N 1],2*e/3,N,N)...
        - sparse(1:N,[3:N 1 2],e/12,N,N);
    D = (D-D')/h;

    % Plot max(abs(D*u-uprime)):
    error = norm(D*u-uprime,inf);
    loglog(N,error,'.','markersize',15), hold on
  end
  grid on, xlabel N, ylabel error
  title('Convergence of fourth-order finite differences')
  semilogy(Nvec,Nvec.^(-4),'--') 
  text(105,5e-8,'N^{-4}','fontsize',18)

SMM0103.png

我们研究了二阶和四阶有限差分。很明显，考虑六阶、八阶和更高阶的方法将导致带宽增加的循环矩阵。谱方法背后的思想是，将这一过程带到极限，并使用无穷阶和无限带宽的微分公式——比如一个稠密矩阵。在下一章中，我们将证明在这个极限下，对于一个无限等距网格，我们得到以下无限矩阵

这是一个反对称（$D=-D$）两边无限的 Toeplitz 矩阵，也称为 Laurent 算子。除主对角线上的项外，它的所有项都不是零。

谱配置方法

实际操作中一个矩阵不可能是无限的。对于有限网格，下面是谱配置方法

  • 令 $p$ 是一个单独的函数（与 $j$ 无关），使得对所有的 $j$ 成立 $p(x_{j})=u_{j}$.
  • 设 $w_{j}=p^{\prime}(x_{j})$.

我们可以自由选择 $p$ 来解决问题。 对于一个周期性的区域，自然的选择是在一个等步长的网格上的三角多项式，由此产生的 “Fourier” 方法将在第4章和后面章节中提到。 对于非周期区域，不规则网格上的代数多项式是正确的选择，我们将从第5章和第6章开始介绍这种类型的 “Chebyshev” 方法。

对于有限的 $N$，为了简单起见直接取 $N$，这是在第3章中为周期的规则网格导出 $N \times N$ 的稠密矩阵

Program 2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
% p2.m - convergence of periodic spectral method (compare p1.m)

% For various N (even), set up grid as before:
  clf, subplot('position',[.1 .4 .8 .5])
  for N = 2:2:100;
    h = 2*pi/N;
    x = -pi + (1:N)'*h;
    u = exp(sin(x)); uprime = cos(x).*u;

    % Construct spectral differentiation matrix:
    column = [0 .5*(-1).^(1:N-1).*cot((1:N-1)*h/2)];
    D = toeplitz(column,column([1 N:-1:2]));

    % Plot max(abs(D*u-uprime)):
    error = norm(D*u-uprime,inf);
    loglog(N,error,'.','markersize',15), hold on
  end
  grid on, xlabel N, ylabel error
  title('Convergence of spectral differentiation')

SMM0104.png
注：谱方法具有 “谱精度”，直到舍入误差超过 $10^{-14}$。这个矩阵是稠密的，但是 $N$ 的值比程序1小得多。

观察余切函数，可以发现这个矩阵确实和 Toeplitz 矩阵一样是循环的。

程序2与程序1相同，只是谱方法替换了差分方法，结果有什么不同！输出2中的误差下降得非常快，直到达到如此高的精度，以至于计算机上的舍入误差阻止了进一步的改进。这种显著的行为称为谱精度。我们将在第4章中给出这个短语的更精确的描述，但是现在，需要注意的是，它与有限差分法和有限元法的收敛速度有多大的不同。随着 $N$ 的增加，有限差分法或有限元法中的误差收敛速率 $\color{red}{O(N^{-m})}$，常数 $\color{red}{m}$ 取决于逼近阶数和解的光滑性。对于谱方法。如果解是无穷可微的，对于每一个 $\color{red}{m}$ 的收敛速度为 $\color{red}{O(N^{-m})}$，如果解是适当解析的，收敛速度甚至比 $\color{red}{O(c^N)}$（$\color{red}{0 < c < 1}$）的速度更快。

我们所描述的矩阵是循环的。循环矩阵的作用是卷积，我们将在第3章中看到，卷积可以通过离散傅里叶变换来计算。从历史上看，正是1965年快速傅里叶变换（FFT）的发现，导致了人们对谱方法的兴趣在1970年代激增。第8章 FFT 不仅适用于等距网格上的三角多项式，也适用于Chebyshev网格上的代数多项式。然而，即使没有 FFT 的情况下谱方法也是很强大的。在许多应用中，使用矩阵计算的结果是非常令人满意的。

舍入误差

关于舍入误差，标准IEEE双精度格式 (Double-precision floating-point format) $\color{red}{\epsilon_{\text{machine}}=2 ^ { - 53 } \approx 1.11 \times 10 ^ {-16}}$, 这意味着每次加法, 乘法, 除法和减法都会产生完全正确的结果乘以某个因子 $1 + \delta$ with $|\delta| \leq \epsilon_{\text{machine}}$.

本章小结

谱配置法的基本原理是在网格上给定离散数据，对数据进行全局插值，然后计算插值在网格上的导数。 对于周期性问题，我们通常在等距点上使用三角插值，对于非周期性问题，我们通常在不等距点上使用多项式插值。