数值分析笔记
求解方程
二分法
时间复杂度：\(O(\log n)\) n步二分法操作后，得到 \[ \lvert x_c-r \rvert < \frac{b-a}{2^{n+1}}\] 如果误差小于\(0.5 \times 10^{-p}\)，精确到小数点后p位。
定义1.5 令\(e_i\)表示迭代过程中第i步时的误差，如果 \[ \lim_{i \to \infty}\frac{e_{i+1}}{e_i}=S<1\] 该方法被称为满足线性收敛，收敛速度为\(S\)。
定理1.6 假设函数\(g\)是连续可微函数，\(g(r)=r\),\(S=\lvert g'(r) \rvert<1\)，则不动点迭代对于一个足够接近\(r\)的初始估计，以速度\(S\)线性收敛到不动点\(r\)。
根据定理1.6，近似误差之间的关系 \[ e_{i+1} \approx Se_i\] 在收敛过程中得到保持。
定义1.7 如果迭代方法对于一个足够接近\(r\)的初值能收敛到\(r\)，该迭代方法被称为局部收敛到\(r\)。
定理1.6的结论是当\(\lvert g'(r) \rvert<1\)时不动点迭代局部收敛。
如果误差小于 \(0.5 \times 10^{-p}\)，解精确到小数点后 \(p\)位.
牛顿法
用 \(x_i\) 点上切线解作为迭代的下一步 \(x_{i+1}\) \[ x_{i+1}=x_i- \frac{f(x_i)}{f'(x_i)} \]
二次收敛性
定义 令 \(e_i\) 表示第 \(i\) 步的误差，若满足条件 \[ M= \lim_{i \to \infty} \frac{e_{i+1}}{e_i^{2}} < \infty \] 则称迭代二次收敛
定理 令 \(f\) 为二阶连续可导函数，满足 \(f(r)=0\) 和 \(f'(r)\neq 0\)，则 Newton方法局部二次收敛，误差满足 \[ \lim_{i \to \infty}= \frac{e_{i+1}}{e_i^{2}}=M=\frac{f''(r)}{2f'(r)} \]
重根问题和改进的Newton方法
Newton方法不能二次收敛到重根，如 \[ f(x)=x^{m} \]
定理 \(f \in C^{m+1}[a,b]\)，\(f\) 在 \(r\) 点有 \(m\) 阶多重根，则Newton方法局部线性收敛到 \(r\)，常数 \(S=(m-1)/m\)
定理 若已知 \(m>1\) 重根，则改进的Newton方法 \[ x_{i+1}=x_i-\frac{mf(x_i)}{f'(x_i)} \] 二次收敛到 \(r\). 设 \(f(x)=(x-r)^{m}g(x)\)，则 \[ \lim_{i \to \infty} \frac{e_{i+1}}{e_{i}^{2}}=\frac{g'(r)}{mg(r)} \]
或者可以通过计算 \(\displaystyle \mu(x)=\frac{f(x)}{f'(x)}\) 的根（注意它没有重根）来得到 \(f(x)\) 的根.
割线法
思想是用差商代替Newton法中的导数. 取 \(x_0\)为初始估计 \[ x_{i+1}=x_i-\frac{f(x_i)(x_i-x_{i-1})}{f(x_i)-f(x_{i-1})} \]
它满足超线性收敛 \[ e_{i+1} \thickapprox \lvert \frac{f''(r)}{2f'(r)} \rvert ^{\alpha-1}e_i^{\alpha} \] 其中 \(\alpha=(1+\sqrt{5})/2\)
一般当 \(1<\alpha<2\) 时都称为超线性收敛.
事实上对于割线法的收敛阶有（当 \(f'(r),f''(r)\neq 0\)） \[ e_{i+1}= \lvert f''(r)/2f'(r) \rvert e_i e_{i-1} \] 成立. 这可以推出上面的 \(\alpha\).
插值
数据和插值函数
对 \(n\) 个不同点，存在唯一一个不高于 \(n-1\) 次的插值多项式
方便的做法是利用Newton差商. \(f[x_1\cdots x_n]\) 表示（唯一）多项式的 \(x^{n-1}\) 项的系数，有递归式 \[ 略 \] 或者 \[ f[x_1\cdots x_k]= \sum_{j=1}^{k} \frac{f(x_j)}{(x_j-1)\cdots (x_j-x_{j-1})(x_j-x_{j+1})\cdots (x_j-x_k)} \] 从而 \[ f[x_1\cdots x_n]=f[\sigma(x_1)\cdots \sigma(x_n)] \] 如果 \(f(x)\) 在区间内有 \(n-1\) 阶导数，则在区间内存在一点 \(c\)使得 \[ f[x_1\cdots x_n]=\frac{f^{(n-1)}(c)}{(n-1)!} \]
Newton插值多项式 \[ P(x)=f[x_1]+f[x_1 x_2](x-x_1)+\cdots +f[x_1\cdots x_n](x-x_1)\cdots (x-x_{n-1}) \]
插值误差
\[ f(x)-P_{n-1}(x)=f[x_1\cdots x_nx](x-x_1)\cdots (x-x_n)= \frac{f^{(n)}(c)}{n!}(x-x_1)\cdots (x-x_n) \]
龙格现象
略，大一上写过文章.
Hermite插值
\(f(x) \in C^{1}[a,b]\)，\(a\leqslant x_0<x_1<\cdots <x_n\leqslant b\). 则 \[ H_{2n+1}(x)= \sum_{i=0}^{n} f(x_i)A_i(x)+ \sum_{i=0}^{n} f'(x_i)B_i(x) \] \[ A_i(x)=[1-2(x-x_i)l_i'(x_i)]l_i^{2}(x) \] \[ B_i(x)=(x-x_i)l_i^{2}(x) \] 其中 \(l_i\) 是Lagrange插值基.
Chebyshev插值
希望找到最佳的一致逼近
取 \[ x_i= \cos \frac{(2i-1)\pi}{2n}, i=1,\cdots,n \] 则此时 \[ \max_{-1\leqslant x\leqslant 1} \lvert (x-x_1)\cdots (x-x_n) \rvert = \frac{1}{2^{n-1}} T_n(x)\leqslant \frac{1}{2^{n-1}} \] 取到最小值
Chebyshev多项式 满足 \(T_n(x)=\cos (n \arccos x)\)，且有递推关系 \[ T_0(x)=1, \quad T_1(x)=x \] \[ T_{n+1}(x)=2x T_n(x)- T_{n-1}(x) \] - \(T_n(x)\) 为 \(n\) 次多项式. - 最高次项系数为 \(2^{n-1}\). \(T_n(1)=1, T_n(-1)=(-1)^{n}\). - 最大绝对值为 \(1\). - \(T_n(x)=\cos (n \arccos x)\) 的极值点为 \[ x_i= \cos \frac{i\pi}{n}, i=0,\cdots ,n \]
令 \(\displaystyle y=\frac{b-a}{2}x+\frac{b+a}{2}, x\in [-1,1]\)，则 \[ \lvert (y-y_1)\cdots (y-y_n) \rvert =(\frac{b-a}{2})^{n}\lvert (x-x_1)\cdots (x-x_n) \rvert \leqslant (\frac{b-a}{2})^{n} \frac{1}{2^{n-1}} \]
三次样条插值
\(n\) 个点 \((x_i,y_i)\) 一共 \(3n-5\) 个方程，\(3n-3\) 个未知数. 需要补充两个条件. 具体计算见p152-153.
- 自然三次样条：\(S_1''(x_1)=S_{n-1}''(x_n)=0\).
- 曲率调整样条：\(S_1''(x_1)=2c_1=v_1, S_{n-1}''(x_n)=2c_n=v_n\).
- 钳制边界条件：\(S_1'(x_1)=v_1, S_{n-1}'(x_n)=v_n\).
- 抛物线端点的三次样条：\(d_1=0=d_{n-1}\).
- 非纽结三次样条：\(S_1=S_2,S_{n-2}=S_{n-1}\)
误差：设 \(f(x) \in C^{4}[a,b]\)，对于曲率调整样条或钳制样条 \(S(x)\)，令 \(\delta = \max_{i} \delta_i\)，有估计式 \[ \max_{a\leqslant x\leqslant b} \lvert f(x)- S(x) \rvert \leqslant \frac{5}{384} \max_{a\leqslant x\leqslant b} \lvert f^{(4)}(x) \rvert \delta^{4} \]
Bezier曲线
最小二乘
离散
给定函数 \(f\)，求 \(P_n(x)\) 使误差 \[ E=\frac{1}{m} \sum_{i=1}^{m} (y_i-P_{n-1}(x_i))^{2} \] 最小 设 \[ A= \begin{bmatrix} 1 & x_1 & \cdots & x_1^{n-1} \\ \vdots & & & \vdots \\ 1 & x_{m} & \cdots & x_m^{n-1} \\ \end{bmatrix} \] 则法线方程（\(a=(a_0,a_1,\cdots ,a_{n-1})^{\mathsf{T}},y=(y_1,\cdots ,y_m)^{\mathsf{T}}\)） \[ A^{\mathsf{T}}A a= A^{\mathsf{T}}y \]
给定 \((t_1,y_1),\cdots ,(t_n,y_n)\)，则平方误差最小逼近的直线 \(y=a_0+a_1t\)满足 \[ \begin{bmatrix} \sum_{i=1}^{m} 1 & \sum_{i=1}^{m} t_i \\ \sum_{i=1}^{m} t_i & \sum_{i=1}^{m} t_i^{2} \\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \end{bmatrix} = \begin{bmatrix} \sum_{i=1}^{m} y_i \\ \sum_{i=1}^{m} t_i y_i \\ \end{bmatrix} \]
有时 \(A^{\mathsf{T}}A\) 的条件数过高，可以考虑QR分解. 设 \(A=QR\)， \[ \left\| QRx-b \right\|_{2}= \left\| Rx-Q^{\mathsf{T}}b \right\|_{2} \] 考虑到 \[ R=\begin{bmatrix} R_1 \\ 0 \\ \end{bmatrix} \] 再令 \[ Q^{\mathsf{T}}b=\begin{bmatrix} b_1 \\ b_2 \\ \end{bmatrix} \]
而 \(R_1\) 是上三角方阵，所以法线方程即为 \[ R_1 x=b_1 \]
数据线性化 \[ y=c_1 \mathrm{e}^{c_2t} \Rightarrow \ln y= \ln c_1+ c_2t \]
连续
在一般的Euclid空间中，给定函数 \(f\)，求 \(P_n(x)\) 使误差 \[ E= \int_{a}^{b} [f(x)-P_n(x)]^{2} \mathrm{d}x \] 最小，称 \(P_n(x)\) 为该函数的最小二乘逼近多项式，假设 \(P_n(x)=a_nx^{n}+\cdots +a_1x+a_0\)，那么法线方程 \[ \frac{\partial E}{\partial a_j}=0 \] 即为 \[ \sum_{k=0}^{n} a_k \int_{a}^{b} x^{k+j} \mathrm{d}x= \int_{a}^{b} f(x) x^{j} \mathrm{d}x \] 如果我们令 \(\varphi_k=x^{k}\) 为基函数，并赋予Euclid空间的典范内积则法线方程变为 \[ \begin{bmatrix} (\varphi_0,\varphi_0) & \cdots & (\varphi_0,\varphi_n) \\ \vdots & & \vdots \\ (\varphi_n,\varphi_0) & \cdots & (\varphi_n,\varphi_n) \\ \end{bmatrix} \begin{bmatrix} a_0 \\ \vdots \\ a_n \\ \end{bmatrix} = \begin{bmatrix} (f,\varphi_0) \\ \vdots \\ (f,\varphi_n) \\ \end{bmatrix} \] 比如如果我们用 \(\varphi_k\) 求 \(f(x)=\sin \pi x\) 的最小二乘逼近多项式，那么得到的系数矩阵是Hilbert矩阵，条件数很大，求解困难. 通常需要采用正交基.
权函数和正交基
权函数 \(w(x)\) 满足 \(\forall x \in I=[a,b],(w(x)\geqslant 0)\)，且对于 \(I\) 上的任意子区间不恒为零. 从而可以定义内积 \[ (f,g)_{w}= \int_{a}^{b} w(x)f(x)g(x) \mathrm{d}x \] 如果 \(\{\varphi_0,\cdots ,\varphi_n\}\) 满足 \[ \int_{a}^{b} w(x)\varphi_j(x)\varphi_k(x) \mathrm{d}x= \begin{cases} 0,\quad j\neq k \\ \alpha_k>0,\quad j=k \end{cases} \] 那么它们就是一组正交基，如果 \(\alpha_k=1\)，则标准正交.
最小二乘问题
给定函数 \(f\)，\(w(x)\) 为权函数，求 \(P_n(x)=\sum_{k=0}^{n} a_k\varphi_k(x)\) 使误差 \[ E(a_0,\cdots ,a_n)=\int_{a}^{b} w(x)[f(x)-P_n(x)]^{2} \mathrm{d}x \] 极小，称 \(P_n(x)\)为该函数以 \(w(x)\) 为权函数的最小二乘函数逼近.
如果 \(\{\varphi_0,\cdots ,\varphi_n\}\) 为一组正交函数系，则 \[ a_k=\frac{(\varphi_k,f)_{w}}{(\varphi_k,\varphi_k)_{w}} \]
Gram-Schmidt 正交化过程
是想把 \(\{1,x,\cdots ,x^{n}\}\rightarrow \{\varphi_0,\cdots ,\varphi _n\}\) 成为正交基.
设 \(w(x)\)为 \([a,b]\)上的权函数. 令 \(\varphi_0(x)=1\)，\(\varphi_1(x)=x-B_1\)，其中 \[ B_1=\frac{(x\varphi_0,\varphi_0)_{w}}{(\varphi_0,\varphi_0)_{w}} \] 对于权函数 \(w(x)\equiv 1\)的情况 \[ B_1=\frac{\int_{-1}^{1} x \mathrm{d}x}{\int_{-1}^{1} 1 \mathrm{d}x}=0 \]
当 \(k\geqslant 2\)时，有 \[ \varphi_k(x)=(x-B_k)\varphi_{k-1}-C_k\varphi_{k-2} \] 其中 \[ B_k=\frac{(x\varphi_{k-1},\varphi_{k-1})_{w}}{(\varphi_{k-1},\varphi_{k-1})_{w}}, \qquad C_k=\frac{(x\varphi_{k-2},\varphi_{k-1})_{w}}{(\varphi_{k-2},\varphi_{k-2})_{w}} \] 则所得到的集合 \(\{\varphi_0,\cdots ,\varphi_n\}\) 是正交的. 按这种办法求出来的多项是首一的，可能与一般定义的Legendre多项式或Chebyshev多项式不同.（标准定义的Legendre多项式要求 \(P_n(1)=1\)）
首一性还可以得到以下性质 \[ (x\varphi_{k-2},\varphi_{k-1})_{w}=(\varphi_{k-1},\varphi_{k-1})_{w} \]
QR分解
完全QR分解 对于 \(A\in \mathbb{F}^{m \times n}, m > n\), \(A=QR\). 其中 \(Q \in O(m)\) \[ R= \begin{bmatrix} \tilde{R} \\ 0 \\ \end{bmatrix} \]
计算的复杂度为 \(O(n^{3})\)
Gram-Schimidt 正交
朴素的和改进的版本p194
过程需要大约 \(m^{3}\) 次乘法和 \(m^{3}\) 次加法：
计算 \(r_{ij}=q_{i} ^{\mathsf{T}}A_j\) 时需要进行 \(m\)次乘法和 \(m-1\)次加法，计算 \(y-y-r_{ij}q_i\)时进行 \(m\)次乘法和 \(m\)次加减法；舍去计算 \(r_{jj}\)和 \(q_j\)的次数. 则加减法次数和乘除法次数都约为 \[ \sum_{j=1}^{n} (j-1)\cdot 2m=m^{3}-m \] 即需要大约 \(m^{3}\)次乘法和 \(m^{3}\)次加法.
Householder 反射
令 \(w \in \mathbb{R}^{n}\)，满足 \(w^{\mathsf{T}}w=1\)，则 \[ H=I-2ww^{\mathsf{T}} \] 为一个 Householder 反射，它对称正交
设 \(\left\| x \right\|_{2}= \left\| y \right\|_{2}\)，取 \(\displaystyle w=\frac{x-y}{\left\| x-y \right\|_{2}}\)，则 \(H=I-2ww^{\mathsf{T}}\) 满足 \(Hx=y\)
其实每步选模长的正负号时应该与原位置的正负号相反，书上忽略了这一点.
Householder变换做QR分解的复杂度
只考虑向量与向量相乘，矩阵与向量相乘，矩阵与矩阵相乘的计算 对 \(m \times n\) 的矩阵 \(A\) 进行QR分解， - 若只需要R，则大约需要 \[ n^{2}(m-\frac{n}{3}) \] 次乘法和加减法.
- 若还需要Q，则大约需要 \[ m^{2}n-mn^{2}+\frac{n^{3}}{3} \] 次乘法和加减法.
QR decompositon - Wikipedia https://en.wikipedia.org/wiki/QR_decomposition
数值微分和积分
数值微分
有限差分公式
二点前向差分公式 \[ f'(x) =\frac{f(x+h)-f(x)}{h} - \frac{h}{2} f''(c) \] 其中 \(c\) 在 \(x\) 和 \(x+h\) 之间
如果误差是 \(O(h^{n})\)，称公式为 \(n\) 阶近似
三点中心差分公式（二阶中心差分公式） \[ f'(x)= \frac{f(x+h)-f(x-h)}{2h}- \frac{h^{2}}{6}f'''(c) \] 其中 \(x-h<c<x+h\)
五点中心差分公式 \[ f'(x)\thickapprox \frac{f(x-2h)-8f(x-h)+8f(x+h)-f(x+2h)}{12h} \] 误差 \(O(h^{4})\)
二阶导数的三点中心差分公式 \[ f''(x)= \frac{f(x+h)-2f(x)+f(x-h)}{h^{2}}-\frac{h^{2}}{12} f^{(4)}(c) \] 其中 \(x-h<c<x+h\)
舍入误差
以三点中心差分公式为例，加上机器误差，总误差的绝对值上界是 \[ \frac{h^{2}}{6} f'''(c)+ \frac{\varepsilon_{mach}}{h} \]
外推
如果有 \[ Q=F_n(h) + Kh^{n} +O(h^{n+1}) \] 那么用 \(\frac{h}{2}\) 代替 \(h\)，重新解出得 \[ Q\thickapprox \frac{2^{n}F(h/2)-F(h)}{2^{n}-1} \] \(F_{n+1}(h)\) 至少是 \(n+1\) 阶近似 \(Q\) 的公式
一个例子是用二阶中心差分公式进行外推，得到的其实是四阶公式. 且误差项仅和 \(h\) 的偶数阶有关.
数值积分
基本上都是用插值公式和误差导出积分的公式和误差.
梯形法则
\[ \int_{x_0}^{x_1} P(x) \mathrm{d}x= \frac{h}{2}(f(x_0)+f(x_1)) \] \[ \int_{x_0}^{x_1} E(x) \mathrm{d}x= -\frac{h^{3}}{12}f''(c) \]
Simpson法则
\[ \int_{x_0}^{x_2} f(x) \mathrm{d}x=\frac{h}{3}(y_0+4y_1+y_2)-\frac{h^{5}}{90}f^{(4)}(c) \] 其中 \(h=x_2-x_1=x_1-x_0\)，\(c\)在 \(x_0\)和 \(x_2\)之间.
误差的证明
首先误差 \[ E_s= \int_{x_0}^{x_2} \frac{(x-x_0)(x-x_1)(x-x_2)}{3!}f'''(c(x)) \mathrm{d}x \] 因为会变号，所以不能使用积分中值定理，尝试改写 \(f'''(c(x))\) 得到
\[ f'''(c(x))=f'''(c(x_1))+(x-x_1)f^{(4)}(c(\bar{x}))c'(\bar{x}) \] 代入 \(E_s\) 后 \(f'''(c(x_1))\) 对应项积分为零，这可以解释为什么误差是跟 \(f^{(4)}(x)\) 有关.
另一方面，用Newton插值，可以认为再在 \(x_1\) 点插值 \[ P_3(x)=P_2(x)+D(x-x_0)(x-x_1)(x-x_2) \] 那么 \[ \int_{x_0}^{x_2} P_2(x) \mathrm{d}x=\int_{x_0}^{x_2} P_3(x) \mathrm{d}x \] 而由插值误差即知 \[ E_s=\int_{x_0}^{x_2} \frac{(x-x_0)(x-x_1)^{2}(x-x_2)}{4!}f^{(4)}(c) \mathrm{d}x \]
这下可以使用积分中值定理了，计算得到 \[ E_s=-\frac{h^{5}}{90} f^{(4)}(c) \]
从该表达式知Simpson法则的代数精度为3.
3阶Newton-Cotes公式（Simpson3/8公式）
\[ \int_{x_0}^{x_3} f(x) \mathrm{d}x \thickapprox \frac{3h}{8}(f_0+3f_1+3f_2+f_3) \] 其误差也为 \(O(h^{5})\)，而因计算更复杂而通常不使用.
复合梯形法则
对于 \(\displaystyle x_j=x_0+jh, x_0=a, x_n=b, h=(b-a)/n\)，如果 \(f''\) 在区间 \([a,b]\) 连续，有 \[ \int_{a}^{b} f(x) \mathrm{d}x=\frac{h}{2}(f_0+2f_1+\cdots +2f_{n-1}+f_n)-\frac{(b-a)h^{2}}{12}f''(c) \]
数值精度是2阶，代数精度是1阶.
复合Simpson公式
如果 \(f^{(4)}\) 在 \([a,b]\) 上连续，有 \[ \int_{a}^{b} f(x) \mathrm{d}x=\frac{h}{3}(y_0+y_{2m}+4\sum_{i=1}^{m} y_{2i-1}+2\sum_{i=1}^{m-1} y_{2i})-\frac{(b-a)h^{4}}{180}f^{(4)}(c) \]
注意使用 \(m\) 次Simpson公式需要 \(2m\) 等分. \(h=(b-a)/2m\)
舍入误差对复合积分公式的影响
\[ \int_{a}^{b} f(x) \mathrm{d}x\thickapprox \frac{h}{2}(f_0+2f_1+\cdots +2f_{n-1}+f_n)+\varepsilon_R \] 有 \[ \varepsilon_R \thickapprox (b-a)\varepsilon_h \]
中点法则
当 \(f \in C^{2}[a,b]\) 时， \[ \int_{x_0}^{x_1} f(x) \mathrm{d}x=hf(w)+\frac{h^{3}}{24}f''(c) \] 其中 \(h=x_1-x_0\)，\(w\) 是中点 \(x_0+h/2\)，\(x_0\leqslant c\leqslant x_1\)
Taylor展开至二阶即可得到误差. 与梯形法则相比，中点法则的误差只有一半.
复合中点法则
\[ \int_{a}^{b} f(x) \mathrm{d}x= h\sum_{i=1}^{m} f(w_i)+\frac{(b-a)h^{2}}{24}f''(c) \] 其中 \(h=(b-a)/m\)，\(w\) 是 \([a,b]\) 中 \(m\) 个相等子区间的中点
另一个开Newton-Cotes法则
\[ \int_{x_0}^{x_4} f(x) \mathrm{d}x=\frac{4h}{3}[2f(x_1)-f(x_2)+2f(x_3)]+\frac{14h^{5}}{45}f^{(4)}(c) \]
Romberg积分
一个表格 \[ \begin{aligned} &R_{11} \\ &R_{21} \quad R_{22} \\ &R_{31} \quad R_{32} \quad R_{33} \\ &\vdots \end{aligned} \] 定义 \[ h_j=\frac{1}{2^{j-1}}(b-a) \]
其中 \[ R_{j1}=\frac{1}{2}R_{j-1,1}+h_j\sum_{i=1}^{2^{j-2}} f(a+(2i-1)h_j) \]
\[ R_{jk}=\frac{4^{k-1}R_{j,k-1}-R_{j-1,k-1}}{4^{k-1}-1} \]
这里 \(R_{jj}\)是 \(2j\)阶的近似.
Romberg积分第二列最后一项和复合Simpson法则的结果一致. 事实上，Romberg积分的第一列定义为连续的复合梯形法则，第二列是复合Simpson项，即复合梯形法则的外推是复合Simpson法则.
Romberg积分的常用停止条件是计算新的一行直到相邻的对角线元素 \(R_{jj}\)差异小于当前的容差.
自适应积分
为什么要自适应： - 高阶导数不易计算 - 有些区域变化缓慢
高斯积分
高斯积分具有 \(2n+1\)阶的代数精度（使用 \(n+1\)个点），该精度是Newton-Cotes公式的两倍.
三点Gauss积分 \[ x_1=-\sqrt{\frac{3}{5}} \quad x_2=0 \quad x_3=\sqrt{\frac{3}{5}} \\ c_1=\frac{5}{9} \quad c_2=\frac{8}{9} \quad c_3=\frac{5}{9} \]
关于正交多项式系，有这样的定理：
如果 \(\{p_0,p_1,p_2,\cdots,p_n\}\)在 \([a,b]\)区间上是多项式的正交集，并且如果deg \(p_i=i\) ，则 \(p_i\)在区间 \((a,b)\)上有 \(i\)个不同的根
证明可以参考
多项式以及相关内容 https://zhuanlan.zhihu.com/p/270620896
| ### 正交多项式 |
| Legendre多项式：在 \([-1,1]\)上以权函数 \(w(x)=1\)的正交多项式. 它们满足递推关系 \[ (n+1)P_{n+1}(x)=(2n+1)x P_n(x)-nP_{n-1}(x) \] 或者 \[ p_i(x)=\frac{1}{2^{i}i!} \frac{\mathrm{d}^{i}}{\mathrm{d}x^{i}}[(x^{2}-1)^{i}] \] 其中 \(P_0(x)=1\)，\(P_1(x)=x\) |
| 这里Legendre多项式非首一，而是满足 \(P_n(1)=1\). |
| 正交性： \[ (P_n,P_m)=\int_{-1}^{1} P_n(x)P_m(x) \mathrm{d}x= \begin{cases} 0, \quad m \neq n \\ \frac{2}{2n+1}, m = n \end{cases} \] |
| 奇偶性：\(P_{2k}(x)\) 只含有偶次幂， \(P_{2k+1}\) 只含有奇次幂 |
| 若令 \(\tilde{P}_n(x)= \frac{2^{n}(n!)^{2}}{(2n)!}P_n(x)\)，则称 \(\tilde{P}_n(x)\) 为首一的Legendre多项式. |
| 更多关于正交多项式还可以参考 |
逆问题解析 https://zhuanlan.zhihu.com/p/352724194
还有一类正交多项式Laguerre多项式，其对应的权函数 $w(x)=^{-x} $
接下去我们要证明插值节点恰为Legendre多项式的根，这跟首项系数无关.
高斯积分
我们想用 \(n\) 阶插值多项式尽可能拟合高阶多项式的积分.
\[ \int_{-1}^{1} x^{k} \mathrm{d}x=\sum_{i=1}^{n} c_i x_i^{k} \] 一共有 \(x_i,c_i,i=1,2,\cdots,n\) 这 \(2n\) 个未知数，因此我们希望 \(k=0,1,\cdots ,2n-1\)，即 Gauss 积分的代数精度为 \(2n-1\). 事实上有
\[ \int_{-1}^{1} 1 \cdot f(x) \mathrm{d}x \thickapprox \sum_{i=1}^{n} c_if(x_i) \] 其中 \[ c_i=\int_{-1}^{1} 1 \cdot L_i(x) \mathrm{d}x, \quad i=1,\cdots ,n \] \(L_i\)是Lagrange插值基. \(x_i\)是Legendre多项式的根. 写出1是为了强调权函数.
定理：Gauss积分方法，在 \([-1,1]\)上使用 \(n\)阶Legendre多项式，具有 \(2n-1\)阶代数精度.
并不需要解那个很复杂的非线性方程组，只要说明Gauss积分方法能够具有 \(2n-1\)阶代数精度，就可由唯一性得知它就是最好的插值法. 然而事实上Lengdre多项式的根确实是那个很复杂的非线性方程组的根，这一点很难直接证明.
在一般区间 \([a,b]\)上近似积分，作变换 \(t=(2x-a-b)/(b-a)\)，容易检验 \[ \int_{a}^{b} f(x) \mathrm{d}x=\int_{-1}^{1} f(\frac{(b-a)t+b+a}{2})\frac{b-a}{2} \mathrm{d}t \]
Gauss积分的误差： \[ R=\int_{-1}^{1} f(x) \mathrm{d}x-\sum_{i=1}^{n} c_if(x_i)=\frac{f^{(2n)}(\xi)}{(2n)!}\int_{-1}^{1} q(x) \mathrm{d}x \] 其中 \(q(x)=[(x-x_1)\cdots (x-x_n)]^{2}\)
可以考虑 \(x_1,x_2,\cdots,x_n\)的Hermite插值.
Gauss-Chebyshev 积分
权函数 \(\displaystyle w(x)=\frac{1}{\sqrt{1-x^{2}}}\)，实际上是换了一个Hilbert空间.
此时希望 \[ \int_{-1}^{1} f(x) \mathrm{d}x= \int_{-1}^{1} w(x)g(x) \mathrm{d}x \thickapprox \sum_{i=1}^{n} c_i g(x_i) \] 其中 \(g(x)= \sqrt{1-x^{2}}f(x)\). 这时我们得到的正交多项式为Chebyshev多项式 \(p_i(x)\)，满足 \[ \int_{-1}^{1} \frac{1}{\sqrt{1-x^{2}}}p_m(x)p_n(x) \mathrm{d}x=0, m\neq n \] 诸零点为 \(\displaystyle x_i=\cos \frac{2i-1}{2n}\pi,i=1,2,\cdots ,n\)，而 \(\displaystyle c_i=\int_{-1}^{1} \frac{L_i(x)}{\sqrt{1-x^{2}}} \mathrm{d}x= \frac{\pi}{n}\)
反常积分的计算
尝试利用数值方法计算 \[ \int_{a}^{b} \frac{g(x)}{(x-a)^{p}} \mathrm{d}x \] 其中 \(0<p<1\)
该积分有一个奇点 \(x=a\). - 梯形和Simpson方法发散 - 开Newton-Cotes公式精度很低
最简单的想法 \[ \int_{a}^{b} \frac{g(x)-g(a)}{(x-a)^{p}} \mathrm{d}x+ \int_{a}^{b} \frac{g(a)}{(x-a)^{p}} \mathrm{d}x \] 这样第一项的被积函数是连续函数. 可以使用之前说到的积分方法. 当然我们可以做得更加精细：Simpson 将 \(g(x)\) Taylor 展开至4阶，这样 \(\displaystyle G(x)=\frac{g(x)-Q_4(x)}{(x-a)^{p}}\in C^{4}[a,b]\) 就可以利用之前说过的公式计算. （ \(Q_4(x)\) 为 \(g(x)\) 在 \(x=a\) 处Taylor展开的前5项. 而 \[ \int_{a}^{b} \frac{Q_4(x)}{(x-a)^{p}} \mathrm{d}x= \sum_{k=0}^{4} \frac{g^{(k)}(a)}{k!}\frac{(b-a)^{k+1-p}}{k+1-p} \] 令 \[ G(x)= \frac{g(x)-Q_4(x)}{(x-a)^{p}} \] 则相应误差 \[ E=-\frac{h^{4}}{180}(b-a) G^{(4)}(\xi) \]
高维积分的计算
这会牵扯到许多极限换序和收敛性的问题.
特征值与奇异值
Gerschgorin 圆盘定理 略
幂迭代法
幂迭代法线性收敛
Rayleigh商 如果有近似特征向量 \(x\) 那么对于特征值的最优近似是 \[ \lambda= \frac{x ^{\mathsf{T}}Ax}{x ^{\mathsf{T}}x} \]
定理：设 \(A\) 是特征值为 $_1 > _2 _m $ 的 \(m \times m\) 矩阵. \(A\) 有完全特征向量系. 则对于几乎所有的初始向量，幂迭代线性收敛到和 \(\lambda_1\) 相关的特征向量，收敛常数为 $S=_2/_1 $.
先取逆再进行幂迭代可以得到绝对值最小的特征值（注意特征向量不变）. 为了避免对矩阵 \(A\) 的逆矩阵的显式求解，将 \(x_{k+1}= A ^{-1} x_k\) 等价为 \(Ax_{k+1}=x_k\) 并使用Gauss消去法求解.
逆向幂迭代
p469
Rayleigh商迭代（RQI）
可利用Rayleigh商估计最大特征值，Rayleigh商迭代对称矩阵三次收敛，一般则二次收敛.
QR算法
- 利用Householder变换 \(B=HAH\) 作相似约化：对称阵约化为三对角，非对称阵约化为上Hessenberg矩阵
- 利用QR分解构造迭代序列. 分解 \(A_j=Q_jR_j\)，计算 \(A_{j+1}=R_{j}Q_{j}\)
约化到上Hessenberg矩阵大约 \(\displaystyle \frac{5}{3} n^{3}\) 次乘法. 对于 \(n\) 阶矩阵 \(A\)需要作 \(n-2\) 次Householder变换将 \(A\) 变为上Hessenberg形式.
SVD
\(A\) 为 \(m \times n\) 矩阵，\(\{u_1,\cdots ,u_m\}\) 和 \(\{v_1,\cdots ,v_n\}\) 是单位正交集合，若存在 \(s_1\geqslant s_2\geqslant \cdots \geqslant s_n\geqslant 0\)，满足 \[ Av_1=s_1u_1 \\ \vdots \\ A v_n= s_n u_n \] 称 \(v_i\) 为 \(A\) 的右奇异向量， \(u_i\) 为 \(A\) 的左奇异向量， \(s_i\) 为 \(A\) 的奇异值.
对称矩阵的SVD简化为 \(A= U ^{\mathsf{T}} SV\)，其中 \(U\) 的每一列与 \(V\) 的对应列相同或相反.
低秩逼近
\(A\) 的最优最小二乘逼近为保留前面 \(p\) 项，\(p<r\).
另一种SVD的方法
考虑 \[ B=\begin{bmatrix} 0 & A ^{\mathsf{T}} \\ A & 0 \\ \end{bmatrix} \] 设有 \(B\) 的特征向量 \[ \begin{bmatrix} A ^{\mathsf{T}}w \\ Av \\ \end{bmatrix} = \begin{bmatrix} 0 & A^{\mathsf{T}} \\ A & 0 \\ \end{bmatrix} \begin{bmatrix} v \\ w \\ \end{bmatrix} =\lambda \begin{bmatrix} v \\ w \\ \end{bmatrix} \] 注意 \(w\) 是 \(A ^{\mathsf{T}}A\) 的特征向量，对应的特征值为 \(\lambda^{2}\).
其他SVD相关性质
Frobenius范数：\(\left\| A \right\|_{F}= (\sum_{i}^{} \sum_{j}^{} a_{ij})^{1/2}\)，则 \[ \left\| A \right\|_{F}^{2}= s_1^{2}+\cdots +s_{r}^{2} \] \[ \left\| A \right\|_{2}=s_1 \]
设 \(U^{\mathsf{T}}AV= diag\{s_1,\cdots ,s_r\}\)，若 \(k<r\)，令 \(A_k=\sum_{i=1}^{k} s_iu_iv_i^{\mathsf{T}}\)，则 \[ \min_{\operatorname{rank}(B)=k}\left\| A-B \right\|_{2}= \left\| A-A_k \right\|_{2}= s_{k+1} \]
数值代数
方程组的求解——高斯消去法和LU分解
朴素的Gauss消去法
\(n\)个方程 \(n\)个未知数的消去计算，可以在 \(\displaystyle \frac{2}{3}n^{3}+\frac{1}{2}n^{2}-\frac{7}{6}n\)次操作后完成. \(\displaystyle \frac{1}{3}n^{3}\) 次操作是加减， \(\displaystyle \frac{1}{3}n^{3}\) 次操作是乘除. 所有操作完成后，变为上三角矩阵.
再经历回代的操作，回代计算复杂度是 \(O(n^{2})\)（实际上就是 \(n^{2}\) 次），加减与乘除也是一半一半.
如果某时主元为0，会导致Gauss消去法终止.
LU分解
\(A=LU\) \[ Ly=b \Rightarrow Ux=y \]
如果 \(L\) 和 \(U\) 的对角线元素满足不同条件，也是不同的分解. - \(l_{ii}=1\) 为 Doolittle 分解 - \(u_{ii}=1\) 为 Crount 分解 - \(l_{ii}=u_{ii}\) 为 Choleski 分解
一般 LU 分解的复杂度是 \(O(n^{3})\)，而有限带宽矩阵的 LU 分解复杂度为 \(O(n)\).
LU分解相比经典的 Gauss 消去法，更适合求解 \(b\) 变化时的问题
矩阵范数
- 算子范数
- 无穷范数 \[ \left\| A \right\|_{\infty}= \max_{\left\| x \right\|_{\infty}=1} \left\| Ax \right\|_{\infty} \] 有 \(\left\| A \right\|_{\infty}= \max_{1\leqslant i\leqslant n} \sum_{j=1}^{n} \lvert a_{ij} \rvert\)
- \(L_2\) 范数 \[ \left\| A \right\|_{2}= \max_{\left\| x \right\|_{2}=1} \left\| Ax \right\|_{2} \] 有 \(\left\| A \right\|_{2}= [\rho(A^{\mathsf{T}}A)]^{1/2}\)；\(\rho(A)\leqslant \left\| A \right\|_{}\) 对任意满足相容性的范数成立
误差放大和条件数
令 \(x_{a}\) 为 \(Ax=b\) 的近似解，后向误差为余项的范数 \(\left\| b-Ax_{a} \right\|_{\infty}\)，前向误差为 \(\left\| x-x_{a} \right\|_{\infty}\).
定义 \[ 误差放大因子Q=\frac{相对前向误差}{相对后向误差}=\frac{\left\| x-x_{a} \right\|_{\infty}/ \left\| x \right\|_{\infty}}{\left\| r \right\|_{\infty}/ \left\| b \right\|_{\infty}} \]
方阵 \(A\) 的条件 cond(\(A\)) 为最大的误差放大因子.
\(n\) 阶矩阵 \(A\) 的条件数为 \[ \text{cond}(A)= \left\| A \right\|_{} \left\| A^{-1} \right\|_{} \] 一般可以取无穷范数
线性方程组-迭代法 0.2：条件数 https://zhuanlan.zhihu.com/p/388053510
Q: 证明矩阵 \(A\)的条件数为 \(\text{cond}(A)=\left\| A \right\|_{}\left\| A ^{-1} \right\|_{}\)
证: 注意误差放大因子 \[ Q=\frac{\left\| x-x_a \right\|_{}/\left\| x \right\|_{}}{\left\| b- Ax_a \right\|_{}/\left\| b \right\|_{}}\leqslant \]
\[ \exists b,x,r \quad \text{s.t.} \quad \left\| A \right\|_{}=\frac{\left\| Ax \right\|_{}}{\left\| x \right\|_{}}=\frac{\left\| b \right\|_{}}{\left\| x \right\|_{}}, \quad \left\| A ^{-1} \right\|_{}=\frac{\left\| A ^{-1} r \right\|_{}}{\left\| r \right\|_{}} \]
PA=LU 分解
淹没问题 p80
部分主元 第一步选择第 \(p\) 行，其中 \[ \lvert a_{p1} \rvert \geqslant \lvert a_{i1} \rvert ,1\leqslant i\leqslant n \] 交换第 \(1\) 行和第 \(p\) 行，然后再作消去. 以后每一步都要比较.
然后 \[ Lc=Pb \Rightarrow Ux=c \]
迭代方法
不动点迭代：\(Ax=b\) 等价于 \(x= Tx+b\)
设 \(A=D+L+U\)，则有 - Jacobi 迭代中的 \(T\): \[ x=D ^{-1} [b-(L+U)x] \] - Gauss-Seidel迭代中的 \(T\)： \[ x=(L+D) ^{-1}(b-Ux) \]
连续松弛迭代（SOR） - 引入参数 \(\omega\)，对上面两种迭代方法加权组合 \[ [\omega(D+L)+(1-\omega)D]x=\omega(-Ux+b)+(1-\omega)Dx \] 当 \(0<\omega<1\) 时称为欠松弛，\(\omega>1\) 称为过松弛
迭代基本定理 设 \(x= Bx+f\)，对任意的 \(x^{(0)}\)和 \(f\)，迭代法 \(x^{(k+1)}=Bx^{(k)}+f\) 收敛 \(\Leftrightarrow\) \(\rho(B)<1\).
收敛速度： \(\varepsilon^{(k)}=B^{k}\varepsilon^{(0)}\) 设 \(B\)有 \(n\)个线性无关的特征向量 \(u_1,\cdots ,u_n\)，对应特征值为 \(\lambda_1,\cdots ,\lambda_n\)，则由 \(\varepsilon^{(0)}=\sum_{i=1}^{n} a_iu_i\),
\[ \varepsilon^{(k)}=\sum_{i=1}^{n} a_i\lambda_i^{k}u_i \] 可以看到 \(\rho(B)<1\) 越小，收敛越快.
定义： \[ R(B)=-\ln \rho(B) \] 称为迭代法的收敛速度
定理： \(x^{(k+1)}=Bx^{(k)}+f\)为迭代公式，存在某一范数（需要满足相容性条件）使得 \(\left\| B \right\|_{}=q<1\)，则 - 迭代法收敛 - \(\displaystyle \left\| x^{(k)}-x^{*} \right\|_{}\leqslant \frac{q}{1-q}\left\| x^{(k)}-x^{(k-1)} \right\|_{}\) - \(\displaystyle \left\| x^{(k)}-x^{*} \right\|_{}\leqslant \frac{q^{k}}{1-q}\left\| x^{(1)}-x^{(0)} \right\|_{}\)
稀疏矩阵计算
稀疏矩阵：非零元密度远小于1的矩阵. 对于很大的 \(N\)，直接法存储 \(N\)阶矩阵数据需要 \(O(N^{2})\)大小；稀疏矩阵只存储非零元位置和数据，需要 \(2N\)（或者 \(3N\)） 大小.
应用例：有限差分离散泊松方程 \[ AU=F \] \(U\) 为 \((n-1)^{2}\) 维向量，\(A\) 的带宽为 \(2(n-1)+1\)
Cholesky 分解
如果 \(A\) 是 \(n \times n\) 对称正定矩阵，则存在上三角 \(n \times n\) 矩阵 \(R\) 满足 \(A= R ^{\mathsf{T}} R\)
证 归纳 \[ A= \begin{bmatrix} a & b^{\mathsf{T}} \\ b & C \\ \end{bmatrix} \] 设 \(\displaystyle u=\frac{b}{\sqrt{a}}\)，令 \(A_1=C-uu ^{\mathsf{T}}\) 定义可逆矩阵 \[ S= \begin{bmatrix} \sqrt{a} & u^{\mathsf{T}} \\ 0 & I \end{bmatrix} \] 得到 \[ S ^{\mathsf{T}} \begin{bmatrix} 1 & 0 \\ 0 & A_1 \\ \end{bmatrix} S=A \] \(A_1\) 也是对称正定的，故 \(A_1=V ^{\mathsf{T}}V\)，最后定义 \[ R= \begin{bmatrix} \sqrt{a} & u ^{\mathsf{T}} \\ 0 & V \\ \end{bmatrix} \] 在寻找 \(R\) 时，继续对 \(A_1\) 作操作即可.
共轭梯度法
等价极值问题
定理 \(A\) 是 \(n\) 阶对称正定矩阵，\(x\) 是方程组 \(Ax=b\)的解 \(\Leftrightarrow\) \(x\) 是二次函数 \(f(x)=\frac{1}{2}x^{\mathsf{T}}Ax-b^{\mathsf{T}}x\) 的极小点. 即 \[ Ax^{*}=b \Leftrightarrow f(x^{*})= \min_{x\in \mathbb{R}^{n}}f(x) \]
证 注意 \(\nabla f(x)= Ax-b\) 且 \(f(x+\alpha y)=f(x)+\alpha (Ax-b,y)+\frac{1}{2}\alpha^{2}(Ay,y)\)
那么我们就把一个求方程根的问题转化为了一个求函数极小值的问题
最速下降法
每次找一个方向 \(d_k\)，使得 \[ x_{k+1}=x_k+\alpha_k d_k \] 其中系数 \(\alpha_k\)通过求一维问题的极小获得.
若 \(d_k=-\nabla f(x_k)=b-Ax_k\)，则该方向为最速下降方向. 有 \(d_k=r_k\)，且 \(d_k\)和 \(d_{k+1}\)正交.
计算可以得到 \[ \alpha_k=\frac{(r_k,r_k)}{(Ar_k,r_k)} \] > 可以直接考虑为什么 \(d_k\) 和 \(d_{k+1}\) 正交，注意并不一定有 \((d_i,d_j)=0, i\neq j\)
误差估计 \[ \left\| x_k-x^{*} \right\|_{A}\leqslant \left(\frac{\lambda_1-\lambda_n}{\lambda_1+\lambda_n}\right)^{k}\left\| x_0-x^{*} \right\|_{A} \]
证 用归纳法，然后用Kantorovich不等式，具体推导略. Kantorovich不等式可以见
康托洛维奇(Kantorovich)不等式的一种初等证明 https://zhuanlan.zhihu.com/p/271983329
最速下降法当矩阵条件数过大，即最大与最小特征值相差极大时收敛非常慢，不具有实际应用意义 ### 共轭梯度法
\(d_k\) 是方向， \(r_k\) 是残差（标量），\(x_k\) 是经过 \(k\) 步得到的估计
\(f(x)\) 极小问题可迭代计算，每次找方向 \(d_k\)，使 \[ x_{k+1}=x_{k}+\alpha_k d_k \] 其中系数 \(\alpha\) 通过求一维问题的极小获得
搜索方向两两共轭，且 \[ d_{k+1}=r_k+\beta_k d_k \] 其中的系数由共轭性质确定.
共轭向量 定义 \(\left\| x \right\|_{A}=(Ax,x)^{\frac{1}{2}}\) 定义 如果两个方向 \(d_i\) 和 \(d_j\) 满足 \[ (d_i,Ad_j)=0 \] 则称这两个向量关于矩阵 \(A\) 共轭. 类似可定义矩阵 \(A\) 的共轭向量组.
共轭梯度法有误差估计，收敛速度为（K为矩阵条件数） \[ \left\| x_k-x^{*} \right\|_{A}\leqslant 2\left(\frac{\sqrt{K}-1}{\sqrt{K}+1}\right)^{k} \left\| x_0-x^{*} \right\|_{A} \]
有限步收敛性定理 设 \({d_1,\cdots ,d_n}\) 为对称正定矩阵 \(A\) 的共轭向量组，令 \[ \alpha_k=\frac{(d_k,r_{k-1})}{(d_k,Ad_k)}, x_k=x_{k-1}+\alpha_k d_k, k=1,\cdots ,n \] 该方向为共轭方向法（比共轭梯度更广泛）. 对于该迭代法， \(n\) 步迭代可得准确解 \[ A x_n=b \]
证 首先 \(Ax_n=Ax_0+\alpha_1 Ad_1+\cdots +\alpha_n A d_n\)
由共轭性质，\((d_k,r_{k-1})=(d_k,b-Ax_0)\).
而 \((b-Ax_n,d_k)=(b-Ax_0,d_k)-\alpha_k(Ad_k,d_k)=0\)（代入 \(\alpha_k\) 的定义即知）
故 \(r_n\) 与所有的 \(d_k\)正交， 故 \(Ax_n=b\).
共轭梯度法 选择的搜索方向为负梯度方向和前一个搜索方向的组合，使得所得到的余项与前面的余项两两正交. 负梯度方向即为残差 \(r_k\)的方向（ \(\beta_{k-1}\)是为了使方向向量两两共轭，\(\alpha_k\) 是为了使残差向量两两共轭） \[ (r_i,r_j)=0, i \neq j \]
令 \[ r_{k-1}=b-Ax_{k-1}, \quad(r_{k-1},d_i)=0, \quad i=1,\cdots ,k-1 \] 则 \[ d_k=r_{k-1}+\beta _{k-1}d_{k-1} \] 选择合适的 \(\beta_{k-1}\) 使得与前面的搜索方向共轭 \[ (d_{k-1},Ad_k)=0 \] 可以得到 \[ \beta_{k-1}=-\frac{(d_{k-1},Ar_{k-1})}{(d_{k-1},Ad_{k-1})} \]
证明实际上 \(d_k\)与 \(d_j\)都共轭（\(j<k\)），归纳法 \[ d_{j-1} ^{\mathsf{T}}A d_k=d_{j-1} ^{\mathsf{T}} A r_{k-1}+\beta_{k-1} d_{j-1} ^{\mathsf{T}}A d_{k-1} \] 那么由归纳法 \(d_{j-1} ^{\mathsf{T}} A d_{k-1}=0\) 且 \[ A d_{j-1}=\frac{1}{\alpha_{j-1}}(r_{j-2}-r_{j-1}) \Rightarrow d_{j-1} ^{\mathsf{T}} Ar_{k-1}=\frac{1}{\alpha_{j-1}}(r_{j-2}-r_{j-1}) ^{\mathsf{T}} r_{k-1}=0 \]
共轭梯度法复杂度还是 \(O(n^{3})\) 甚至比直接法还要慢，但是在处理一些特殊问题（大型稀疏矩阵、良态即条件数很小）时有非常高的效率.
总结算法如下 初始：\(r_0=b-Ax_0, \quad d_1=r_0\) For \(k=1,\cdots ,n\) \(\displaystyle \alpha_k=\frac{(r_{k-1},r_{k-1})}{(d_k,Ad_k)}\) \(x_k=x_{k-1}+\alpha_k d_k\) \(r_k=r_{k-1}-\alpha_k A d_k\) （新残差） \(\displaystyle \beta_k=\frac{(r_k,r_k)}{(r_{k-1},r_{k-1})}=\frac{-(d_k, A r_k)}{(d_k,A d_k)}\) \(d _{k+1}=r_k+\beta_k d_k\) End For
（下标可能与书本有所不同）
预条件 若 \(A\)为病态矩阵，过程易受舍入误差的影响，使用预条件可以在 \(\sqrt{n}\) 步得到收敛解.
目标是选择合适的矩阵 \(M\)，使得方程 \[ M ^{-1} Ax= M ^{-1}b \] 更“好”. \(M\) 当然可逆，称为预条件子. 我们希望 - \(M\)和 \(A\) 足够接近 - \(M\) 容易求逆
雅可比预条件子 \(M=D\) ，其中 \(D\) 是 \(A\) 的对角线矩阵. 用 \((v,w)_{M}\) 替换欧几里得内积. 则 \(M ^{-1}A\) 相对于新的内积仍然是对称正定矩阵.
高斯-赛德尔预条件子 对称连续过松弛（SSOR），中预条件子定义为 \[ M=(D+\omega L)D ^{-1}(D+\omega U) \] 其中 \(A=L+D+U\) 被分为下三角部分、对角线以及上三角部分. \(\omega \in [0,2]\)，当 \(\omega=1\) 时称为高斯-赛德尔预条件子
经过一些化简可以得到有预条件的算法 \(x_0=\) 初始条件 初始：\(r_0=b-Ax_0, \quad d_1=z_0=M ^{-1}r_0\) For \(k=1,\cdots ,n\) \(\displaystyle \alpha_k=\frac{(r_{k-1},z_{k-1})}{(d_k,Ad_k)}\) \(x_k=x_{k-1}+\alpha_k d_k\) \(r_k=r_{k-1}-\alpha_k A d_k\) （新残差） \(z_k=M ^{-1}r_k\) \(\displaystyle \beta_k=\frac{(r_k,z_k)}{(r_{k-1},z_{k-1})}\) \(d _{k+1}=z_k+\beta_k d_k\) End For
在SSOR下，不需要求 \(M ^{-1}\)，在计算 \(z_{k}= M ^{-1} r_k\) 时，由于 \(M=(I+\omega LD ^{-1})(D+\omega U)\) 是下上三角矩阵的乘积，故该方程可以通过两次回代求解. \[ (I+\omega LD ^{-1})c=v \\ (D+ \omega U)z=c \\ \]
可以将具备SSOR的共轭梯度法与一般的迭代法作比较.