数值分析

2020-11-20 本文已影响0人 Zeabin

1. 数学基础

泰勒公式

If $f \in C^n[a, b]$ and if $f^{(n+1)}$ exists on the open interval $(a, b)$ , then for any points $c$ and $x$ in the closed interval $[a, b]$ ,
$f(x) = \sum^n_{k=0} \frac{1}{k!}f^{(k)}(c)(x-c)^k + E_n(x)$
for some point $\xi$ between $c$ and $x$
$E_n(x) = \frac{1}{(n+1)!}f^{(k+1)}(\xi)(x-c)^{n+1}$

收敛阶数 order of convergence

描述一个序列的收敛速度
如果 $|x_{n+1} - x^*| \le C |x_n - x^*|^{\alpha}, (n \ge N)$ ，那么序列 $[x_n]$ 的收敛阶数为 $\alpha$

使用迭代公式 $x_{k+1} = \frac{x_k(x_k^2 + 3a)}{3x_k^2+a}$ 来计算 $\sqrt{a}$ ，求 $[x_k]的收敛阶数$
$\begin{align*} |x_{k+1} - \sqrt{a}| & = |\frac{x_k(x_k^2 + 3a)}{3x_k^2+a} - \sqrt{a} | \\ & = |\frac{x_k^3 + 3a x^k - 3x_k\sqrt{a} + a\sqrt{a} }{3x_k^2+a}| \\ & = |\frac{(x_k - \sqrt{a})^3}{3x_k^2+a}| \\ & \approx \frac{1}{4a}|(x_k - \sqrt{a})|^3 , (x_k \to \sqrt{a}) \end{align*}$
$\therefore$ 收敛阶数为3

定理：对于迭代公式 $x_{k+1} = \varphi(x_k)$ ，如果 $\varphi^{(p)}(x)$ 在 $x^*$ 附近收敛，并且
$\varphi'(x^*) = \varphi''(x^*) = \cdots = \varphi^{(p-1)}(x^*) = 0, \varphi^p(x^*) \ne 0$ 那么迭代公式的收敛阶数为 $p$

2. 计算机算法

计算机使用二进制表示数字，由于有限的数字位数，无法精确地表示每一个数字，可能需要进行舍入

绝对误差： $|x - x^*|$
相对误差： $\left| \frac{x - x^*}{x^*} \right|$

两个非常接近的数字相减会带来较大地误差，通过分式上下同乘两数相加来避免

$9 - \sqrt{80} = \frac{1}{9 + \sqrt{80}}$

一个数值计算过程是不稳定的，如果在这个过程中一个阶段的误差被下一个阶段放大

3. 解线性方程组

解 $Ax = b$ ，先将 $A$ 化成两个简单的上三角、下三角矩阵 $A = LU$
解 $Lz = b$ ，得到 $z$
解 $Ux = z$ ，得到原方程的解 $x$

Doolittle分解

令 $L_{ii} = 1$ ，依次更新 $U$ 的行， $L$ 的列

LDU分解

对Doolittle分解的 $U$ 提取对角线元素，分解成 $U = DU'$
$A = LDU$

Crout分解

令 $U_{ii} = 1$ ，依次更新 $L$ 的列， $U$ 的行

Cholesky分解

$A = L L^T$
$A$ 需是实对称正定矩阵

高斯消元法 with Scaled Row Pivoting

$PA = LU$
$s_i = \max_{ 1 \le i \le n } |a_{ij}| = \max \{ |a_{i1}|, |a_{i2}|, \cdots, |a_{in}| \}$ 取每一行绝对值最大的数
取 $|a_{ik}|/s_i$ 最大的行 $i$ 作为第 $k$ 次消元的pivot row， $a_{ik}$ 为第 $k-1$ 次更新后的值，前面已选为pivot row的行不参与

范数

矩阵 $A$ 的范数
$\infty$ 范数： $\max_{1 \le i \le n } \sum_{j=1}^{n}|a_{ij}|$ 绝对值行求和取最大
1范数： $\max_{1 \le j \le n } \sum_{i=1}^{n}|a_{ij}|$ 绝对值列求和取最大
2范数： $\sqrt{\rho(A^T A)}$ $A^T A$ 的谱半径（最大特征值）开平方

用迭代法解方程

$Ax = b$
$Qx = (Q-A)x+b$
The equation suggests an iterative process
$Qx^{(k)} = (Q-A)x^{(k-1)}+b (k \ge 1)$
$x^{(k)} = (I-Q^{-1}A)x^{(k-1)}+Q^{-1}b$

$x = (I-Q^{-1}A)x+Q^{-1}b$
$x^{(k)} - x = (I-Q^{-1}A)(x^{(k-1)} - x)$
$||x^{(k)} - x|| \le ||I-Q^{-1}A|| \ ||x^{(k-1)} - x||$
如果 $||I-Q^{-1}A|| < 1$ ，则迭代法收敛

$A = D-L-U$

Richardson: $Q = I, x^{k} = (I-A) x^{k-1}+ b$

Modified Richardson: $Q = \frac{1}{ \omega } I, x^{k} = (I-\omega A) x^{k-1}+\omega b = x^{k-1} + \omega (b - A x^{k-1})$

Jacobi: $Q = D$

Gauss-Seidel: $Q = D-L$

SOR(successive over relaxation): $Q = \frac{1}{\omega}(D-\omega L)$

Consider a system
$\left( \begin{matrix} 1 & 2 & -2 \\ 1 & 1 & 1 \\ 2 & 2 & 1 \end{matrix} \right) \left( \begin{matrix} x_1 \\ x_2 \\ x_3 \end{matrix}\right) = \left(\begin{matrix} 1 \\ 2 \\ 3 \end{matrix}\right)$
discuss its convergence when using Jacobi and G-S methods.
Solution: decomposing the coefficient matrix such that $A = D - L - U$
Jacobi:
$G_J = D^{-1}(L+U)$
解特征方程 $|\lambda I - G_J|$ 得到特征值 $\lambda_1 = \lambda_2 = \lambda_3 = 0$
$\rho(G_J) = 0 < 1$ , Jacobi方法收敛
Gauss-Seidel:
$G_{G-S} = (D-L)^{-1}U$
解特征方程 $|\lambda I - G_{G-S}|$ 得到特征值 $\lambda_1 = 0, \lambda_{2,3} = 2$
$\rho(G_{G-S}) = 2 > 1$ , Gauss-Seidel方法不收敛

4. 函数逼近

多项式插值

给定一组点 $(x_0, y_0), \cdots, (x_n, y_n)$ ，找到一个插值多项式 $p(x)$ 使得 $p(x_i) = y_i \ (i = 0, \cdots, n)$
误差项 $R_n(x) = f(x) - p_n(x) = \frac{f^{(n+1)}(\xi) }{(n+1)! } \omega_{n+1}(x), \ w_{n+1}(x) = \prod_{i=0}^{n}(x - x_i)$
$| R_n(x) | \le \frac{ M_{n+1 } } {(n+1)! } |\omega_{n+1} (x)|$

拉格朗日插值法

$l_i(x_j) = \begin{cases} 0, & i \ne j \\ 1, & i = j \end{cases}$
$p_n(x) = \sum_{i=0}^{n} l_i(x) y_i$
$l_i(x) = \prod_{ j=0\\ j \ne i }^{ n } \frac{ x - x_j }{ x_i - x_j }$

牛顿插值法

$f[x_0, \cdots, x_k] = \frac{f[x_0, \cdots, x_{k-1}] - f[x_{1}, \cdots, x_k]}{x_0 - x_k}$
$f(x) = N_n(x) + R_n(x)$
$N_n(x) = f(x_0) + (x - x_0)f[x_0, x_1] + (x - x_0)(x - x_1)f[x_0, x_1, x_2] + \cdots$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ + (x - x_0)(x - x_1) \cdots (x-x_{n-1})f[x_0, x_1, \cdots, x_n ]$
$R_n(x) = \omega_{n+1}(x) f[x, x_0, \cdots, x_n]$

等间距： $x_n = x_0 + n h$
迎风差： $\Delta f_k = f_{k+1} - f_k, \ \Delta^m f_k = \Delta^{m-1}f_{k+1} - \Delta^{m-1}f_k$
$f[x_i, \cdots, x_{i+m}] = \frac{1}{m! h^m} \Delta^m f_i$
$N_n(x_0 + t h) = f(x_0) + t \Delta f_0 + \frac{ t(t-1) }{ 2 } \Delta^2 f_0 + \cdots + \frac{ t(t-1) \cdots (t-n+1) }{ n! } \Delta^n f_0$
$R_n(x_0 + th) = \frac{t(t-1) \cdots (t-n)}{(n+1)!} h^{ n+1 } f^{ (n+1) } (\xi)$

厄米特插值法

插值函数在给定点的函数值、导数都需与给定值相等
$a \le x_0 \le \cdots \le x_n \le b, \ y_i = f(x_i), \ m_i = f'(x_i)$
$H(x_i) = y_i, H'(x_i) = m_i$
2n+2个条件，可以确定一个2n+1阶多项式

构造基函数
$\begin{cases} \alpha_j(x_k) = \delta_{jk}, \alpha'_j(x_k) = 0 \\ \beta_j(x_k) = 0, \beta'_j(x_k) = \delta_{jk} \end{cases} \ \ \delta_{jk} = \begin{cases} 0, j \ne k \\ 1, j = k \end{cases} \ (j, k = 0, 1, \cdots, n)$
$H_{2n+1}(x) = \sum_{j = 0}^n [\alpha_j(x) y_j + \beta_j(x) m_j]$

最小二乘逼近

$\left(\begin{matrix} \sum_{i=1}^{m} x_i^0 & \sum_{i=1}^{m} x_i^1 & \cdots & \sum_{i=1}^{m} x_i^n \\ \sum_{i=1}^{m} x_i^1 & \sum_{i=1}^{m} x_i^2 & \cdots & \sum_{i=1}^{m} x_i^{n+1} \\ \vdots & \vdots & \ddots & \vdots \\ \sum_{i=1}^{m} x_i^n & \sum_{i=1}^{m} x_i^{n+1} & \cdots & \sum_{i=1}^{m} x_i^{2n} \\ \end{matrix}\right) \left(\begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{matrix}\right) = \left(\begin{matrix} \sum_{i=1}^{m} y_i x_i^0 \\ \sum_{i=1}^{m} y_i x_i^1 \\ \vdots \\ \sum_{i=1}^{m} y_i x_i^n \\ \end{matrix}\right)$

函数逼近
$\left(\begin{matrix} \int_{a}^{b} x^0 dx & \int_{a}^{b} x^1 dx & \cdots & \int_{a}^{b} x^n dx \\ \int_{a}^{b} x^1 dx & \int_{a}^{b} x^2 dx & \cdots & \int_{a}^{b} x^{n+1} dx \\ \vdots & \vdots & \ddots & \vdots \\ \int_{a}^{b} x^n dx & \int_{a}^{b} x^{n+1} dx & \cdots & \int_{a}^{b} x^{2n} dx \\ \end{matrix}\right) \left(\begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{matrix}\right) = \left(\begin{matrix} \int_{a}^{b} f(x) x^0 dx \\ \int_{a}^{b} f(x) x^1 dx \\ \vdots \\ \int_{a}^{b} f(x) x^n dx \\ \end{matrix}\right)$

正交多项式

如果 $f(x), g(x) \in C[a, b]$
$(f(x), g(x)) = \int_a^b \rho(x)f(x)g(x) dx = 0$ 则称 $f(x), g(x)$ 在权函数 $\rho(x)$ 下正交

勒让德正交多项式

区间为 $[-1, 1], \rho(x) = 1$
$P_0 (x) = 1, \ P_1 (x) = x, \ P_2 (x) = (3x^2 - 1)/2, P_3(x) = (5x^3 - 3x)/2, \cdots$

切比雪夫正交多项式

区间为 $[-1, 1], \rho(x) = \frac{1}{\sqrt{1-x^2 }}$
$T_n(x) = \cos(n \arccos x)$
$T_0(x) = 1, T_1(x) = x, T_2(x) = 2x^2-1, T_3(x) = 4x^3-3x, \cdots$

5. 数值微分与积分

数值微分

$f'(x) = \frac{1}{h}[f(x+h) - f(x)], E = -\frac{h}{2}f''(\xi )$
$f'(x) = \frac{ 1 }{ 2h }[ f(x+h) - f(x-h) ], E = -\frac{ h^2 }{ 6 } f''( \xi )$

Richardson外推法

反复将 $h, \frac{h}{2 }$ 带入泰勒展开式、组合，消去低阶项

数值积分

Newton-Cotes

使用拉格朗日插值多项式近似函数，对多项式积分
$\int_a^b f(x) dx \approx \sum_{i=0}^n A_i f(x_i), \ A_i = \int_a^b l_i(x) dx$

Trapezoid

使用一次函数近似，计算梯形面积
$\int_a^b f(x) dx \approx \frac{b-a}{2}[f(b) + f(a)]$

Composite Trapezoid

将区间 $[a, b]$ n等分， $h = \frac{b-a}{n}, x_i = a+ih$ , 计算n个梯形面积
$\int_a^b f(x) dx \approx \frac{h}{2}[f(a) + 2\sum_{i=1}^{n-1}f(a+ih) + f(b)]$
$R = - \frac{1}{12}(b-a)h^2 f''(\xi)$

Simpson

$\int_a^b f(x) dx \approx \frac{b-a}{6}[f(a) + 4f(\frac{b+a}{2}) + f(b)]$
$R_S = - \frac{1}{ 90 }[\frac{b-a}{2} ]^5 f^{(4)}(\xi)$

Composite Simpson

n等分区间 $[a, b]$ ， $h = \frac{b-a}{n}, x_i = a+ih$
$\int_a^b f(x) dx \approx \frac{h}{3} \sum_{i=1}^{n/2} [f(x_{2i-2}) + 4f(x_{2i-1}) + f(x_{2i}) ]$
$R = - \frac{1}{180 }(b-a) h^4 f^{(4)} (\xi)$

Gaussian Quadrature

寻找最合适的点和系数，使得近似的精度最高
$n+1$ 个系数 $A_i$ ， $n+1$ 个点 $x_i$ ，可以确定精度最高 $2n+1$ 阶多项式

选择的 $n$ 个点为 $n$ 阶正交多项式的零点
选择正交多项式要根据权重函数 $w(x)$ ，区间 $[a, b]$ ，若区间不是 $[-1, 1]$ ，则需要作变量代换化成 $[-1, 1]$

当权重函数 $w(x) = \frac{1}{\sqrt{1-x^2 }}$ ，选择切比雪夫正交多项式
$\int_{-1}^1 \frac{f(x)}{\sqrt{1-x^2}} dx = \frac{\pi}{n} \sum_{k=1}^n f(x_k), \ x_k = \cos \frac{2k-1}{2n} \pi \ (k = 1, 2, \cdots, n)$

6. 非线性方程求根

Bisection Method

$r = \lim_{n \to \infty} c_n, \ c_n = \frac{1}{2}(a_n + b_n)$ , then
$|r - c_n| \le \frac{1} {2^{n+1}}(b_0 - a_0)$

$|e_n| = \frac{1}{2} |e_n|$
收敛阶数： $1$

Newton Method

$r = x_n + h$
$0 = f(r) \approx f(x_n) + hf'(x_n)$
$h = - \frac{ f(x_n) }{ f'(x_n) }$
$x_{n+1} = x_n + h = x_n - \frac{f(x_n)}{f'(x_n)}$

$|e_{n+1}| \approx C|e_n|^2$
收敛阶数： $2$

Secant Method

$x_{n+1} = x_n - \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})} f(x_n)$

使用 $\frac{f(x_n) - f(x_{n-1}) }{x_n - x_{n-1} }$ 近似 $f'(x_n)$

$e_{n+1} = C e_n e_{n-1}$
$|e_{n+1}| \approx A|e_n|^{\frac{1+\sqrt{5}}{2}}$
收敛阶数： $\frac{1+\sqrt{5}}{2}$

7. 常微分方程的数值解法

在给定曲线上一个点的前提下，解一阶常微分方程
不直接求原函数，而是求某个点 $x_i$ 的近似函数值 $y(x_i)$

局部截断误差： $T(n) = O(h^{p+1})$ ，则称该方法是 $p$ 阶的

Euler’s Method

$\frac{y_{n+1} - y_n}{x_{n+1} - x_n} = f(x_n, y_n)$
$y_{n+1} = y_n + f(x_n, y_n)(x_{n+1} - x_n) = y_n + hf(x_n, y_n)$

$p = 1, \text{order } 1$

Trapezoidal Method

implicit Euler formula
$\begin{cases} y_{n+1}^{(0)} = y_n + h f(x_n, y_n) \\ y_{n+1}^{(k+1)} = y_n + \frac{h}{2}[f(x_n, y_n) + f(x_{n+1}, y_{n+1}^{(k)})] \end{cases}$
第二步采用迭代的方法

$p = 2, \text{order } 2$

Modified Euler Method

$\begin{cases} \overline{y}_{n+1} = y_n + h f(x_n, y_n) \\ y_{n+1} = y_n + \frac{h}{2}[f(x_n, y_n) + f(x_{n+1}, \overline{y}_{n+1} )] \end{cases}$
不用迭代

Runge-Kutta Methods

$y_{n+1} = y_n + h \varphi(x_n, y_n, h)$
$\varphi(x_n, y_n, h) = \sum_{i=1}^r c_i K_i$
$K_1 = f(x_n, y_n)$
$K_i = f(x_n + \lambda_i h, y_n + h \sum_{j = 1}^{i-1} u_{ij} K_j), \ i = 2, 3, \cdots, r$

上面的方式是单步(single-step methods)的，因为 $y_{n+1}$ 只依赖 $y_n$

Adams methods

$y_{n+k} = y_{n+k-1} + \sum_{i=0}^{k} \beta_i f_{n+i}$