数值分析笔记

这里是数值分析的部分笔记，(对前三章的算法都用python做了实验，程序见具体内容处)
PS：其中比较不重要的地方偷懒用了一点点numpy的API。内容写的比较简略，仅供参考，望见谅（更详细的数值分析Java实现 请见清华相关实验网站）

Chap 1:误差

$a=\underline{+}{a}_0{a}_1{a}_2...a_m.a_{m+1}...a_n$

def:

绝对误差(界)

相对误差(界)

有效数字

calculate：

误差界:注意近似时向上取整

有效数字的计算(以舍入法为基准)：

n+1-s :s是第一位非零数字角标，n的算法:$|a-\bar{a}|\le \frac{1}{2}\times 10^{-(n-m)}$

（1）真实值的有效数字

（2）已知实际值和近似值，求近似值的有效数组

thm(1.2.1):a的近似数$\bar{a}$有效数字有n+1-s位，相对误差有估计式:

$|\frac{a-\bar{a}}{\bar{a}}|\le \frac{1}{2a_s}\times 10^{-(n-s)}{a}_s\ne 0$是$\bar{a}$的第一位数字

浮点数

def:

p进制表示：$x=\underline{+}0.d_1{d}_2...d_t\times p^J,0\le d_i\le p-1$

尾数:$0.d_1{d}_2...d_t$

计算机表示中为$\underline{+}{d}_1{d}_2...d_{7}J$(阶码)

上溢下溢

舍入法&断尾法

calculate

（1）相加相减:先对位再加减

（2）避免减法:分子有理化；三角变换；Taylor展开

（3）乘除:双精度运算

（4）避免大数小数运算(不满足结合律)

def:向前误差分析；向后误差分析

Chap 2:解方程组

解线性方程组的Gauss消去法

操作:(1)消元(消为上三角矩阵)

​ (2)回代

原理:右乘矩阵行变换

优化方法1:列主元Gauss消去法

原理:避免在较小值为主元时，后面的行减第一行乘一个因数(乘子)，的乘数过大引起误差。(列中元素差距大)

操作:在一般Gauss消去法中对每列消元时把该列绝对值最多的数所在的行换到上面进行消元。

优化方法2：按比例列主元消去法

原理:为避免每一行中元素大小差距过大带来误差

操作:（1）在消第i列时取系数矩阵中$\frac{|a_{ji}|}{|第j行元素最大的绝对值|}$对应的$a_{ij}$为主元，(即绝对值最接近所在行最大绝对值的 元素)。

​ （2）不用交换行，但是改变回代顺序

remark:上面两种方法都是在规避用较小的数求比例得到乘子乘到较大的数上带来误差

优化方法3：Gauss-Jordan 消去法

原理:为避免回代过程，直接将矩阵化为对角矩阵。

操作:仍然采取列主元的方法，但是在消元过程中，消第k列时，将1,2,…,k-1行第k列也消去。

#Gauss-Jordan消去法的python实现(因为电脑没有装java和C++的编译器orz)
import numpy as np
Abb=np.array([(1,2,3,4),(2,12,3,5),(2,10,2,4)])
Abb = Abb.astype(np.float64) #转化为浮点型(否则运算时会自动化为int)

def GaussJordan(Ab,n): 
    row_rec=[];#记录已经有主元的行
    for k in range(n):
        max_k=0
        for i in range(n):
            if ((i in row_rec)==False)&(abs(Ab[i,k])>max_k):#从剩下行中寻找主元
                max_k=abs(Ab[i,k])
                index=i#记录主元所在行数
        row_rec.append(index)
        for j in range(n):
            if (j!=row_rec[-1]):#消去剩余行
                m=-(Ab[j,k]/Ab[row_rec[-1],k])#乘子
                for i in range(n-k+1):
                    Ab[j,i+k]+=m*Ab[row_rec[-1],i+k]
    #计算结果
    res=[]
    for i in range(n):
        res.append(Ab[row_rec[i],n]/Ab[row_rec[i],i])
    return(res)

print(GaussJordan(Abb,3))
res:[1.0, -0.0, 1.0]

直接三角分解法

1.LDR分解

$A=LDR$,L是单位上三角矩阵，D是对角矩阵，R是单位下三角矩阵

thm:矩阵A存在唯一LDR分解$\iff$A的顺序主子矩阵均非奇异

$A=LDR=LU$

从而直接解方程:

$\{\begin{array}{c}Ly=b\\ y=Ux\end{array}$

2.Doolittle分解

def :A=LU,L是单位下三角矩阵，U是上三角矩阵

原理:Gauss消元法的矩阵表示

表示消去第2列的右乘矩阵如下

$L_2=\left[\begin{array}{cccc}1& 0& 0& ...\\ 0& 1& 0& 0\\ 0& -l_{3,2}& ...& \\ 0& -l_{n,2}& 0...& 1\end{array}\right]$

记${\vec{l}}_k=\left[\begin{array}{cc}& 0\\ & ...\\ & l_{k+1,k}\\ & ...\\ & l_{m,k}\end{array}\right]$

从而$L_k=I-\overrightarrow{l_k}{e}_k,e_k$是第k个坐标是1，其余分量是0的n维向量

$(I-\overrightarrow{l_k}{e}_k)(I+\overrightarrow{l_k}{e}_k)=I\to L_k^{-1}=I+\overrightarrow{l_k}{e}_k$

$i<j$时有$L_i^{-1}{L}_j^{-1}=I+\overrightarrow{l_i}{e}_i+\overrightarrow{l_j}{e}_j$

$a_{ii}^{(i-1)}=0$(即主元为0)时要交换行

此时非下三角矩阵

要是每一个主元都不为0的充要条件是:A的每一个顺序主子矩阵非奇异(这一点根据Gauss消去法的操作过程是显然的)

3.Crout分解

原理:类比Gauss消元法的矩阵表示，按列消元从左向右

原理:$A=LU,U$是单位上三角矩阵，$i<j,a_{ij}={\sum }_{k=1}^i{a}_{i,k}{u}_{k,j}$

先求L中第k列，再求U中第k行

操作:将$L,U-I$存入同一个矩阵。对矩阵A进行操作:

​ (1)第一列不变。

​ (2)第一行除对角线以外除以$a_{11}$

​ 一行一列边框over

​ (3)第二列中元素从减去所在行和所在列边框互乘元素

​ (4)第二列减去互乘元素除以所在对角线元素

​ …

Crout分解存在唯一的充要条件依旧是顺序主子式非奇异

优化方法1:紧凑的crout分解

在求第i列之后可以直接求得$Ly=b$中的y

#Abbb顺序主子式非奇异
Abbb=np.array([(3,2,1,2),(2,4,1,-1),(1,2,4,3)])
Abbb = Abbb.astype(np.float64) 

def Crout_jin(Ab,n):
    y=[Ab[0,-1]/Ab[0,0]]#求L(1,1)
    for i in range(n-1):Ab[0,i+1]=Ab[0,i+1]/Ab[0,0]#求U中第1行
    for i in range(n-1):
        for j in range(n-i-1):#求U第i+2行
            minus=0
            for t in range(i+1):
                minus+=Ab[j+i+1,t]*Ab[t,i+1]
            Ab[j+i+1,i+1]-=minus##如果是按列选主元的话在这里进行行变换
        for j in range(n-i-2):#求L第i   +2列
            minus=0
            for t in range(i+1):
                minus+=Ab[t,i+2+j]*Ab[i+1,t]
            Ab[i+1,i+2+j]=(Ab[i+1,i+2+j]-minus)/Ab[i+1,i+1]
        y_i=Ab[i+1,-1]
        for j in range(i+1):#求y   紧凑的crout分解在于在这个循环里计算y的第i+2个分量
            y_i-=y[j]*Ab[i+1,j]
        y_i=y_i/Ab[i+1,i+1]
        y.append(y_i)
    x=[]#求解x,从最后一个分量向前求解
    for j in range(n):
        minus=0
        for i in range(j):
            minus+=Ab[-(j+1),-(i+2)]*x[i]
        x_j=y[-(j+1)]-minus;
        x.append(x_j)
    return(x[::-1])#最后reverse

Crout_jin(Abbb,3)
res:[0.9999999999999998, -0.9999999999999998, 1.0]

优化方法2:按列选主元的Crout方法

和doolittle分解一样在求完一列之后将绝对值最大的所在行提到上面来

remark:考试时不要求按照书上步骤，只需要写出最终形式即可

4.Cholesky 分解（对实对称矩阵）

thm:n阶实对称矩阵A是正定的$\iff$$\exists$n阶非奇异下三角矩阵$L,s.t.:$
$$A=L{L}^T$$
且L的主对角元均为正数时该分解唯一

pf:$\Rightarrow$A是正定的，则顺序主子式非奇异有$A=L_{1}DR,A^T=R^T{D}^T{R}^T$，且$\forall i,D_{ii}>0$由分解式的唯一性得$L_1=R^T,A=L_1{\sqrt{D}}_1(L_1{\sqrt{D}}_2{)}^T$,显然L主对角元为正数时(${\sqrt{D}}_1={\sqrt{D}}_2$)分解唯一。

$\Leftarrow$$xA{x}^T=xL(xL{)}^T\ge 0$,由于$xL\ne 0$等号不能取到，从而A正定。

平方根法

如果直接将$A=L{L}^T$利用Crout算法的步骤计算则对角线元素$l_{ii}$需要开方求解，由于开方计算复杂度大，且会带来更多误差，则考虑下面的分解。

5.$LD{L}^T$分解

$A=LD{L}^T$其中L是单位下对角矩阵，D是对角矩阵。则根据$L{L}^T$分解的唯一性，可知A是对称正定时$LD{L}^T$也是存在唯一的。

对此分解同样利用Crout分解的方法分成$LD$和$L^T$，再转换成$L,D{L}^T$求解但是在每次计算第i行第i列的时候加上$d_{ii}$,来避免开方运算。也就是

改进的平方根算法

计算过程中
$$A=LD{L}^T,Ax=b\to \{\begin{array}{c}Ly=b\\ y=D{L}^{T}x\end{array}$$
计算上和一般的Crout算法有以下不同:

1. 求完L第i列的第k行($n-i\le k\le n$)个元素之后，第i行的第k个元素$D{L}_{ik}^T=\frac{L_{ik}}{D{L}_{ii}^T}$
2. 并非先算DL第i列的后n-i个元素,再算$L^T$的第i行的后n-i个元素，而是依次算DL和$L^T$的第i个顺序主子式
3. 且不需要同时储存$DL和L^T$,只需要储存$L$和D的对角线即可

6.解三对角矩阵方程组的三对角算法

$$\left[\begin{array}{cccc}{d}_1& c_1& & \\ a_2& d_2& c_2& \\ & ...& ...& ...\\ & a_{n-1}& d_{n-1}& c_{n-1}\\ & & a_n& d_n\end{array}\right]=\left[\begin{array}{cccc}{p}_1& & & \\ a_2& p_2& & \\ & ...& ...& ...\\ & a_{n-1}& p_{n-1}& \\ & & a_n& p_n\end{array}\right]\left[\begin{array}{cccc}1& q_1& & \\ & 1& q_2& \\ & ...& ...& ...\\ ...& & 1& q_{n-1}\\ & & & 1\end{array}\right]$$

对半带宽稀疏矩阵，往往有更简单的分解方法。这里对于对角占优的三对角矩阵:

1. $|d_1|>|c_1|>0$

2. $|d_k|\ge |a_k|+|c_k|$且$a_k{c}_k\ne 0,k=2,3,...,n-1$

3. $|d_n|>|a_n|>0$

   在上述三个条件下$p_1,p_2,...,p_n$皆非0。

   pf:先归纳法证明$\forall i,|q_i|<1$:由$(1)$知$|q_1|<1$

   $a_i{q}_{i-1}+p_i=d_i,p_i{q}_i=c_i$由此得到$|q_i|$的递推进行证明。

   同样利用上面两个等式得到$|p_i|>0$

对系数矩阵为上述矩阵的方程的解法就是Crout算法，但是计算过程中隐含$p_i,q_i$的交错递归所以被称为追赶法

Appendix

1.QR分解

$A=QR$,Q为正交矩阵，R为上三角矩阵

分解过程:Gram-Schmidt过程即列变换反复生成正交向量最后得到标准正交集

2.矩阵对角化

$S=[a_1,...,a_n],a_i$是特征向量，$AS=ST,T=(t_{ij})=({\lambda }_i{)}_{ii}\to A=ST{S}^{-1}$

要求:S有n个独立的特征向量

要求:$T\in C^{n\times n}:$正规矩阵,$T\in R^{n\times n}:$自伴矩阵

例子:实对称矩阵:特征值相互正交(易证)

3.Schur分解

$A\in C^{n\times n},U^{H}AU=T$,T是上三角矩阵，特征值在T的对角线上。

再进一步，进行Jordan分解

pf:1.需要明确一点:对上三角矩阵A，$[x_1,...x_n]A$使得$x=[x_1,...,x_n]$中任意分量张成的空间在变化前后不变。

​ 即$\forall j=1,2,....,n,span(v_1,...,v_j)$均在T下不变

2.对任意矩阵都可以表示为某个基上的上三角矩阵(构造过程由归纳0法，通过张成特征向量张成的空间以外的空	  间$(T-\lambda I)V$进行降维，详见done right P114)

​ 3.表示为某个基上的上三角矩阵则可以表示为某个规范正交基上的上三角矩阵(根据Gram-Schmit过程由1.中的 陈述或者用QR分解均可说明)

4.极分解

$T\in \mathcal{L}(V)$,存在等距同构(Hermit矩阵)S，使得$T=S\sqrt{T^{\ast }T}$

5.满秩分解

$rank(A)=r,A=BC,B\in C^{s\times r},C\in C^{r\times n},B,C$秩都为r

$A=BC$矩阵的秩为r，线性无关的r行主子式为B，C对应的主子式是单位矩阵

具体例子见

矩阵论有关知识

def:矩阵相似$A=SB{S}^{-1},S$可逆(即同一线性变换在不同基上的矩阵)

def:矩阵合同$A=S^{T}BS$(同一双线性映射在不同基上的矩阵)

def:伴随算子(adjoin):对T，T的伴随算子$T^{\ast },s.t.<Tv,w>=<v,T^{\ast }w>$,可知$T^{\ast }$是T的共轭转置。

def:自伴算子(adjoint):$T=T^{\ast }$

def:正规算子(normal):$T{T}^{\ast }=T^{\ast }T\iff ||Tx||=||T^{\ast }X||$

def:对任意矩阵$A,A=P^{-1}JP,J$是jodan矩阵，P可逆

prop:正规算子对应不同特征值的特征向量是正交的:

(pf:$(a-b)<u,v>=a<u,v>-b<u,v>=<Tu,v>-<u,T^{\ast }v>=0$)

prop:$rank(A{A}^H)=rank(A)$

prop:$A,A^T$有相同的特征值

prop:可对角化矩阵迹等于特征值之和

prop:AB,BA有相同特征值

thm:复谱定理:正规矩阵能做对角分解(S为正交矩阵)

thm:实谱定理:在实数域上对角分解要求自伴矩阵(S为Hermite矩阵)

6.广义矩阵

def:$A\in C^{s\times n},ifG\in C^{n\times s}$,满足下面M-P方程:则称G是A的广义逆矩阵:

• AGA=A

• GAG=G

• $(AG{)}^H=AG$

• $(GA{)}^H=GA$

prop：零矩阵的广义逆是零矩阵

thm:设$A\in C^{s\times n}$,则A的广义逆是存在唯一的

pf:设$G_1,G_2$是A的广义逆则:$G_1=G_{1}A{G}_1=G_1(A{G}_1{)}^H=G_1{G}_1^H{A}^H=G_1{G}_1^H(A{G}_{2}A{)}^H$

​ $=G_1{G}_1^H{A}^H(A{G}_2{)}^H=G_1{G}_1^H{A}^{H}A{G}_2=G_1(A{G}_1{)}^{H}A{G}_2$

​ $=G_{1}A{G}_{1}A{G}_2=G_{1}A{G}_2$

​ 类似地$G_2=G_{1}A{G}_2$，唯一性得证

​ 利用矩阵的满秩分解$Y(A)=r,A=BC,B\in C^{s\times r},C\in C^{r\times n},B,C$秩都为r

​ 构造广义逆矩阵:$G=C^H(C{C}^H{)}^{-1}(B^{H}B{)}^{-1}{B}^H$

remark:A的广义矩阵记为$A^{+}$

remark:$(AB{)}^{+}$与$B^{+}{A}^{+}$一般不相等

prop:

• $(A^{+}{)}^{+}=A$

• $(A^H{)}^{+}=(A^{+}{)}^H$

• $(A^T{)}^{+}=(A^{+}{)}^H$

• $k\in R,(kA{)}^{+}=k^{+}{A}^{+}$

• $A^H=A^{h}A{A}^{+}=A^{+}A{A}^H$

• $(A{A}^H{)}^{+}=A^{+}(A^H{)}^{+};(A{A}^H{)}^{+}=(A^H{)}^{+}{A}^{+}$

• $A^{+}=(A^{H}A)A^H=A^H(A{A}^H{)}^{+}$

• $U,V$是酉矩阵，则$(UAV{)}^{+}=V^H{A}^{+}{U}^H$

• $A^{+}AB=A^{+}AC\iff AB=AC$

example:

• $A$是Hermite矩阵，$A^{+}$也是Hermite矩阵
• A是正规矩阵，$(A^2{)}^{+}=(A^{+}{)}^2$

apply:

线性方程组$Ax=b$无解时，找到x使得$||A\hat{x}-b|{|}_2$最小。$\eta$是最小二乘解$\iff \eta 是A^{H}Ax=A^{H}b$的解

attention:$x\perp span(a_1,...,a_n),A=(a_1,...,a_n),\to \forall y:(x,Ay)=0\iff (A^{H}x,y)=0\iff A^{H}x=0$

通解$x=A^{+}b+(I-A^{+}A)y$，其中$A^{+}b$是极小最小二乘解

行列式的计算

（1）Gauss消去法

对于顺序的Gauss消去法$det(A)=a_{11}{a}_{22}^{(1)}...a_{nn}^{(n-1)}$即各主元的乘积

对于中间进行行变换的Gauss消去法则需要乘以$(-1{)}^s$,s为消元中行变换的次数

（2）Crout分解

$A=LU,det(A)=det(L)$

（3）$L{L}^T,LD{L}^T$消去法

$det(A)=$对角线元素的乘积

逆矩阵的计算

对于非奇异n阶矩阵A，用Gauss-Jordan消去法对A进行消元时$M_n...M_{1}A=I$，其中
$$M_k=\left[\begin{array}{ccccc}1& & m_{1k}& & \\ & 1& ...& & \\ & & m_{kk}& ...& \\ & & ...& 1& \\ & & m_{nk}& ...& 1\end{array}\right]$$
从而$A^{-1}=M_n...M_1$

对于$M_k$的求和:$M_n...M_2{M}_1$存在特殊解法:

$T_1=M_1(e_{11}+A(e_{22}+…+e_{nn})) $ $e_{ii}{e}_{jj}=0,i\ne j,e_{ii}{e}_{ii}=e_{ii}$

$T_2=M_2(e_{22}+T_1(e_{11}+e_{22}+...e_{nn}))$

$A^{-1}=M_n...M_1=T_n$

范数

下面讨论的都是有限维的情况，在无限维的情况下，max=sup,min=inf其中多数叙述仍然成立。

（一）向量的范数

线性赋范空间

def:(1)非负性(正性&定性)

​ (2)齐次性(线性性)

​ (3)三角不等式

remark:$l_p$范数的三角不等式证明由Minkowski不等式可得

example:

1.$\infty$范数 def:$||x|{|}_{\infty }=li{m}_{p\to \infty }({\sum }_{i=1}^{\infty }|x_i{|}^p{)}^{\frac{1}{p}}$

​ $||x|{|}_{\infty }=ma{x}_{1\le i\le n}|x_i|$

​ pf: 记$max|x_i|=|x_j|,|x_j|\le (\sum |x_i{|}^p{)}^{\frac{1}{p}}\le n^{\frac{1}{p}}|x_j|,n\to \infty ,||x|{|}_{\infty }=|x_j|$

​

2.$||x|{|}_2$:Euclid范数，相关的有Cauchy-Schwarz不等式:$|(x,y)|\le ||x|{|}_2||y|{|}_2$(实质上是Euclid内积空间的性质)

3.A是给定的n阶实对称正定矩阵，$x\in R^n$,则$||x|{|}_A=(x^{T}Ax{)}^{\frac{1}{2}}$是$R^n$中的一种范数

pf:对于三角不等式的证明:$||x+y|{|}_A=||L^T(x+y)|{|}_2=||L^{T}x+L^{T}y|{|}_2\le ||L^{T}x|{|}_2+||L^{T}y|{|}_2=||x|{|}_A+||y|{|}_A$

prop:

1.对任意向量范数，$\forall x,y\in R^n$,恒有:$|||x||-||y||\le ||x-y||$,(在一维情况下即为绝对值不等式)​

pf:$||x||\ge ||y||$时,$|||x||-||y|||=||x||-||y||\le ||x-y||$

​ $||x||\le ||y||$时，同理。

2.$R^n$空间中，任意向量范数$||x|{|}_{\alpha }$关于$l_2$范数是x的一致连续函数 (以欧氏空间上的拓扑为基础拓扑)

pf:$\forall x=(x_1,...,x_n),y=(y_1,...,y_n),||x-y|{|}_2=({\sum }_{i=1}^n(x_i-y_i{)}^2{)}^{\frac{1}{2}}<\delta$

​ $|||x|{|}_{\alpha }-||y|{|}_{\alpha }|\le ||x-y|{|}_{\alpha }$

​ $||x|{|}_{\alpha }=||\sum x_i{e}_i|{|}_{\alpha }\le \sum |x_i|||e_i|{|}_{\alpha }\le ||x|{|}_2(\sum ||e_i|{|}_{\alpha }^2{)}^{\frac{1}{2}}$(第一个不等式由三角不等式可得)

​ $({\sum }_{i=1}^n||e_i|{|}_{\alpha }^2{)}^{\frac{1}{2}}=const$

​ 从而 $|||x|{|}_{\alpha }-||y|{|}_{\alpha }|\le ||x-y|{|}_{\alpha }\le ||x-y|{|}_2(\sum ||e_i|{|}_{\alpha }^2{)}^{\frac{1}{2}}<M\delta$

3.$R^n$空间中的一切向量范数都是等价的，即对于任意两种范数$||x|{|}_{\alpha },||\beta |{|}_{\alpha }$总存在与x无关的正常数$c_1,c_2\in R,s.t.$(利用赋范线性空间特点)
$$c_1||x|{|}_{\beta }\le ||x|{|}_{\alpha }\le c_2||x|{|}_{\beta },\forall x\in R^n$$
pf:对于这点的证明,我们先将x单位化、标准化再做其它考虑。

​ $||x|{|}_{\alpha }=||x|{|}_2||\frac{x}{||x|{|}_2}|{|}_{\alpha }=||x|{|}_2||y|{|}_{\alpha }$其中 y ∈ { y : ∣ ∣ y ∣ ∣ 2 = 1 , y ∈ R n } , y\in \lbrace y:||y||_2=1 ,y\in R^n\rbrace, y∈{y:∣∣y∣∣2​=1,y∈Rn},又$t(y)=||y|{|}_{\alpha }$是连续函数从而在闭球 面上有最大值和最小值$c_1,c_2$，从而$c_2||x|{|}_2\le ||x|{|}_{\alpha }\le c_1||x|{|}_2$利用$||x|{|}_2$做过渡可轻松进行证明。

4.$R^n$空间中的$l_p(p=1,2,\infty )$范数满足以下关系式:

• $||x|{|}_{\infty }\le ||x|{|}_1\le n||x|{|}_{\infty }$
• $||x|{|}_{\infty }\le ||x|{|}_2\le \sqrt{n}||x|{|}_{\infty }$
• $\frac{1}{\sqrt{n}}||x|{|}_1\le ||x|{|}_2\le ||x|{|}_1$ 由Cauchy-Schwarz不等式可得

remark:这一步的证明过程适用于证明上一步，3强于2

example:设$A\in R^{n\times n}$,在$R^n$中给定范数$||\cdot |{|}_{\alpha }$,证明$g(x)=||Ax|{|}_{\alpha }$关于范数$||\cdot |{|}_{\alpha }$是x的连续函数。

​ remark:从线性泛函的角度上看上述叙述是显然的。

​ 证明: ∣ g ( x ) − g ( y ) ∣ ≤ ∣ ∣ A ( x − y ) ∣ ∣ α ≤ c 1 ∣ b ⃗ ⋅ ( x − y ⃗ ) ∣ ≤ c 1 m a x { ∣ b i ∣ } ∣ ∣ x − y ∣ ∣ 1 ≤ c 1 m a x { ∣ b i ∣ } c 2 ∣ ∣ x − y ∣ ∣ α |g(x)-g(y)|\leq ||A(x-y)||_{\alpha}\leq c_1 |\vec{b}\cdot(\vec{x-y})|\leq c_1max\lbrace|b_i|\rbrace ||x-y||_1\leq c_1max\lbrace|b_i|\rbrace c_2||x-y||_{\alpha} ∣g(x)−g(y)∣≤∣∣A(x−y)∣∣α​≤c1​∣b ⋅(x−y ​)∣≤c1​max{∣bi​∣}∣∣x−y∣∣1​≤c1​max{∣bi​∣}c2​∣∣x−y∣∣α​ 其中 m a x { ∣ b i ∣ } max\lbrace|b_i|\rbrace max{∣bi​∣}不超过A中每列绝对值最大的数的绝对值的求和。

(二)、矩阵的范数

首先对方阵进行定义:

$R^{n\times n}$表示全体$n\times n$阶实矩阵构成的线性空间，设$A\in R^{n\times n}$若在$R^{n\times n}$中定义实值函数$||A||$,满足:

1. 正定性
2. 线性性
3. 三角不等式
4. 相容性 $||AB||\le ||A||\cdot ||B||,\forall A,B\in R^{n\times n}$

称$||A||$为矩阵A的一种范数或者模或矩阵范数，相容范数。

def:矩阵范数$||\cdot |{|}_{\beta }$和向量范数$||\cdot |{|}_{\alpha }$是相容的:$||Ax|{|}_{\alpha }\le ||A|{|}_{\beta }||x|{|}_{\alpha },\forall A\in R^{n\times n},x\in R^n$

thm:设$||\cdot |{|}_{\beta }$是$R^{n\times n}$中的任意一种矩阵范数，则在$R^n$中至少存在一种向量范数$||\cdot |{|}_{\alpha }$使得$||\cdot |{|}_{\alpha },||\cdot |{|}_{\beta }$是相容 的。

   pf：构造$||\cdot||_{\beta}=||[x,0,...,0]||_\beta,\qquad [x,0,...,0]\in R^{n\times n}$

def:假设矩阵范数$||A|{|}_{\beta }$和向量范数$||x|{|}_{\alpha }$相容且$\forall A\in R^{n\times n}$都存在一个非零向量$x_0\in R^n$,使得
$$||A{x}_0|{|}_{\alpha }=||A|{|}_{\beta }||x_0|{|}_{\alpha }$$
则说$||A|{|}_{\beta }$是从属于向量范数$||x|{|}_{\alpha }$的矩阵范数

remark:从属的必要条件是$||I|{|}_{\beta }=1$

thm:对于$R^n$中的向量范数$||\cdot |{|}_{\alpha }$都存在从属于它的矩阵范数$||\cdot |{|}_{\beta }$

​ 其中一种为$||A||=ma{x}_{||x|{|}_{\alpha }=1}||Ax|{|}_{\alpha }$，次范数被定义为算子范数。

pf:由范数等价性&$||Ax||$是有界闭集上的连续函数知最大值能取到。(1)非负性(正且定)(2)齐次性(3)三角不等式(3)相容性 $||AB||=max||ABx||\le max||A||||Bx||\le max||A||||B||||x||=||A||||B||$；证明相容性$||A|{|}_{\alpha }\ge ||\frac{Ax}{||x|{|}_{\alpha }}|{|}_{\alpha }$

prop:

同样有$|||A||-||B|||\le ||A-B||$

example:

• 1-范数 $||A|{|}_1=ma{x}_{1\le j\le n}{\sum }_{i=1}^n|a_{ij}|$

• 2-范数(谱范数) $||A|{|}_2={\rho }^{\frac{1}{2}}(A^{T}A),\rho$表示矩阵谱半径(即特征值最大值)

  注意这里是$A^{T}A$本质上$AB,BA$具有相同的特征值，但是其特征向量不一定相同；同样P和$P^T$也有相同的特征值但是特征向量不一定相同(由特征多项式的形式可看出)

  pf:$ABx=\lambda x,BABx=\lambda Bx,Bx$是$BA$关于$\lambda$的特征向量。

• $\infty$-范数 $||A|{|}_{\infty }=ma{x}_{1\le i\le n}{\sum }_{i=j}^n|a_{ij}|$

• F-范数 $||A|{|}_F=({\sum }_i{\sum }_j|a_{ij}{|}^2{)}^{\frac{1}{2}}||Ax|{|}_2\le ||A|{|}_F||x|{|}_2，由于||I|{|}_F\ne 1$知F-范 数不从属于任何范数

$||A||_2\leq ||A||_F\leq \sqrt{n}||A||_2 $

pf:$||A||_2^2=max_i \lambda_i\geq \frac{1}{n}(\lambda_1+…+\lambda_n )=\frac{1}{n}tr(A^{TA)=\frac{1}{n}||A||_F}2 $实对称矩阵的迹等于特征值的和，且是半正定矩阵（特征值不小于0）。

上述定义在共轭转置的基础上可推广到复数域上的矩阵，在考虑$C^n,C^m$上的范数后可推广到$C^{n\times n}$非方阵上，详见P69.

thm:

对于任意矩阵范数$||\cdot ||$,都有$\rho (A)\le ||A||,\forall A\in C^{n\times n}$

pf:由于对于所有矩阵范数都存在向量范数$||\cdot |{|}_{\alpha }$使得两者相容，即$||Ax|{|}_{\alpha }\le ||A||\cdot ||x|{|}_{\alpha }$,取$x$为A的特征向量，则有$\lambda ||x|{|}_{\alpha }\le ||A||\cdot ||x|{|}_{\alpha }$,QED.

补充Jordan矩阵幂运算

thm:设$A\in C^{n\times n},\forall ϵ>0,\exists ||\cdot |{|}_{\beta },s.t.||A|{|}_{\beta }\le \rho (A)+ϵ$

pf:A必与一个Jordan标准型$J$相似，即存在非奇异矩阵P使得$P^{-1}AP=J，D=diag(1,ϵ,...,{ϵ}^{n-1}),D^{-1}JD=\tilde{J}$将J的每一个非对角1换成$ϵ,\tilde{J}=Q^{-1}AQ,Q=PD,||Q^{-1}AQ|{|}_1=||\tilde{J}|{|}_1\le \rho (A)+ϵ||A|{|}_{\beta }=||Q^{-1}AQ|{|}_1$是一种矩阵范数 QED.

remark:上面两个定理说明$\rho (A)$是A的范数的下确界。

thm:Banach Lemma

$B\in C^{n\times n},\rho (B)<1$,则矩阵$I\pm B$都是非奇异的，而且对任何$||I||=1$的矩阵范数$||\cdot ||$，若有$||B||<1$,则
$$||(I\pm B^{-1})||\le \frac{1}{1-||B||}$$
pf:因为$\rho (B)<1,I\pm B$的特征值的模非0(做jordan分解)。

​ 又$(I+B)(I+B{)}^{-1}=I,(I+B{)}^{-1}=I-B(I+B{)}^{-1}$
$$||(I+B{)}^{-1}||\le ||I||+||B||||(I+B{)}^{-1}||\to (1-||B||)||(I+B{)}^{-1}||\le ||1||$$
QED

(三)、向量列与矩阵列的收敛

def:向量列的收敛:各分量收敛

def:矩阵列的收敛:各元素收敛

矩阵级数收敛:部分和收敛

thm:$C^{n\times n}$空间中的矩阵序列$A_1,...$收敛于A的充要条件是存在一种矩阵范数$||\cdot ||,s.t.li{m}_{k\to \infty }||A_k-A||=0$

pf:(1)对于向量范数有，从而进行类比$lim |B|=lim |C|\leftarrow lim B=limC $当且仅当$lim |C|=0$从左到右可推，从而考虑$lim|A_k-A|\rightarrow 0$

​ (2)证明有限维空间上任何范数等价。所以只对$|A{|}_M=nmax|a_{ij}|$,$lim||A_K-A|{|}_M=limnmax||a_{ij}-a_j||$这样各元素收敛等价于M范数收敛。

thm:$C^n$空间中的数列$a_1,...$收敛于a的充要条件是存在一种矩阵范数$||\cdot ||,s.t.li{m}_{k\to \infty }||a_k-a||=0$

thm:设$A\in C^{n\times n},li{m}_{k\to \infty }{A}^k=0\iff \rho (A)<1$

pf:$\Rightarrow$反设$\rho (A)\ge 1,\rho (A^k)\ge 1\to ||A^K||\ge 1$,从而$lim||A^k||\ne 0\to A\ne 0$

$\Leftarrow$ $\rho (A)<1,{ϵ}_0=\frac{1-\rho (A)}{2}>0$存在矩阵范数

$||\cdot |{|}_{\beta },s.t.||A|{|}_{\beta }\le \rho (A)+{ϵ}_0<1\to ||A^k|{|}_{\beta }\le ||A|{|}_{\beta }^k\le q^k,q<1,||A^k||\to 0,A^k\to 0$

thm:$A\in C^{n\times n},{\sum }_{k=0}^{\infty }{A}^k$收敛$\Leftrightarrow $$\rho (A)<1$, 且若$\rho (A)<1$,则$(I-A{)}^{-1}$存在 , ${\sum }_{k=0}^{\infty }{A}^k=(1-A{)}^{-1}$

pf:$\Rightarrow$令$S_n={\sum }_{k=0}^m{A}^k$,则由矩阵幂级数收敛知$li{m}_{m\to \infty }{S}_m=S$,$lim{A}^k=lim(S_k-S_{k-1})=0$由上述定理知$\rho (A)<1$

$\Leftarrow$设$\rho (A)<1$,$(I-A{)}^{-1}$存在

$(I-A)S_k=(I-A)(1+A+A^2+...+A^k)=I-A^{k+1}$

$S_k=(I-A{)}^{-1}(I-A^{k+1})$

$=(I-A{)}^{-1}-(I-A{)}^{-1}{A}^{k+1}$

从而$||S_k-(I-A{)}^{-1}||\le ||(I-A{)}^{-1}||||A^{k+1}||$由$\rho (A)<1\to li{m}_{k\to \infty }||A^{k+1}||=0$从而$S_n\to (I-A{)}^{-1}$

thm:设$A\in C^{n\times n}$,若对$C^{n\times n}$中的某种范数$||\cdot ||$有$||A||<1$,则有
$$||(1-A{)}^{-1}-(I+A^2+...+A^k)||\le \frac{||A|{|}^{k+1}}{1-||A||}$$
pf：$LHS=||A^{k+1}+A^{k+2}+...||\le ||A^{k+1}||||1+A+A^2+...||\le ||A^{k+1}||||(1-A{)}^{-1}||\stackrel{BanachLemma}{\le }\frac{||A^{k+1}||}{1-||A||}=RHS$

(banach引理需要加上条件$||I||=1$)

如果没有$||I||=1:||(I-A{)}^{-1}-(I+A+A^2+...+A^k)||$

​ $\le ||A|{|}^{k+1}+...+||A|{|}^{k+m}+||A|{|}^{k+m+1}(I-A{)}^{-1},(m\to \infty )$

​ $\le {\sum }_{i=k+1}^{\infty }||A|{|}^i=RHS$

thm:酉矩阵的保范性质

$A\in C^{n\times n}，Q{Q}^H=I,\Rightarrow ||QA|{|}_2=||AQ|{|}_2=||A|{|}_2$

pf:$||Q^H|{|}_2=||Q|{|}_2=1$

$||QA|{|}_2\le \Vert |Q|{|}_2||A|{|}_2=||A||,||A|{|}_2=||Q^{H}QA|{|}_2\le ||Q^H|{|}_2||QA|{|}_2=||QA|{|}_2$QED

误差分析

条件数$Cond(A)=||A||||A^{-1}||$

谱条件数$K(A)=Cond(A{)}_2=||A|{|}_2||A^{-1}|{|}_2$

条件数大:坏条件的，病态方程组

摄动问题

1.左右同时摄动

thm:设$A\in R^{n\times n}$非奇异，$||\cdot ||$是任一从属矩阵范数且假设误差$\delta A\in R^{n\times n},||A^{-1}||||\delta A||<1$若X是$Ax=b$的解，$x+\delta x$为$(A+A\delta )y=b+b\delta$的解，则(X(A)=Cond(A))
$$\frac{||\delta x||}{||x||}\le \frac{X(A)}{1-X(A)\frac{||\delta A||}{||A||}}(\frac{||\delta b||}{||b||}+\frac{||\delta A||}{||A||})$$
pf:$||A^{-1}||||\delta A||<1,A+\delta A=A(I+A^{-1}\delta A)$可逆

​ 且$(\delta A+A{)}^{-1}=(I+A^{-1}\delta A{)}^{-1}{A}^{-1}$

​ $(A+\delta A)(x+\delta x)=b+b\delta$

​ $\delta x$

$ =(A+\delta A)^{-1}(b+\delta b)-x$

$=(A+A\delta {)}^{-1}[b+\delta b-(A+\delta A)x]$

$=(I+A^{-1}\delta A{)}^{-1}{A}^{-1}(\delta b-\delta Ax)$

取范数可得上述不等式

remark:条件数对应的摄动问题实质上表现的是矩阵在各维度的拉伸比例，因为比例不均匀而带来误差。

2.舍入误差

避免大数吞小数一类的误差

def:数值稳定

列主元Gauss消去法不是数值稳定的

thm:用Gauss列主元消去法解n阶线性方程得到的计算解x是系数矩阵$A$有摄动$\delta A$的摄动方程组$(A+\delta A)x=b$的准确解，这里摄动矩阵(或称误差矩阵$\delta A$)满足估计式:
$$||\delta A|{|}_{\infty }\le 1.01(n^3+3n^2)\rho ||A|{|}_{\infty }eps$$
$\rho =\frac{ma{x}_{i,j,k}|a_{ij}^{(k)}|}{||A|{|}_{\infty }},eps$按进制估计舍入误差

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-KIYdoZ8s-1615550215966)(C:\Users\10539\Desktop\课程\数值分析\2.2矩阵和向量范数.jpg)]

Chap 3:解线性方程的迭代法

迭代法的基本理论

迭代法基本概念

解方程组$Ax=b$

$x_0,...,x_k,...$称为迭代序列

一阶线性定常迭代法:$x_k=G{x}_{k-1}+g$:迭代公式； $G:$迭代矩阵

迭代法收敛&发散

$e^{(k)}=x^{\ast }-x_k$为第k步的误差向量，若迭代法收敛，则当k充分大时可以把$x_m$当作解$x^{\ast }$的 近似向量，即方程的近似解。

$x^{\ast }$是方程$x=Gx+g$的解。若$x=Gx+g$与$Ax=b$同解则称迭代法与方程组是完全相容的。

构造迭代矩阵：$A=Q-R$,Q非奇异，称Q为分裂矩阵，从而$Ax=b\Rightarrow (Q-R)x=b$

​ $\Rightarrow x=Q^{-1}Rx+Q^{-1}b\Rightarrow G=1-Q^{-1}A;b=Q^{-1}b$

迭代法收敛条件

thm:上述迭代收敛的充要条件:$li{m}_{k\to \infty }{G}^k=O$

​ pf:$e^k=x^{\ast }-x_k=G{x}^{\ast }+g-(G{x}_{k-1}+g)=G(x^{\ast }-x_{k-1})=G^k(x^{\ast }-x_0)$

​ $e^k\to 0$等价于题目条件。

thm:上述迭代收敛等价于G的谱半径小于1，即$\rho (G)<1$ 等价于存在范数使得$||G||<1$

误差估计

对近似解$x^k$有误差估计式:$||x_k-x^{\ast }||\le ||G^k||||x^{\ast }-x_0||$

以及$||x_k-x^{\ast }||\le \frac{||G|{|}^k}{1-||G||}||x_1-x_0||$

pf:第一个不等式容易由$x^{\ast }-x_k=G^k(x^{\ast }-x_0)$得出

​ 对第二个不等式:$x_k-x^{\ast }={\sum }_{i=k}^{\infty }(x_i-x_{i+1})={\sum }_{i=k}^{\infty }G(x_{i-1}-x_i)={\sum }_{i=k}^{\infty }{G}^i(x_0-x_1)={\sum }_{i=1}^{\infty }{G}^i(x_0-x_1)$

从而得到不等式

收敛速度

$||x_k-x^{\ast }||\sim [\rho (B){]}^k||x_0-x^{\ast }||,\rho (B)$越小收敛速度越快。

确定迭代次数k使得$[\rho (G){]}^k\le 10^{-t}\Rightarrow k\ge \frac{tln10}{-ln\rho (G)}$

称$-ln\rho (G)$为相应的(渐进)收敛速度

Jacobi迭代法

原理

$Ax=b，A$非奇异，对角线无0元素，取对角线元素为Q,$A=D-(D-A),x=(1-D^{-1}A)x+D^{-1}b=Bx+g$

矩阵B恰好是对角线为0，每行其它元素是原来元素除以该行对角线元素的相反数

收敛条件

Jacobi迭代法收敛的充分必要条件是:$\rho (B)<1$

Jacobi迭代法收敛的充分条件是:存在范数使得$||B||<1$，由矩阵范数的等价性，此条件等价于对于A矩阵矩阵行/列，对角占优。

Python代码:

输入A,b,TOL(误差容限)，计算相应的B,g,进行迭代

这里D取A的对角值是为了方便计算取逆

import numpy as np
def Jacobi(A,b,tol):#A是n阶方阵array,b是n*1 array
    n=b.shape[0]
    x0=np.zeros((n,1))
    D=np.zeros((n,n))
    for i in range(n):
        D[i,i]=A[i,i]
    B=np.eye(n)-np.dot(np.linalg.pinv(D),A)
    g=np.dot(np.linalg.pinv(D),b)
    e=10
    num=0
    while e>tol:
        num+=1
        x=np.dot(B,x0)+g
        e=np.linalg.norm(x-x0)
        x0=x
        print("第",num,"次迭代结果为",'\n',x)
    return x    

b=np.array([[9],[7],[6]])
A=np.array([[10,-1,0],[-1,10,-2],[0,-4,10]])
Jacobi(A,b,0.001)
'''
第 1 次迭代结果为 
 [[0.9]
 [0.7]
 [0.6]]
第 2 次迭代结果为 
 [[0.97]
 [0.91]
 [0.88]]
第 3 次迭代结果为 
 [[0.991]
 [0.973]
 [0.964]]
第 4 次迭代结果为 
 [[0.9973]
 [0.9919]
 [0.9892]]
第 5 次迭代结果为 
 [[0.99919]
 [0.99757]
 [0.99676]]
第 6 次迭代结果为 
 [[0.999757]
 [0.999271]
 [0.999028]]
第 7 次迭代结果为 
 [[0.9999271]
 [0.9997813]
 [0.9997084]]
array([[0.9999271],
       [0.9997813],
       [0.9997084]])'''

Gauss-Seidel迭代法

原理

$Ax=b,A=D(I-L)-DU$即对Jacobi迭代法中的B做LU分解，目的是这样在每次迭代计算$x_i^k$时能利用$x_j^k,j<i$的数据。

从而有$D(I-L)x-DUx=b\Rightarrow x=(I-L{)}^{-1}Ux+(I-L{)}^{-1}{D}^{-1}bor{x}_k=L{x}_k+U{x}_{k-1}+D^{-1}b$

由于$L^{-1}$是对角线为0的下三角矩阵，从而$L{x}_k$在计算x的第i个分量时只需要$x_j^k,j<i$

收敛条件

thm:Gauss-Seidel迭代法收敛的充要条件是迭代矩阵的谱半径小于1，即
$$\rho ((I-L{)}^{-1}U)<1$$
thm:充分条件:$||B|{|}_{\infty }=ma{x}_{1\le i\le n}{\sum }_{j=1,j\ne i}^n|\frac{a_{ij}}{a_{ii}}|<1$

​ 且若记 μ = m a x 1 ≤ i ≤ n { ( ∑ j = i + 1 n ∣ a i j a i i ∣ ) / ( I − ∑ j = 1 i − 1 ∣ a i j a i i ∣ ) } \mu=max_{1\leq i\leq n}\lbrace(\sum_{j=i+1}^n|\frac{a_{ij}}{a_{ii}}|)/(I-\sum_{j=1}^{i-1}|\frac{a_{ij}}{a_{ii}}|) \rbrace μ=max1≤i≤n​{(∑j=i+1n​∣aii​aij​​∣)/(I−∑j=1i−1​∣aii​aij​​∣)} U的每行和除以1-L的每行和的最大值小于1

​ 则$\mu \le ||B|{|}_{\infty }<1$ (易得)

​ $||e^{(k)}|{|}_{\infty }=||x_k-x^{\ast }|{|}_{\infty }\le \frac{{\mu }^k}{1-\mu }||x_1-x_0|{|}_{\infty }$

​ 此外同样可证若$||B|{|}_1<1$,则Gauss-Seidel迭代法收敛

​ ( 这里$\mu$近似地看作$(I-L{)}^{-1}U$的无穷范数，但是对此待证)

python代码

实际计算中为了减小计算量实际上（***）通过向量中元素相乘求和得到的而不是向量的乘法得到的

import numpy as np
def G_S(A,b,tol):#A是n阶方阵array,b是n*1 array
    n=b.shape[0]
    x0=np.zeros((n,1))
    D=np.zeros((n,n))
    for i in range(n):
        D[i,i]=A[i,i]
    B=np.eye(n)-np.dot(np.linalg.pinv(D),A)
    g=np.dot(np.linalg.pinv(D),b)
    L=np.tril(B)
    U=np.triu(B)
    e=10
    num=0
    while e>tol:
        num+=1
        x=np.zeros((n,1))
        for i in range(n):
            x_i=np.dot(L[i],x)+np.dot(U[i],x0)+g[i]##     （***）
            x[i]=x_i
        e=np.linalg.norm(x-x0)
        x0=x
        print("第",num,"次迭代结果为",'\n',x)
    return x    

b=np.array([[9],[7],[6]])
A=np.array([[10,-1,0],[-1,10,-2],[0,-4,10]])
G_S(A,b,0.001)
'''
第 1 次迭代结果为 
 [[0.9  ]
 [0.79 ]
 [0.916]]
第 2 次迭代结果为 
 [[0.979  ]
 [0.9811 ]
 [0.99244]]
第 3 次迭代结果为 
 [[0.99811  ]
 [0.998299 ]
 [0.9993196]]
第 4 次迭代结果为 
 [[0.9998299 ]
 [0.99984691]
 [0.99993876]]
第 5 次迭代结果为 
 [[0.99998469]
 [0.99998622]
 [0.99999449]]
array([[0.99998469],
       [0.99998622],
       [0.99999449]])
'''

Jacobi迭代法和Gauss-Seidel迭代法的比较

（1）收敛速度比较

​ $||B|{|}_{\infty }<1$,则Jacobi迭代法和Gauss-Seidel迭代法都收敛，根据误差比较式看出G-S算法比Jacobi收敛速度更快。这一点从上面的例子也可以看出。

（2）收敛相关关系

两者收敛之间没有包含关系

（3）严格对角占优均收敛

point1:广义对角占优矩阵非奇异:反证存在x使得$Ax=0$，$a_{kk}{x}_k=-{\sum }_{j\ne k,j=1}^n{a}_{kj}{x}_j$取绝对值利用绝对值不等式得证

point2:迭代收敛则对于Jacobi算法:A对角占优则$||B|{|}_{\infty }<1$,收敛；对于Gauss-Seidel算法:证明:

​ pf:$G=(I-L{)}^{-1}U$,对于G的特征值$\lambda$有

​ $det(\lambda I-(I-L{)}^{-1}U)=0=|(I-L{)}^{-1}||\lambda I-\lambda L-U|=|\lambda (I-L)-U|$
$$\left[\begin{array}{cccc}\lambda & \frac{a_{12}}{a_{11}}& \frac{a_{12}}{a_{11}}& ...\\ \lambda \frac{a_{21}}{a_{22}}& \lambda & \frac{a_{23}}{a_{11}}& ...\\ \lambda \frac{a_{31}}{a_{11}}& \lambda \frac{a_{32}}{a_{11}}& \lambda & ...\\ ...& ...& ...& ...\end{array}\right]$$
​ 从而只需证明$\lambda <1$。现在反证:假设$\lambda \ge 1$,则由于A对角占优，则$|a_{ii}|\ge {\sum }_{j\ne i}|a_{ij}|$ ，从而 $ \lambda>\lambda(\sum_{j \neq i}\frac{|a_{ij}|}{|a_{ii}|})\geq \lambda(\sum_{j=1}^{{n-1}\frac{|a_{ij}|}{|a_{ii}|})+\sum_{j=i+1}}n\frac{|a_{ij}|}{|a_{ii}|}$

​ 即$(\lambda (I-L)-U)$对角占优，$|\lambda (I-L)-U|\ne 0$从而$\lambda <1$ QED

逐次超松弛迭代法(SOR方法)

原理

$Ax=b$,首先用Gauss-Seidel迭代法定义中间量${\tilde{x}}_i^{(k)}$:

$${\tilde{x}}_i=\frac{1}{a_{ii}}(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^k-\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k-1)}),i=1,2,...,n.$$

此时引入一个加速收敛的参数$\omega \in R$,称为松弛因子。然后我们取$x_i^{(k)}$为$x_i^{(k-1)}$和${\tilde{x}}_i^{(k)}$的加权平均:
$$x_i^{(k)}=x_i^{(k-1)}+\omega ({\tilde{x}}_i^k-x_i^{(k-1)})=(1-\omega )x_i^{(k-1)}+\omega {\tilde{x}}_i^k$$
def:$\omega =1$时为Gauss-Sediel迭代法

def:$\omega <1$逐次低松弛迭代法

def:$\omega >1$逐次高松弛迭代法

从而有$x_k=\omega (L{x}_k+U{x}_{k-1}+D^{-1}b)+(1-\omega )x_{k-1}\Rightarrow x_k=T_{\omega }{x}_{k-1}+\omega (I-\omega L{)}^{-1}{D}^{-1}b$

$T_{\omega }=(I-\omega L{)}^{-1}((1-\omega )I+\omega U)$

收敛条件

前提：$\forall i,a_{ii}\ne 0$

thm:充要条件:$\rho (T_{\omega })<1$

thm:必要条件:$0<\omega <2$

pf:$det(T_{\omega })=det((I-\omega )I+\omega U)=(1-\omega {)}^n$

​ $\rho (T_{\omega }{)}^n>(1-\omega {)}^n\to \rho (T_{\omega })>|1-\omega |$

thm:若A是对称正定矩阵，则当$0<\omega <2$时，解方程组$Ax=b$的SOR方法收敛

pf:$\lambda$是迭代矩阵$T_{\omega }=(I-\omega L{)}^{-1}((1-\omega )I+\omega U)$的任一特征值，从而有:$(I-\omega L{)}^{-1}((1-\omega )I+\omega U)x=\lambda x$

​ 从而$((1-\omega )I+\omega U)x=\lambda (I-\omega L)x$

​ 左右乘$x^{H}D$,记$x^{H}Dx=q,x^{H}DLx=\alpha +i\beta$又因为$A=D-DL-DU$又因为A对称正定，从而

​ $DU=(DL{)}^T,x^{H}DUx=x^H(DL{)}^{H}x=(x^{H}DLx{)}^H=\alpha -i\beta$

​ 从而有$x^{H}Ax=x^H(D-DL-DU)x=q-2a$

​ 从而$\lambda =\frac{1-\omega q+\omega \alpha -i\omega \beta }{q-\omega \alpha -i\omega \beta }$

在$0<\omega <2$时$|\lambda {|}^2<1$,从而$\rho (T_{\omega })<1$

remark:若A是对称正定矩阵Gauss-Sediel算法收敛。

最佳松弛因子

${\rho }_{{\omega }_0}(T_{\omega })=mi{n}_{\omega \in R}\rho (T_{\omega })$

python代码

import numpy as np
def SOR(A,b,w,tol,m):#A是n阶方阵array,b是n*1 array,w是松弛因子(0-2),tol误差容限，m是最大迭代次数,防止迭代不收敛
    n=b.shape[0]
    x0=np.zeros((n,1))
    D=np.zeros((n,n))
    for i in range(n):
        D[i,i]=A[i,i]
    B=np.eye(n)-np.dot(np.linalg.pinv(D),A)
    g=np.dot(np.linalg.pinv(D),b)
    L=np.tril(B)
    U=np.triu(B)
    e=10
    num=0
    while (e>tol)& (num<m):
        num+=1
        x=np.zeros((n,1))
        for i in range(n):
            x_i=w*(np.dot(L[i],x)+np.dot(U[i],x0)+g[i])+(1-w)*x0[i]
            x[i]=x_i
        e=np.linalg.norm(x-x0)
        x0=x
        print("第",num,"次迭代结果为",'\n',x)
    return x    

b=np.array([[9],[7],[6]])
A=np.array([[10,-1,0],[-1,10,-2],[0,-4,10]])
SOR(A,b,1.2,0.001,15)
'''
第 1 次迭代结果为 
 [[1.08    ]
 [0.9696  ]
 [1.185408]]
第 2 次迭代结果为 
 [[0.980352  ]
 [1.04822016]
 [0.98606408]]
第 3 次迭代结果为 
 [[1.00971602]
 [0.98817727]
 [0.99711227]]
第 4 次迭代结果为 
 [[0.99663807]
 [1.00126806]
 [1.00118621]]
第 5 次迭代结果为 
 [[1.00082455]
 [1.00013003]
 [0.99982517]]
第 6 次迭代结果为 
 [[0.99985069]
 [0.99991412]
 [0.99999374]]
第 7 次迭代结果为 
 [[1.00001956]
 [1.00001802]
 [1.0000099 ]]
array([[1.00001956],
       [1.00001802],
       [1.0000099 ]])
'''

共轭斜量法

原理

对实对称正定的矩阵A，解方程$Ax=b$可以转化为求二次函数$f(x)=\frac{1}{2}{x}^{T}Ax-b^{T}x$的极小值点的问题

thm:$u\in R^n$是方程组的解的 充要条件是u是二次函数$f(x)=\frac{1}{2}{x}^{T}Ax-b^{T}x$$f(x)=\frac{1}{2}{x}^{T}Ax-b^{T}x$的极小值点。

pf:KaTeX parse error: Undefined control sequence: \part at position 27: …wn f(x)=[\frac{\̲p̲a̲r̲t̲ ̲f}{\part x_1},.…

​ 而且$\forall p\in R^n,p\ne 0$有:$f(x+tp)-f(x)=tg(x{)}^{T}p+\frac{1}{2}{t}^2{p}^{T}Ap$,$g(u)=0,令$$x=u$可知u是f(x)的极小点

即用函数的最值来代替方程的解

共轭斜量法算法

找到初值点$x_0$

下一步$x=x_0+t{p}_0$

之后$x=x_k+t{p}_k$,其中$p_k$为寻查方向

记${\phi }_k(t)=f(x_k+t{p}_k)=\frac{1}{2}(x_k+t{p}_k{)}^{T}A(x_k+t{p}_k)-b^T(x_k+t{p}_k)$,要使得固定$p_k$时${\phi }_k$最小，对t求导得

remark:其中scalar、vector、matrix之间的求导有两种形式，具体见:

• 5.1Numerator-layout notation

• 5.2Denominator-layout notation

${\phi }_k^{\prime }(t)=t{p}_k^{T}A{p}_k+p_k^T(A{x}_k-b)$从而${\phi }_k^{\prime }(t)=0\Rightarrow t={\alpha }_k=-\frac{p_k^T(A{x}_K-b)}{p_k^{T}A{p}_k}$由于${\phi }_k^{\prime \prime }(t)>0$

从而得到对应的t的值。

这样得到的$x_{k+1}$就是直线$x=x_k+t{p}_k$上的极小值点。

剩余向量$r_k=A{x}_k-b$,有$p_{k}TA{p}_k$ $r_k=grad(f(x_k))$从而$t={\alpha }_k=-\frac{r_k^T{p}_k}{p_{k}TA{p}_k}$

显然$f(x_{k+1})\le f(x_k)$

为了保证等号不成立，选择寻查方向$p_0,p_1,...,p_{n-1}$为$R^n$中A共轭向量系，即:
$$p_i^{T}A{p}_j=0,i\ne j$$
即在内积$<Lx,Ly>=x^{T}Ay,A=L{L}^T$意义下正交,所以每次要做正交化:$p_0=-r_0;p_1=-r_1-{\beta }_0{r}_0,s.t.p_1^{T}A{p}_0=0\Rightarrow {\beta }_0=\frac{r_1^{T}A{p}_0}{p_0^{T}A{p}_0}$

${\beta }_k=\frac{r_{k+1}^{T}A{p}_k}{p_k^{T}A{p}_k}$

且可以证明这样得到的$p_0,...,p_k,(k\le n-1)$是共轭的。

remark:整个操作类似Schmidt正交化，但是需要注意的是，在求下一步共轭向量时只需要在对应x的梯度上减去$\beta$乘以上一个向量。

thm:解n阶方程组的共轭方向法至多进行n步便可得到方程的解(共轭斜量法实际上是解线性方程组的直接方法)

thm:若矩阵A只有m个不同特征值，则共轭斜量法至多进行m步即可得到方程组的解

共轭斜量法计算误差估计

$r=Ax-b=Ax-A{x}^{\ast }=A(x-x^{\ast }),x-x^{\ast }=A^{-1}r$

$||x-x^{\ast }||=||A^{-1}r||\le ||A^{-1}||||r||,||b||\le ||A||||x^{\ast }||$

$\frac{||x-x^{\ast }||}{||x^{\ast }||}\le ||A||||A^{-1}||\frac{||r||}{||b||}=cond(A)\frac{||r||}{||b||}$

从而坏条件方程组剩余向量的大小不能说明近似解的精确程度。

最速下降算法

相比于共轭斜量法需要保证向量的共轭来达到把正交维度的向量抵消的目的，

最速下降法不用着眼于全局而只需要在目前的点处进行梯度下降。

直接取$p_k=-r_k=b-a{x}_k$

​ $x_{k+1}=x_k+{\alpha }_k{p}_k,r_k=▽f(x_k)$

​ ${\alpha }_k=-\frac{r_k^T{p}_K}{p_k^{T}A{p}_k}=\frac{r_k^T{r}_k}{r_k^{T}A{r}_k}$

thm:设$A\in R^{n\times n}$对称正定且其特征值为${\lambda }_1\ge {\lambda }_2\ge ...\ge {\lambda }_n>0$则对于任意初始向量$x_0$最速下降法都收敛，且$A{x}^{\ast }=b$,$||x_k-x^{\ast }|{|}_A\le |\frac{{\lambda }_1-{\lambda }_n}{{\lambda }_1+{\lambda }_n}{|}^k||x_0-x^{\ast }|{|}_A$

其中$||x|{|}_A=(Ax,x{)}^{\frac{1}{2}}$

最速下降算法python实现

import numpy as np
def RD(A,b,tol,m,x0):#A是n阶方阵array,b是n*1 array,w是松弛因子(0-2),tol误差容限，m是最大迭代次数,防止迭代不收敛
    e=10
    num=0
    x=x0
    while (e>tol)& (num<m):
        num+=1
        r=np.dot(A,x)-b
        a=(np.dot(r.T,r))/np.dot(np.dot(r.T,A),r)
        x=x0-a*r
        e=np.linalg.norm(x-x0)
        x0=x
        print("第",num,"次迭代结果为",'\n',x)
    return x 

b=np.array([[9],[7],[6]])
A=np.array([[10,-1,0],[-1,10,-2],[0,-4,10]])
x0=np.array([[0],[0],[0]])
RD(A,b,0.001,15,x0)
'''
第 1 次迭代结果为 
 [[1.16536661]
 [0.90639626]
 [0.77691108]]
第 2 次迭代结果为 
 [[0.97702454]
 [0.97702454]
 [0.97702454]]
第 3 次迭代结果为 
 [[1.00379937]
 [0.99784941]
 [0.99487443]]
第 4 次迭代结果为 
 [[0.99947213]
 [0.99947213]
 [0.99947213]]
第 5 次迭代结果为 
 [[1.00008729]
 [0.99995059]
 [0.99988224]]
array([[1.00008729],
       [0.99995059],
       [0.99988224]])
'''

共轭斜量法python实现

import numpy as np
def HG(A,b,tol,m,x0):#A是n阶方阵array,b是n*1 array,w是松弛因子(0-2),tol误差容限，m是最大迭代次数,防止迭代不收敛
    n=b.shape[0]
    e=10
    num=0
    x=x0
    r=np.dot(A,x)-b;p=-r
    while (e>tol)& (num<m):
        q=np.dot(A,p)#预先计算
        pAp=np.dot(p.T,q)#减少反复运算计算量
        a=-(np.dot(r.T,p))/pAp
        x=x0+a*p
        r=r+a*q#r改为迭代形式
        beta=np.dot(r.T,q)/pAp
        p=-r+beta*p
        
        e=np.linalg.norm(x-x0)
        x0=x 
        num+=1
        print("第",num,"次迭代结果为",'\n',x)
    return x 

b=np.array([[9],[7],[6]])
A=np.array([[10,-1,0],[-1,10,-2],[0,-4,10]])
x0=np.array([[0],[0],[0]])
HG(A,b,0.001,15,x0)
'''
第 1 次迭代结果为 
 [[1.16536661]
 [0.90639626]
 [0.77691108]]
第 2 次迭代结果为 
 [[1.01404056]
 [1.01092044]
 [1.00936037]]
第 3 次迭代结果为 
 [[0.99938234]
 [1.00042663]
 [1.00094877]]
第 4 次迭代结果为 
 [[0.99994017]
 [0.99995717]
 [0.99996567]]
第 5 次迭代结果为 
 [[1.00000228]
 [0.99999808]
 [0.99999599]]
array([[1.00000228],
       [0.99999808],
       [0.99999599]])
'''

Chap 4:插值法

overview

$f(x)={\sum }_{k=0}^n\frac{f^{(k)}(x_0)}{k!}(x-x_0{)}^k+R_n(x)=P_n(x)+R_n(x),R_n(x)=\frac{f^{(n+1)(\xi )}}{(n+1)!}(x-x_0{)}^{(n+1)}$

如果$P_n(x_0)=f(x_0),f_n^{\prime }(x_0)=f^{\prime }(x_0),...$则是Hermite插值

如果$x_0,...,x_n$为n+1个互异点(插值基点/节点)找函数Z$p_n(x)\in P_n(x)=R_{n+1}[x]$

其中 P n ( x ) = { ∑ k = 0 n a k x k ∣ a k ∈ R } P_n(x)=\lbrace \sum_{k=0}^n a_kx^k|a_k\in R \rbrace Pn​(x)={∑k=0n​ak​xk∣ak​∈R},s.t. $p_n(x_i)=f(x_i)$则是Lagrange插值

[a,b]为插值区间$a=x_0\le ...\le x_n=b,x\in [a,b]$:内插，$x\notin [a,b]$:外插

选择多项式为插值函数:$g(x)$为插值多项式，代数插值 g ( x ) ∈ s p a n { x 0 , x 1 , . . . , x n } g(x)\in span\lbrace x^0,x^1,...,x^n\rbrace g(x)∈span{x0,x1,...,xn}

${\gamma }_n(x)=f(x)-p_n(x)$:插值余项

代数插值的存在唯一性:

$p_n(x_i)=f(x_i),i=0,1,2,...,n$

存在唯一性由范德蒙行列式的存在唯一性保证

Method1:Lagrange插值公式

构造Lagrange多项式

形式:

recite

$p_n=y_0{l}_0+y_1{l}_1+y_2{l}_2+...+y_n{l}_n{y}_i=f(x_i),l_i={\Pi }_{j=0,j\ne i}^n\frac{(x-x_j)}{(x_i-x_j)}$

或者:$w_{n+1}(x)=(x-x_0)(x-x_1)...(x-x_n),$

​ $w_{n+1}^{\prime }(x_i)$

​ $=li{m}_{x\to x_i}\frac{w_{n+1}(x)-w_{n+1}(x_i)}{x-x_i}$

​ $=(x_i-x_0)(x_i-x_1)...(x_i-x_{i-1})(x_i-x_{i+1})...(x_i-x_n)$

​ $l_i=\frac{w_{n+1}(x)}{(x-x_i)w_{n+1}^{\prime }(x_i)}$

插值公式的余项

thm:f(x)包含n+1个互异基点$x_0,x_1,...,x_n$在区间[a,b]之间且有n阶连续导数，在(a,b)有n+1阶有界导数，则$\forall x\in [a,b]$存在$\xi \in (a,b),s.t.$
$$r_n(x)=\frac{f^{n+1}(\xi )}{(n+1)!}{w}_{n+1}(x)$$
pf:引入辅助函数$g(t)=f(t)-f_n(t)-K(x)w_{n+1}(t),r_n(x)=K(x)w_{n+1}(x)$

​ 从而$t=x,x_0,x_1,...,x_n$时$g(t)=0$反复使用Rolle定理得证。

remark:其中$w_{n+1}(x)$还可以通过不等式进一步放缩。例如 r 1 ( x ) = 1 2 f ′ ′ ( ξ ) ( x − a ) ( x − b ) ≤ ( a − b ) 2 8 m a x { ∣ f ′ ′ ( x ) ∣ } r_1(x)=\frac{1}{2}f''(\xi)(x-a)(x-b)\leq \frac{(a-b)^2}{8}max\lbrace |f''(x)|\rbrace r1​(x)=21​f′′(ξ)(x−a)(x−b)≤8(a−b)2​max{∣f′′(x)∣}

n取的大不一定能保证插值公式的逼近效果好。

python实现插值求x处函数值

import numpy as np
def Lagrange(X,Y,x):
    l=[]
    k=len(X)
    for i in range(k):
        temp=1
        for j in range(k):
            if j!=i:
                temp=temp*(x-X[j])/(X[i]-X[j])
        l.append(temp)
    y=sum(list(map(lambda x,y:x*y,l,Y)) )
    return y

a=lambda x:x**3+3*x**2+38
X=[3,4,-2,1]
Y=list(map(a,X))
a(0)# =38

Lagrange(X,Y,0)
'''37.99999999999999'''

Method2:Newton插值公式

由于Lagrange插值公式需要一次性知道n+1个点的信息，这样在逐步得到点的信息并输入时会造成计算的复杂度的增加，从而采用Newton插值公式，这样增加一个点，插值多项式的次数+1

$N_n(x)=a_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+...+a_n(x-x_0)(x-x_1)...(x-x_{n-1})$

均差

$a_0=f(x_0),a_1=\frac{f(x_1)-f(x_0)}{x_1-x_0},a_2=[\frac{f(x_2)-f(x_1)}{x_2-x_1}-\frac{f(x_1)-f(x_0)}{x_1-x_0}]\frac{1}{x_2-x_0}...$

差商、一阶均差:$f[x_i,x_j]=\frac{f(x_j)-f(x_i)}{x_j-x_i}{a}_1=f[x_0,x_1]$

二阶均差:$f[x_i,x_j,x_k]=\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k}{a}_2=f[x_0,x_1,x_2]$

…

零阶均差:$f[x_i]=f(x_i)$

可以表示为$f[x_0,...,x_n]={\sum }_{i=0}^n\frac{f(x_i)}{{\Pi }_{j=0,i\ne j}^n(x_i-x_j)}$

info:1.均差有对称性

​ 2.$x_1,...x_n$是n次多项式的n个互异实根，从而$f[x_1,...,x_k]=0,(1\le k\le n)$

​ 3.$f[x,x_0,...,x_k]$是x的m次多项式，则 $f[x,x_0,...,x_k,x_{k+1}]$是x的m-1次多项式

Newton均差插值多项式

(要求记忆)

$f(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)$

​ $+...+f[x_0,...,x_n](x-x_0)...(x-x_{n-1})+f[x,x_0,...,x_n](x-x_0)...(x-x_n)$

​ $=f[x_0]+f[x_0,x_1]w_1(x)+...+f[x_0,...,x_n]w_n(x)+f[x,x_0,...,x_n]w_{n+1}(x)$

​ $=N_n(x)+f[x,x_0,...,x_n]w_{n+1}(x)$

由代数插值的唯一性$f[x,x_0,...,x_n]=\frac{f^{n+1}(\xi )}{(n+1)!},\xi \in (a,b)$

从而有$f[x,x,x]=\frac{f^{(2)}(x)}{2!}$

import numpy as np
def Newton(X,Y,x):
    l=[]
    k=len(X)
    Q=np.zeros((k,k))
    Q[:,0]=Y
    for i in range(k-1):
        for j in range(k-i-1):
            Q[j+i+1,i+1]=(Q[j+i,i]-Q[j+i+1,i])/(X[j]-X[j+i+1])
    y=0
    for i in range(k):
        y=Q[k-i-1,k-i-1]+y*(x-X[k-1-i])
    return y

Newton(X,Y,0)
'''38.0'''

Method2*:有限差

在Newton插值公式中因为$x_0,x_1,...,x_n$之间的距离不同而需要多次对距离差做商，可以运用有限差使得样本x等距，这样分子上$w=(x-x_0)(x-x_1),...$能够与分母$h^k$消元，用以减小计算量。

一阶向前差分$△f(x)=f(x+h)-f(x)$

一阶向后差分$▽f(x)=f(x)-f(x-h)$

一阶中心差分$\delta f(x)=f(x+\frac{h}{2})-f(x-\frac{h}{2})$

高阶情况:${△}^{n}f(x)=△({△}^{n-1}f(x))$

prop:

(1)有限差算子为线性算子

(2)用函数值表示高阶有限差 P134

(3)用有限差表示函数值 P134

均差与有限差的关系

$x_i-x_{i-1}=h$

$f[x_0,x_1,...,x_n]=\frac{{△}^{n}f(x_0)}{n!h^n}=\frac{{▽}^{n}f(x_n)}{n!h^n}$

用归纳法证明

Newton前差和后差插值公式

计算$x=x_0+sh$处的函数值(s不一定为整数)

利用Newton插值公式和均差与有限差之间的关系有:

（1）带余项的Newton前差插值公式

$f(x)=f(x_0+sh)=f(x_0)+s△f(x_0)+\frac{s(s-1)}{2!}{△}^{2}f(x_0)+...+\frac{s(s-1)...(s-n+1)}{n!}{△}^{n}f(x_0)+\frac{s(s-1)...(s-n)}{(n+1)!}{h}^{n+1}{f}^{(n+1)}(\xi )$

$\xi \in (x_0，x_n)$

计算$x=x_n+th,t\le 0$处的函数值(s不一定为整数)

利用Newton插值公式和均差与有限差之间的关系有:

（2）带余项的Newton后差插值公式

$f(x)=f(x_n+th)=f(x_n)+t▽f(x_n)+\frac{t(t+1)}{2!}{▽}^{2}f(x_n)+...+\frac{t(t+1)...(t+n+1)}{n!}{▽}^{n}f(x_n)+\frac{t(t+1)...(t+n)}{(n+1)!}{h}^{n+1}{f}^{(n+1)}(\xi )$

$\xi \in (x_0，x_n)$

Python计算前差插值公式求解x处近似值

import numpy as np
def Newton_forward(h,X,Y,x):
    s=(x-X[0])/h
    k=len(X)#k=n+1
    Q=np.zeros((k,k))
    Q[:,0]=Y
    for i in range(k-1):
        for j in range(k-i-1):
            Q[j+i+1,i+1]=(-Q[j+i,i]+Q[j+i+1,i])
    print(Q)
    y=0
    for i in range(k):
        y=Q[k-i-1,k-i-1]+y*(s-k+i+1)/(k-i)
    return y

a=lambda x:x**3+3*x**2+38
X=[-0.3,0.2,0.7,1.2]
Y=list(map(a,X))
a(0.5) #a(0.5)=38.875

Newton_forward(0.5,X,Y,0.5)
'''
[[38.243  0.     0.     0.   ]
 [38.128 -0.115  0.     0.   ]
 [39.813  1.685  1.8    0.   ]
 [44.048  4.235  2.55   0.75 ]]
38.875
'''

Python计算后差插值公式求解x处近似值

import numpy as np
def Newton_backward(h,X,Y,x):
    t=(x-X[-1])/h
    k=len(X)#k=n+1
    Q=np.zeros((k,k))
    Q[:,0]=Y
    for i in range(k-1):
        for j in range(k-i-1):
            Q[j+i+1,i+1]=(-Q[j+i,i]+Q[j+i+1,i])
    y=0
    for i in range(k):
        y=Q[-1,k-i-1]+y*(t-(-k+i+1))/(k-i)
    return y
   
Newton_backward(0.5,X,Y,0.5)
'''
38.875
'''

Method3:Hermite插值公式

原理

已知某点函数值和若干阶导数情况下的插值

下面只考虑:
已知$x_0,...,x_n$处的函数值和一阶导数值时的插值，此时用2n+1阶多项式$H_{2n+1}$进行逼近

使得$H_{2n+1}$满足:
$$\{\begin{array}{c}{H}_{2n+1}(x_j)=f(x_j)\\ H_{2n+1}^{\prime }=f^{\prime }(x_j)\end{array}$$
构造:$H_{2n+1}(x)={\sum }_{i=0}^{n}f(x_i)A_i(x)+{\sum }_{i=0}^n{f}^{\prime }(x_i)B_i(x)$

使得$A_i(x_j)={\delta }_{ij},A^{\prime }(x_j)=0$

​ $B_i(x_j)=0,B_i^{\prime }(x_j)={\delta }_{ij}$

从而$B_i(x)={\beta }_i(x-x_i)(x-x_0{)}^2...(x-x_n{)}^2\to B_i(x)=(x-x_i)l_i^2(x)$

从而$A_i(x)={\alpha }_i(ax+b)(x-x_0{)}^2...(x-x_{i-1}{)}^2(x-x_{i+1}^2)...(x-x_n{)}^2$

​ $\to A_i(x)=[1-2(x-x_i){\sum }_{j=0,j\ne i}^n\frac{1}{x_i-x_j}]l_i^2(x)$

从而得到Hermite插值公式(不需要记)

余项

假设$f$在$[a,b]$上有2n+1阶连续导数，在$(a,b)$上存在$2n+1$阶导数，则$forx\in [a,b],\exists \xi \in (a,b),s.t.$
$$f(x)-H_{2n+1}(x)=\frac{w_{2n+1}^2(x)}{(2n+2)!}{f}^{(2n+2)}(\xi )$$
有n+1个二重零点，引入辅助函数$F(t)=f(t)-H_{2n+1}(t)-K(x)w_{2n+1}^2(t)$

已知一些点的值&一些点的导数，求多项式逼近

n个点的值知道n-1次多项式，加一个值的导数+一次

一般先以$Lagrange$作为基础，用待定系数法求。

Method4:样条插值

了解概念、构造方法

分段插值

一次:$g_{1,i}(x)=f(x_i)+f[x_i,x_{i+1}](x-x_i)$

分段抛物线插值:$g_{2,i}(x)=f(x_i)+f[x_i,x_{i+1}](x-x_i)+f[x_i,x_{i+1},x_{i+2}](x-x_i)(x-x_{i+1})$

求对应的误差限（利用Taylor展开，详见P143）

三次样条插值

P142 分段插值，且每段是不超过三次的多项式。

满足$(1)S(x_j)=f(x_j),(2)$ 连接点处光滑$S_{i-1}^{\prime }(x_i-0)=S_i^{\prime }(x_i+0)$

​ $(3)$边界条件/端点条件

先由$S^{\prime \prime }(x)=m_i\frac{x_{i+1}-x}{h_i}+m_{i+1}\frac{x-x_i}{h_i},x\in [x_i,x_{i+1}]$

积分由$(1)$得到两次积分常数C的值，由$(2)$在n-1个连接处得到$m_1,...,m_{n+1}$的n-1个方程

加上$(3)$求解三对角矩阵

(1)已知边界点的二阶导数

$f^{\prime \prime }(a)=m_1,f^{\prime \prime }(b)=m_{n+1}$

如果$S^{\prime \prime }(x_1)=S^{\prime \prime }(x_{n+1})=0$则为自然三次样条插值

(2)已知边界点的导数

$f^{\prime }(a)=S^{\prime }(a),f^{\prime }(b)=S^{\prime }(b)$

这样解$m_1,...,m_{n+1}$的方程是对角占优的，从而有唯一解，这样的插值函数称为完备三次样条插值

(3)对Newton插值多项式

由$x_1,...,x_4$做牛顿插值$N_a(x)$

由$x_{n+1},x_n,...,x_{n-2}$做牛顿插值$N_b(x)$

取$S^{\prime }(a)=N^{\prime }(a),S^{\prime }(b)=N^{\prime }(b)$

这样的插值叫做Lagrange三次样条插值

Chap 5:函数逼近

5.1 基本概念

def:广义多项式，最小多项式

def:$C[a,b]$上的函数范数

(1)非负性 (2)齐次性 (3)三角不等式

def:距离$||f(x)-g(x)||$

权函数

$W(x)\in C[a,b]$:

(1)$W(x)\ge 0,\forall x\in [a,b]$且最多有有限个零点

(2)${\int }_a^b{x}^{n}W(x)dx,n=0,1,...$存在

则称$W(x)$是关于区间$[a,b]$的权函数

remark:对应Euclid范数:$||f(x)|{|}_2=({\int }_a^b{f}^2(x)W(x)dx{)}^{\frac{1}{2}}$

逼近

$\phi \in span({\phi }_0,...,{\phi }_n)\subset C[a,b]$

(1)一致逼近:$li{m}_{n\to \infty }||f(x)-{\phi }_{(}x)|{|}_{\infty }=0$

(2)平方逼近:$li{m}_{n\to \infty }||f(x)-{\phi }_{(}x)|{|}_2=0$

​ 即$li{m}_{n\to \infty }{\int }_a^b[f(x)-g(x){]}^{2}W(x)dx=0$

​ 带权$W(x)$的平方逼近

5.2最佳一致逼近

Bernstein 多项式,但是收敛慢。

Chebyshev:在不超过n的多项式中找最佳逼近

remark: H n = s p a n { 1 , x , x 2 , . . . , x n } H_n=span\lbrace 1,x,x^2,...,x^n\rbrace Hn​=span{1,x,x2,...,xn}

def:

bias:$||f(x)-p_n(x)|{|}_{\infty }=ma{x}_{a\le x\le b}|f(x)-p_n(x)|$

n次最佳逼近/最小偏差:$E_n=mi{n}_{p_n\in H_n}ma{x}_{a\le x\le b}|f(x)-p_n(x)|$

​ 其中使得bias取得$E_n$的多项式$p_n(x)$称为$[a,b]$上的n次最佳一致逼近多项式(简称最佳逼近多项式)

​ remark:由Bernstein定理知$n\to \infty ,E_n\to 0$

交错点组:$g(x)\in C[a,b],a\le x_1<x_2<...<x_n\le b,|g(x_i)|=ma{x}_{a\le x_i\le b}|g(x)|$称$x_1,...,x_n$为[a,b]上的交错点组。

thm:

(1)Chebyshev定理:设$f\in C[a,b],p(x)\in H_n$,则$p(x)$是$f(x)$的最佳一致逼近多项式的充要条件是$f(x)-p(x)$在$[s,b]$上至少有n+2个点组成的交错点组

(2)$f(x)\in C[a,b]$,则在$H_n$中，$f(x)$有唯一的一个最佳一致逼近多项式$p(x)$

(3)$f\in C[a,b]$有n+1阶导数，且$f^{(n+1)}(x)$在$(a,b)$内保持定号，$p(x)\in H_n$是$f(x)$的最佳一致逼近多项式，则区间$[a,b]$的端点属于$f(x)-p(x)$的交错点组 （这一点可以由Rolle定理证明）

apply:

最佳一致逼近的计算:利用上面定理。方程依据:

(1)交错点组上的点是区间内的极值点。

​ 从而如果高一阶导数定号，则已确定两点，其它内点为极值点$g^{\prime }(x_i)=0$

​ 如果高一阶导数不定号再做讨论

(2)n+2个点的函数值相同，从而可以列出n-2个等式

5.3 最佳平方逼近

原理

关于权函数$W(x)$的最佳平方逼近/最小二乘逼近$\phi (x):ar{g}_{\phi (x)}min{\int }_a^b[f(x)-\phi (x){]}^{2}W(x)dx$

其中$\phi (x)={\sum }_{j=0}^n{a}_j{\phi }_j(x)$

$T={\int }_a^b[f(x)-\phi (x){]}^{2}W(x)dx$对$a_i$求偏导等于0，从而$f(x)-\phi (x)$与${\phi }_k(x)$在$W(x)$二范数意义下正交。

即$f(x)-\phi (x)$与其它所有 s p a n { φ 1 , . . . φ n } span\lbrace \varphi_1,...\varphi_n\rbrace span{φ1​,...φn​}上的函数逼近正交

从而$({\phi }_i,{\phi }_0)a_0+({\phi }_i,{\phi }_1)a_1+...+({\phi }_i,{\phi }_n)a_n=({\phi }_i,f)$

得到Gram矩阵

矩阵的解就是最佳平方逼近(相关证明可由函数内积给出)，Gram矩阵形成的 方程称为法方程。P162

具体计算细节

为了方便解法方程，取正交函数系，从而$a_j=\frac{({\phi }_j,f)}{||{\phi }_j||}$

称$a_j$为$f(x)$关于正交函数系的Fourier系数，级数${\sum }_{j=0}^{\infty }{a}_j{\phi }_j(x)$称为f(x)的广义Fourier级数

5.4 直交多项式

def：直交多项式系$p_1,...,p_n$,在范数$<p_1,p_2>=\int p_1{p}_{2}W(x)dx$下直交

首一直交多项式系算法

$q_0=1,q_{k+1}=(x-{\alpha }_k)q_k-{\beta }_k{q}_{k-1}$

${\alpha }_k=\frac{\int x[q_k(x){]}^{2}W(x)dx}{\int [q_k(x){]}^{2}W(x)dx}$

$\beta =\frac{\int [q_k(x){]}^{2}W(x)dx}{\int [q_{k-1}(x){]}^{2}W(x)dx},k\ge 1$

证明:$q_{k+1}=(x-c_k)q_k+c_{k-1}{q}_{k-1}+...+c_0{q}_0$

PS:其中$q_k$与任何次数低于k的多项式都直交

直交多项式$P_n$在区间有n个互异的根

因为和常数$p_0$积分为0,如果r<n,从而在$(a,b)$有奇数重实根$x_1,...,x_r,g=(x-x_1)...(x-x_r)$从而$p_n(x)g_r(x)W(x)\ge 0,\int p_n(x)g_r(x)W(x)dx>0$，从而矛盾，所以$r=n$4

5.4* Chebyshev多项式

$==T_n(x)=\frac{1}{2}[(x+\sqrt{x^2-1}{)}^n+(x-\sqrt{x^2-1}{)}^n]=cos(narccosx)(x=cos\theta 可证)$==

递推关系

$T_{n+1}=2x{T}_n-T_{n-1},T_0=1,T_1=x$

prop:

• 首项系数$2^{n-1}$
• 奇偶性:n为奇数是奇函数，n为偶数是偶函数
• $|x|\le 1,|T_n(x)|\le 1;|x|>1,|T_n(x)>1|;x\in [-1,1],max|T|=1$
• 存在n+1个点$x_k=cos(\frac{k\pi }{n}),T_n(x_n)=(-1{)}^k$交错点组
• 存在n个根$x_k^n=cos\frac{(2k-1)\pi }{2n},k=1,...,n$
• { T k ( x ) } \lbrace T_k(x)\rbrace {Tk​(x)}是$W(x)=(1-x^2{)}^{-\frac{1}{2}}$的直交多项式

对应最佳一致逼近

对所有次数不超过n$q_n(z)=1，ma{x}_{x\in [-1,1]}|q_n(z)|$最小时${\bar{T}}_n(x)=\frac{T_n(x)}{T_n(z)}$

从而对于首一多项式$ma{x}_{x\in [-1,1]}|p_n(x)|$最小时对应多项式为$\frac{T_n(x)}{2^{n-1}}$偏差为$\frac{1}{2^{n-1}}$

求1/2/3次最佳一致逼近多项式

5.4**Legendre多项式

$P_n(x)=\frac{d^n(x^2-1{)}^n}{2^{n}n!d{x}^n}$

prop:

• 首项系数$\frac{(2n)!}{2^n(n!{)}^2}$
• ${\int }_{-1}^1{P}_m(x)P_n(x)dx=0(m\ne n),\frac{2}{2m+1}(m=n)$积分过程见P175 ps:${\int }_0^{\frac{\pi }{2}}co{s}^{2m+1}\theta d\theta =\frac{(2m)!!}{(2m+1)!!}$
• $P_n(-x)=(-1{)}^n{P}_n(x)$

递推关系

$P_{n+1}=\frac{2n+1}{n+1}x{P}_n(x)-\frac{n}{n+1}{P}_{n-1}(x),P_0=1,P_1=x$

最佳平方逼近

在首一的n次实系数多项式中首一Legendre多项式${\tilde{P}}_n(x)=\frac{2^n(n!{)}^2}{(2n)!}{P}_n(x)$在[-1,1]上和0的平方误差最小，即内积最小

5.4***其它

对权函数$W(x)=e^{-x}$对应直交多项式:Laguerre多项式

对权函数$W(x)=e^{-x^2}$对应直交多项式:Hermite多项式

5.5近似最佳一致逼近

Lagrange插值多项式

$r_n(x)=f(x)-p_n(x)=\frac{f^{(n+1)}(\xi )}{(n+1)!}{w}_{n+1}(x)$

$w_{n+1}(x)=(x-x_0)(x-x_1)...(x-x_n)$

选取Cheyshev多项式零点作为插值基点减小$w_{n+1}(x)$

插值区间的转换&余项的求法P180

幂级数项数缩短

$S_n=b_0{T}_0+b_1{T}_1+...+b_n{T}_n$

5.6 函数按照直交多项式展开

chebyshev级数$\frac{1}{2}{a}_0+{\sum }_{j=1}^{\infty }{a}_j{T}_j(x),a_j=\frac{(T_j,f)}{(T_j,T_j)}=\frac{2}{\pi }{\int }_0^{\pi }f(cos\theta )cosj\theta d\theta$

类似在权函数为1时用Legendre做展开但是$P(x)=a_0{P}_0(x)+a_1{P}_1(x)+...$

Chap 6:最小二乘拟合

离散的最佳二次逼近

Chebyshev多项式作为拟合基函数

Chap 7 数值积分

利用$x_1,x_2,...,x_{n+1}$作为求积基点(节点)，$A_1,...,A_{n+1}$作为求积系数，$I(f)\simeq I_n(f)={\sum }_i{A}_{i}f(x_i)$

余项/离散误差$E_n(f)=I(f)-I_n(f)$

7.1 Newton-Cotes

插值求积公式

用Lagrange插值多项式$f_n(x)$逼近$f(x)$,对$f_n(x)$积分代替$\int f(x)dx$

求法&误差

Newton-Cotes求积公式

等步长求积

求法+误差

梯形公式和Simpson公式

梯形公式

一次函数拟合

Simpson公式

二次函数拟合

误差分析

• 离散误差
• 计算误差($f(x_i)$理论值与实际值的误差)

7.2 复合求积公式

把区间分段，每段分别求积分

复合梯形公式

误差

变步长梯形公式

误差，公式，精度

复合Simpson公式

7.5 Gauss 型数值求积公式

代数精确度

$\forall j:I_n(x_j^k)=I(x_j^k),k\le n$求满足一定代数精确度的积分的系数与基点

最小是n，最大是2n+1（证明上界:取$f(x)=w_{n+1}^2(x)$）

Gauss型求积公式

先做插值展开，选择基点为Legendre多项式零点

误差P239

Gauss-Legendre求积公式

W=1 计算积分

Gauss-Chebyshev求积公式

代数精确度2n-1

n[-1,1]}|q_n(z)|$最小时$\bar{T}_n(x)=\frac{T_n(x)}{T_n(z)}$

从而对于首一多项式$ma{x}_{x\in [-1,1]}|p_n(x)|$最小时对应多项式为$\frac{T_n(x)}{2^{n-1}}$偏差为$\frac{1}{2^{n-1}}$

求1/2/3次最佳一致逼近多项式

5.4**Legendre多项式

$P_n(x)=\frac{d^n(x^2-1{)}^n}{2^{n}n!d{x}^n}$

prop:

• 首项系数$\frac{(2n)!}{2^n(n!{)}^2}$
• ${\int }_{-1}^1{P}_m(x)P_n(x)dx=0(m\ne n),\frac{2}{2m+1}(m=n)$积分过程见P175 ps:${\int }_0^{\frac{\pi }{2}}co{s}^{2m+1}\theta d\theta =\frac{(2m)!!}{(2m+1)!!}$
• $P_n(-x)=(-1{)}^n{P}_n(x)$

递推关系

$P_{n+1}=\frac{2n+1}{n+1}x{P}_n(x)-\frac{n}{n+1}{P}_{n-1}(x),P_0=1,P_1=x$

最佳平方逼近

在首一的n次实系数多项式中首一Legendre多项式${\tilde{P}}_n(x)=\frac{2^n(n!{)}^2}{(2n)!}{P}_n(x)$在[-1,1]上和0的平方误差最小，即内积最小

5.4***其它

对权函数$W(x)=e^{-x}$对应直交多项式:Laguerre多项式

对权函数$W(x)=e^{-x^2}$对应直交多项式:Hermite多项式

5.5近似最佳一致逼近

Lagrange插值多项式

$r_n(x)=f(x)-p_n(x)=\frac{f^{(n+1)}(\xi )}{(n+1)!}{w}_{n+1}(x)$

$w_{n+1}(x)=(x-x_0)(x-x_1)...(x-x_n)$

选取Cheyshev多项式零点作为插值基点减小$w_{n+1}(x)$

插值区间的转换&余项的求法P180

幂级数项数缩短

$S_n=b_0{T}_0+b_1{T}_1+...+b_n{T}_n$

5.6 函数按照直交多项式展开

chebyshev级数$\frac{1}{2}{a}_0+{\sum }_{j=1}^{\infty }{a}_j{T}_j(x),a_j=\frac{(T_j,f)}{(T_j,T_j)}=\frac{2}{\pi }{\int }_0^{\pi }f(cos\theta )cosj\theta d\theta$

类似在权函数为1时用Legendre做展开但是$P(x)=a_0{P}_0(x)+a_1{P}_1(x)+...$

Chap 6:最小二乘拟合

离散的最佳二次逼近

Chebyshev多项式作为拟合基函数

Chap 7 数值积分

利用$x_1,x_2,...,x_{n+1}$作为求积基点(节点)，$A_1,...,A_{n+1}$作为求积系数，$I(f)\simeq I_n(f)={\sum }_i{A}_{i}f(x_i)$

余项/离散误差$E_n(f)=I(f)-I_n(f)$

7.1 Newton-Cotes

插值求积公式

用Lagrange插值多项式$f_n(x)$逼近$f(x)$,对$f_n(x)$积分代替$\int f(x)dx$

求法&误差

Newton-Cotes求积公式

等步长求积

求法+误差

梯形公式和Simpson公式

梯形公式

一次函数拟合

Simpson公式

二次函数拟合

误差分析

• 离散误差
• 计算误差($f(x_i)$理论值与实际值的误差)

7.2 复合求积公式

把区间分段，每段分别求积分

复合梯形公式

误差

变步长梯形公式

误差，公式，精度

复合Simpson公式

7.5 Gauss 型数值求积公式

代数精确度

$\forall j:I_n(x_j^k)=I(x_j^k),k\le n$求满足一定代数精确度的积分的系数与基点

最小是n，最大是2n+1（证明上界:取$f(x)=w_{n+1}^2(x)$）

Gauss型求积公式

先做插值展开，选择基点为Legendre多项式零点

误差P239

Gauss-Legendre求积公式

W=1 计算积分

Gauss-Chebyshev求积公式

代数精确度2n-1

$W=\frac{1}{\sqrt{1-x^2}}$计算${\int }_{-1}^{1}f(x)(1-x^2{)}^{-\frac{1}{2}}dx\simeq \frac{\pi }{n}\sum f(cos\frac{2i-1}{2n}\pi )$