计算物理复习

计算物理复习

写在前面,本文纯属是自己复习的时候整理的,其中可能会有错误,仅供参考。
大佬勿喷,谢谢。

同时,本文可能会边更新,边发布,所以请各位阅读的时候注意。

最后,希望计算物理能考个好成绩吧。加油

插值

插值是根据已有的数据,拟合出一个过所有已知点的曲线,拟合出这种曲线就相当于在已知的数据点之间进行添加数据点。
插值的一个比较明显的特征就是:

插值过所有的已知数据点。

插值的这个特点意味着,如果使用插值法来拟合的话,那么所给的数据点必须全部都是正确的,不存在误差。

拉格朗日插值

已知 f ( x ) f(x) f(x) n + 1 n+1 n+1个任意分布的结点 x = x 0 , x 1 , … , x n x=x_0,x_1,\dots,x_n x=x0,x1,,xn上的值为 y = y 0 , y 1 , … , y n y=y_0,y_1,\dots,y_n y=y0,y1,,yn

构造拉格朗日多项式 P n ( x ) P_n(x) Pn(x),使其在 n + 1 n+1 n+1个点上和 f ( x ) f(x) f(x)有相同的值。显然,受到所给数据个数的限制, P n ( x ) P_n(x) Pn(x)是一个项数不超过 n n n的多项式。
P n ( x ) = L 0 ( x ) y 0 + L 1 ( x ) y 1 + ⋯ + L n ( x ) y n (1) P_n(x) = L_0(x)y_0 + L_1(x)y_1+\cdots+L_n(x)y_n\tag{1} Pn(x)=L0(x)y0+L1(x)y1++Ln(x)yn(1)
其中的, L 0 ( x ) , L 1 ( x ) , … , L n ( x ) L_0(x),L_1(x),\dots,L_n(x) L0(x),L1(x),,Ln(x)是待定多项式(称为 插值基函数

这里可以注意到,在求解拉格朗日多项式的时候我们主要就是求解 L n ( x ) L_n(x) Ln(x)的通项。

L i ( x j ) = { 0 , j ≠ i 1 , j = i (2) L_i(x_j) = \begin{cases} 0,\text{$j \neq i$}\\ 1,\text{$j = i$}\\ \end{cases}\tag{2} Li(xj)={0,j=i1,j=i(2)

满足式(2),则在 x 0 , x 1 , … , x n x_0,x_1,\dots,x_n x0,x1,,xn处, P n ( x ) = f ( x ) P_n(x)= f(x) Pn(x)=f(x)

由于 x = x 1 , x 2 , … , x n x=x_1,x_2,\dots,x_n x=x1,x2,,xn处, L 0 ( x ) = 0 L_0(x) = 0 L0(x)=0,而在 x = x 0 x=x_0 x=x0处, L 0 ( x ) = 1 L_0(x) = 1 L0(x)=1,所以
L 0 ( x ) = C ( x − x 1 ) ( x − x 2 ) ⋯ ( x − x n ) (3) L_0(x) = C(x-x_1)(x-x_2)\cdots(x-x_n)\tag{3} L0(x)=C(xx1)(xx2)(xxn)(3)
式(3)满足,只有当 x = x 0 x=x_0 x=x0的时候, L 0 ≠ 0 L_0\neq 0 L0=0,其他情况 L 0 = 0 L_0 = 0 L0=0

x = x 0 x=x_0 x=x0带入得到,
C = 1 ( x 0 − x 1 ) ( x 0 − x 2 ) ⋯ ( x 0 − x n ) (4) C = \frac{1}{(x_0-x_1)(x_0-x_2)\cdots(x_0-x_n)}\tag{4} C=(x0x1)(x0x2)(x0xn)1(4)
将式(4)回带到式(3)中,进而得到,
L 0 ( x ) = ( x − x 1 ) ( x − x 2 ) ⋯ ( x − x n ) ( x 0 − x 1 ) ( x 0 − x 2 ) ⋯ ( x 0 − x n ) = ∏ j = 1 n ( x − x j x 0 − x j ) (5) \begin{aligned} L_0(x) &= \frac{(x-x_1)(x-x_2)\cdots(x-x_n)}{(x_0-x_1)(x_0-x_2)\cdots(x_0-x_n)}\\ &=\prod_{j=1}^{n}\Big(\frac{x-x_j}{x_0-x_j}\Big) \end{aligned}\tag{5} L0(x)=(x0x1)(x0x2)(x0xn)(xx1)(xx2)(xxn)=j=1n(x0xjxxj)(5)
很容易从式(5)中的递推到所有的拉格朗日多项式 L i ( x ) L_i(x) Li(x),如下
L i = ∏ j = 0 j ≠ i n ( x − x j x i − x j ) (6) L_i = \prod_{j=0\\j\neq i}^{n}\Big(\frac{x-x_j}{x_i-x_j}\Big) \tag{6} Li=j=0j=in(xixjxxj)(6)
将式(6)回带到式(1)中得到:
P n ( x ) = ∑ i = 0 n L i ( x ) y i = ∑ i = 0 n { f i ( x ) ∏ j = 0 j ≠ i n ( x − x j x i − x j ) } (7) P_n(x) = \sum_{i=0}^{n}L_i(x)y_i = \sum_{i=0}^{n}\Big\{f_i(x)\prod_{j=0\\j\neq i}^{n}\Big(\frac{x-x_j}{x_i-x_j}\Big)\Big\}\tag{7} Pn(x)=i=0nLi(x)yi=i=0n{fi(x)j=0j=in(xixjxxj)}(7)

一次插值和二次插值

一次插值就是当 n = 1 n=1 n=1的时候,二次插值就是 n = 2 n=2 n=2的时候,这个直接套用上面的公式不就行了吗?为什么还非要拿出来单独讲?

n等于几就是几次插值。

余项

R n ( x ) = f ( x ) − P n ( x ) (8) R_n(x) = f(x) - P_n(x)\tag{8} Rn(x)=f(x)Pn(x)(8)

表示插值的精度,就是 P ( x ) P(x) P(x) f ( x ) f(x) f(x)的逼近程度。
R n ( x ) = f ( n + 1 ) ( ξ ) ( n + 1 ) ! ∏ i = 0 n ( x − x i ) , ξ ∈ ( a , b ) R_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^n(x-x_i),\xi\in(a,b) Rn(x)=(n+1)!f(n+1)(ξ)i=0n(xxi),ξ(a,b)
其中, a = m i n { x i , x } , b = m a x { x i , x } a = min\{x_i,x\},b = max\{x_i,x\} a=min{xi,x},b=max{xi,x}

上面的 f ( ξ ) f(\xi) f(ξ)是要拟合的函数,这里有点不太明白,既然已经知道了原本的函数表达式了,为什么还要进行拟合?当x落在插值中部的时候, ∣ ∏ ( i = 0 ) n ( x − x i ) ∣ \Big|\prod_{(i=0)}^{n}(x-x_i)\Big| (i=0)n(xxi)比较小。

余项这里,考试不怎么考,现在就是稍微提一下,以后如果用到了的话,再查也来得及。

例题

x x x1-12
f ( x ) f(x) f(x)0-34

P n ( x ) = ∑ i = 0 n L i ( x ) f ( x i ) = ∑ i = 0 n f ( x i ) ∏ j = 0 , j ≠ i n { x − x j x i − x j } \begin{aligned} P_n(x) &= \sum_{i=0}^{n}L_i(x)f(x_i)\\ &=\sum_{i=0}^{n}f(x_i)\prod_{j=0,j\neq i}^{n}\Big\{ \frac{x - x_j}{x_i - x_j} \Big\} \end{aligned} Pn(x)=i=0nLi(x)f(xi)=i=0nf(xi)j=0,j=in{xixjxxj}

x 0 = 1 , x 1 = − 1 , x 2 = 2 x_0 = 1,x_1 = -1,x_2 = 2 x0=1,x1=1,x2=2, f ( x 0 ) = 0 , f ( x 1 ) = − 3 , f ( x 2 ) = 4 f(x_0) = 0,f(x_1) = -3,f(x_2) = 4 f(x0)=0,f(x1)=3,f(x2)=4

同时根据给定的数据可以知道 n = 2 n=2 n=2

将这些数带到上面的式子中得到
P 2 ( x ) = 5 x 2 + 9 x − 14 6 P_2(x) = \frac{5x^2+9x-14}{6} P2(x)=65x2+9x14
余项:
R 2 ( x ) = f 3 ( ξ ) 3 ! ( x − 1 ) ( x + 1 ) ( x − 2 ) , ξ ∈ ( − 1 , 2 ) R_2(x) = \frac{f^3(\xi)}{3!}(x-1)(x+1)(x-2),\xi \in (-1,2) R2(x)=3!f3(ξ)(x1)(x+1)(x2),ξ(1,2)

解题技巧:

  1. 先看n等于几
  2. 将给定的x和f(x)标号,最好下标从0开始
  3. 计算 L i ( x ) L_i(x) Li(x) f ( x i ) = 0 f(x_i) = 0 f(xi)=0的话,就不用算这一项了
  4. L i ( x ) , f ( x i ) L_i(x),f(x_i) Li(x),f(xi)回带到 P n ( x ) P_n(x) Pn(x)

小结

这个拉格朗日多项式写程序最重要的就是最后 L i ( x ) L_i(x) Li(x)的递推公式,而且仔细观察可以发现 L i ( x ) L_i(x) Li(x)的递推公式和 f ( x ) f(x) f(x)无关,只和 x x x有关,这个应该有一些道道,

如果是理解这个拉格朗日怎么来的话,主要就是弄明白式(3)就行了。

均差法

起因

改进拉格朗日插值多项式的缺点(结点增加时,整个公式必须重新计算)。

定义均差

已知函数 f ( x ) f(x) f(x),和函数曲线上的一些点
f ( x k ) − f ( x j ) x k − x j , k ≠ j (1) \frac{f(x_k) - f(x_j)}{x_k-x_j},k\neq j\tag{1} xkxjf(xk)f(xj),k=j(1)
式(1)叫做 f ( x ) f(x) f(x)的一阶均差,记作 f ( x k , x j ) f(x_k,x_j) f(xk,xj)表示 f ( x ) f(x) f(x)在区间 [ x k , x j ] [x_k,x_j] [xk,xj]的平均变化率

f ( x ) f(x) f(x)的二阶均差如下
f ( x i , x j , x k ) = f ( x j , x k ) − f ( x j , x i ) x k − x i = f ( x j ) − f ( x k ) x j − x k − f ( x j ) − f ( x i ) x j − x i x k − x i (2) \begin{aligned} f(x_i,x_j,x_k) &= \frac{f(x_j,x_k) - f(x_j,x_i)}{x_k-x_i}\\ &=\frac{\frac{f(x_j) - f(x_k)}{x_j - x_k} - \frac{f(x_j) - f(x_i)}{x_j - x_i}}{x_k - x_i} \end{aligned}\tag{2} f(xi,xj,xk)=xkxif(xj,xk)f(xj,xi)=xkxixjxkf(xj)f(xk)xjxif(xj)f(xi)(2)
f ( x ) f(x) f(x)的n阶均差如下,
f ( x 0 , x 1 , ⋯   , x n ) = f ( x 0 , x 1 , ⋯   , x n − 1 ) − f ( x 1 , x 2 , ⋯   , x n ) x 0 − x n (3) f(x_0,x_1,\cdots,x_n) = \frac{f(x_0,x_1,\cdots,x_{n-1}) - f(x_1,x_2,\cdots,x_n)}{x_0 - x_n}\tag{3} f(x0,x1,,xn)=x0xnf(x0,x1,,xn1)f(x1,x2,,xn)(3)
通过观察可以发现,求解 f ( x ) f(x) f(x)的n阶均差的时候,可以利用 f ( x ) f(x) f(x)的n-1阶均差快速计算。

均差表

x i x_i xi f ( x i ) f(x_i) f(xi)一阶均差二阶均差三阶均差 ⋯ \cdots n阶均差
x 0 x_0 x0 f ( x 0 ) f(x_0) f(x0)
x 1 x_1 x1 f ( x 1 ) f(x_1) f(x1) f ( x 0 , x 1 ) f(x_0,x_1) f(x0,x1)
x 2 x_2 x2 f ( x 2 ) f(x_2) f(x2) f ( x 1 , x 2 ) f(x_1,x_2) f(x1,x2) f ( x 0 , x 1 , x 2 ) f(x_0,x_1,x_2) f(x0,x1,x2)
x 3 x_3 x3 f ( x 3 ) f(x_3) f(x3) f ( x 2 , x 3 ) f(x_2,x_3) f(x2,x3) f ( x 1 , x 2 , x 3 ) f(x_1,x_2,x_3) f(x1,x2,x3) f ( x 0 , x 1 , x 2 , x 3 ) f(x_0,x_1,x_2,x_3) f(x0,x1,x2,x3)
⋮ \vdots ⋮ \vdots ⋮ \vdots ⋮ \vdots ⋮ \vdots ⋯ \cdots
x n x_n xn f ( x n ) f(x_n) f(xn) f ( x n − 1 , x n ) f(x_{n-1},x_n) f(xn1,xn) f ( x n − 2 , x n − 1 , x n ) f(x_{n-2},x_{n-1},x_{n}) f(xn2,xn1,xn) f ( x n − 3 , x n − 2 , x n − 1 , x n ) f(x_{n-3},x_{n-2},x_{n-1},x_{n}) f(xn3,xn2,xn1,xn) ⋯ \cdots f ( x 0 , x 1 , … , x n ) f(x_0,x_1,\dots,x_n) f(x0,x1,,xn)

这个均差表非常的好用,基本上就说明了均差法的全部计算原理。

  • 使用n个节点,就说明使用了n-1次插值
  • 均差表中的 x x x节点好像都必须是连续的。
  • 再有增加的节点的话,增加在最下面计算效率最高

上面其实还应该有一个最最主要的性质,就是每个均差的计算都可以通过上一阶的均差计算出来,但是由于需要除以x的差值,这个差值不知道怎么描述的好,所以就被我删掉了。

这里在唠叨唠叨,就是当如果计算均差的时候,你写的记号满足顺序的规则 f ( x 0 , x 1 , x 2 , ⋯   , x n ) f(x_0,x_1,x_2,\cdots,x_n) f(x0,x1,x2,,xn),或者逆序 f ( x n , x n − 1 , x n − 2 , ⋯   , x 0 ) f(x_n,x_{n-1},x_{n-2},\cdots,x_0) f(xn,xn1,xn2,,x0)的规则(总之就是不能随便的写 f ( x 2 , x 1 , x n , x 5 , ⋯   ) f(x_2,x_1,x_n,x_5,\cdots) f(x2,x1,xn,x5,)),那么计算均差的通项就能固定下来了
f ( x 0 , x 1 , x 2 , ⋯   , x n ) = f ( x 0 , x 1 , x 2 , ⋯   , x n − 1 ) − f ( x 1 , x 2 , x 3 , ⋯   , x n ) x 0 − x n f(x_0,x_1,x_2,\cdots,x_n) = \frac{f(x_0,x_1,x_2,\cdots,x_{n-1}) - f(x_1,x_2,x_3,\cdots,x_n)}{x_0-x_n} f(x0,x1,x2,,xn)=x0xnf(x0,x1,x2,,xn1)f(x1,x2,x3,,xn)
或者
f ( x n , x n − 1 , x n − 2 , ⋯   , x 0 ) = f ( x 1 , x 2 , x 3 , ⋯   , x n ) − f ( x 0 , x 1 , x 2 , ⋯   , x n − 1 ) x n − x 0 f(x_n,x_{n-1},x_{n-2},\cdots,x_0) = \frac{f(x_1,x_2,x_3,\cdots,x_{n}) - f(x_0,x_1,x_2,\cdots,x_{n-1})}{x_n-x_0} f(xn,xn1,xn2,,x0)=xnx0f(x1,x2,x3,,xn)f(x0,x1,x2,,xn1)
总之就是用均差值相减,然后除以两边相减。

均差多项式

P n ( 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 1 , x 2 , ⋯   , x n ) ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x n − 1 ) = P n − 1 + f ( x 0 , x 1 , x 2 , ⋯   , x n ) ( x − x 0 ) ( x − x 1 ) ( x − x 2 ) ⋯ ( x − x n − 1 ) P_n(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)+\cdots+f(x_0,x_1,x_2,\cdots,x_n)(x-x_0)(x-x_1)\cdots(x-x_{n-1}) \\=P_{n-1} + f(x_0,x_1,x_2,\cdots,x_n)(x-x_0)(x-x_1)(x-x_2)\cdots(x-x_{n-1}) Pn(x)=f(x0)+f(x0,x1)(xx0)+f(x0,x1,x2)(xx0)(xx1)++f(x0,x1,x2,,xn)(xx0)(xx1)(xxn1)=Pn1+f(x0,x1,x2,,xn)(xx0)(xx1)(xx2)(xxn1)

这里,如果之前计算出来了均差表的话,这里就是非常好计算的,当然这里用到的均差只是上方表格的对角线上的均差,但是我们想要计算对角线上的均差的话,就必须计算出前面全部的均差值。

例题1

在这里插入图片描述

解题技巧:

  1. 先设出来 x 0 , x 1 , x 2 , ⋯   , x n x_0,x_1,x_2,\cdots,x_n x0,x1,x2,,xn f ( x 0 ) , f ( x 1 ) , f ( x 2 ) , ⋯   , f ( x n ) , f(x_0),f(x_1),f(x_2),\cdots,f(x_n), f(x0),f(x1),f(x2),,f(xn),
  2. 然后开始写均差表
  3. 写出均差表之后,开始写多项式的通项

例题2

在这里插入图片描述

在这里插入图片描述

解题技巧:

  1. 使用 n n n个节点就说明使用的是 n − 1 n-1 n1次插值
  2. n n n还在拟合函数 P n ( x ) P_n(x) Pn(x)的下标那里出现过

差分

起因

自变量是等距结点时, 函数的变化率与自变量区间无关,使用有限差分建立插值公式。

定义

课件上的定义:
在这里插入图片描述
在这里插入图片描述

不得不吐槽一下,这个定义,真的是老太太的裹脚布又臭又长!!!

从上面的定义中完全看不出来二阶差分具体怎么来的,当时我就是上课的时候,看到了这里,结果一时脑子没转过来,导致后面没跟上,下课又自己看的。

下面仔细来说说这个定义:

首先看通项:
△ n f k = △ ( △ n − 1 f k ) = [ △ n − 1 f ( x k + h ) − △ n − 1 f ( x k ) ] \bigtriangleup^nf_k = \bigtriangleup(\bigtriangleup^{n-1}f_k) = [\bigtriangleup^{n-1}f(x_k+h) - \bigtriangleup^{n-1} f(x_k)] nfk=(n1fk)=[n1f(xk+h)n1f(xk)]
这里就不会再多写上一步?
△ n f k = △ ( △ n − 1 f k ) = [ △ n − 1 f ( x k + 1 ) − △ n − 1 f ( x k ) ] \bigtriangleup^nf_k = \bigtriangleup(\bigtriangleup^{n-1}f_k) = [\bigtriangleup^{n-1}f(x_{k+1}) - \bigtriangleup^{n-1} f(x_k)] nfk=(n1fk)=[n1f(xk+1)n1f(xk)]
这样写的话,一下子就清楚明了了啊!!当前的差分可以由之前求出来的差分直接相减求出来。

明确指出一个点:

f k = f ( x k ) f_k = f(x_k) fk=f(xk)

上面一阶中心插值基本不考,主要就是向前差分和向后差分。

并且从上面的定义中也能看出来,n阶差分可以由n-1阶差分求出来,所以之前算出来的差分,在后面求差分的时候能直接拿来用。

并且注意,向前差分和向后差分的不同,虽然两者看上去非常的像,但是可以看出来向后差分是知道了前面的,然后 k k k是慢慢递增的,同理,向前差分是知道了后面的然后 k k k是慢慢递减的。这一点在后面的差分表中可以非常的明显看出来。

差分和均差的关系

f ( x 0 , x 1 ) = f 1 − f 0 h = △ f 0 h f ( x 0 , x 1 , x 2 ) = f ( x 0 , x 1 ) − f ( x 1 , x 2 ) x 0 − x 2 = 1 2 h ( △ f 0 h − △ f 1 h ) = △ 2 f 0 2 h 2 f ( x 0 , x 1 , x 2 , ⋯   , x n ) = △ n f 0 n ! h n \begin{aligned} f(x_0,x_1)& = \frac{f_1-f_0}{h} = \frac{\bigtriangleup f_0}{h}\\ f(x_0,x_1,x_2)& = \frac{f(x_0,x_1) - f(x_1,x_2)}{x_0-x_2} = \frac{1}{2h}(\frac{\bigtriangleup f_0}{h} - \frac{\bigtriangleup f_1}{h}) = \frac{\bigtriangleup^2f_0}{2h^2}\\ f(x_0,x_1,x_2,\cdots,x_n) &= \frac{\bigtriangleup^nf_0}{n!h^n} \end{aligned} f(x0,x1)f(x0,x1,x2)f(x0,x1,x2,,xn)=hf1f0=hf0=x0x2f(x0,x1)f(x1,x2)=2h1(hf0hf1)=2h22f0=n!hnnf0

差分表

在这里插入图片描述

从上面的差分表是标准的差分表,其实我们一般写成下面的形式

x i x_i xi f ( x ) f(x) f(x)一阶差分二阶差分三阶差分
x 0 x_0 x0 f 0 f_0 f0
x 1 x_1 x1 f 1 f_1 f1 △ f 0 \bigtriangleup f_0 f0
x 2 x_2 x2 f 2 f_2 f2 △ f 1 \bigtriangleup f_1 f1 △ 2 f 0 \bigtriangleup ^2 f_0 2f0
x 3 x_3 x3 f 3 f_3 f3 △ f 2 \bigtriangleup f_2 f2 △ 2 f 1 \bigtriangleup ^2 f_1 2f1 △ 3 f 0 \bigtriangleup^3 f_0 3f0

课件上不知道为什么,从 x 0 x_0 x0后面紧接的跟着的是 f 1 f_1 f1这和之前我们定义的 f k = f ( x k ) f_k = f(x_k) fk=f(xk)不一致,所以我在上面的表格中改成了 f 0 f_0 f0并且后面的所有的 f f f都将下标减一。

这里由于课件老师给的是pdf版本的,个人无法修改,就没有改,各位(也不知道有没有人)看的时候,注意一下。

这里只是列举到了三阶差分,其实考试最多也就到四阶差分,这里就先列出这些了,然后我们可以看出来,如果是向前差分的话,其实就是主对角线上的那些数据,就是 f 0 , △ f 0 , △ 2 f 0 , △ 3 f 0 f_0,\bigtriangleup f_0,\bigtriangleup^2 f_0,\bigtriangleup^3 f_0 f0,f0,2f0,3f0,如果是向后差分的话,就是最底下的那一行,感觉有点怪其实就是这样的,和课件上的对比也会发现,其实主对角线上的数据和课件上的红线标出出来的数据一样,最后一行的数据和课件上的蓝色标出来的数据是一样的。

牛顿向前插值

通项公式

均差法插值多项式
P n ( 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 1 , x 2 , ⋯   , x n ) ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x n − 1 ) P_n(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)+\cdots+f(x_0,x_1,x_2,\cdots,x_n)(x-x_0)(x-x_1)\cdots(x-x_{n-1}) Pn(x)=f(x0)+f(x0,x1)(xx0)+f(x0,x1,x2)(xx0)(xx1)++f(x0,x1,x2,,xn)(xx0)(xx1)(xxn1)
其实牛顿向前插值就是 x k + 1 − x k = h ( 常 数 ) x_{k+1}-x_k = h(常数) xk+1xk=h()的均差插值多项式,所以
P n ( x ) = f ( x 0 ) + ( x − x 0 ) △ f ( x 0 ) h + ( x − x 0 ) ( x − x 1 ) △ 2 f ( x 0 ) 2 ! h 2 + ⋯ + ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x n − 1 ) △ n f ( x 0 ) n ! h n P_n(x) = f(x_0) + (x-x_0)\frac{\bigtriangleup f(x_0)}{h}+(x-x_0)(x-x_1)\frac{\bigtriangleup^2 f(x_0)}{2!h^2}+\cdots+(x-x_0)(x-x_1)\cdots(x-x_{n-1})\frac{\bigtriangleup^nf(x_0)}{n!h^n} Pn(x)=f(x0)+(xx0)hf(x0)+(xx0)(xx1)2!h22f(x0)++(xx0)(xx1)(xxn1)n!hnnf(x0)
再令 x = x 0 + s h x=x_0 + sh x=x0+sh 0 ≤ s ≤ 1 0\leq s \leq 1 0s1,则 x − x i = ( s − i ) h x-x_i = (s-i)h xxi=(si)h,所以

在这里插入图片描述

牛顿向后插值

牛顿向后插值多项式和牛顿向前插值多项式类似,最终的插值公式如下

在这里插入图片描述

例题

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

解题技巧:

  1. 首先标定 x 0 , x 1 , x 2 , ⋯   , x n x_0,x_1,x_2,\cdots,x_n x0,x1,x2,,xn f ( x 0 ) , f ( x 1 ) , ⋯   , f ( x n ) f(x_0),f(x_1),\cdots,f(x_n) f(x0),f(x1),,f(xn)
  2. 列出差分表
  3. 看看是要用向前插值还是向后插值
  4. 选取基点
    1. 向前插值:从基点开始画一条平行于主对角线的直线,直线经过的数据就是我们需要带入到差值公式中去的差分
    2. 向后插值:选取基点所在的那行,就是我们需要带入到插值公式中的差分
  5. 然后将差分表中的数据带入到插值公式中

曲线平滑和数据拟合

插值得到的拟合曲线,经过所有给出的点。拟合得到的函数不一定经过所有的点。

在这里插入图片描述

做数据拟合首先需要有一组基函数。

注意上面出现的几个常量:

  • n:n对数据点
  • m:m个基函数

这两个常量在下面的时候需要用到。

上面的课件中又有一点我感觉不对的地方!!!!

为什么 c i ( i = 1 , 2 , … , m ) c_i(i=1,2,\dots,m) ci(i=1,2,,m)这里 i i i不是从0开始的?!!!

线性拟合的最小二乘法

最小二乘法记住一个最最主要的公式:正则公式
P C = Q PC = Q PC=Q

推导正则公式

  • 假设我们已经确定了参数 c i c_i ci,那么对于给定的 x k x_k xk,我们就能计算出线性模型的 y k y_k yk的估计值 y k ^ ,   ( k = 0 , 1 , 2 , ⋯   , n ) \hat{y_k},\ (k=0,1,2,\cdots,n) yk^, (k=0,1,2,,n)

  • 首先我们拟合出来的数据 y ^ k \hat{y}_k y^k和给定的值 y k y_k yk不完全相同,他们之间的差叫做残差,记作

    e k = y k ^ − y k = ∑ j = 0 n c j ϕ j ( x k ) − y k , k = 0 , 1 , 2 , ⋯   , n e_k = \hat{y_k} - y_k=\sum_{j=0}^{n}c_j\phi_j(x_k) - y_k,k =0,1,2,\cdots,n ek=yk^yk=j=0ncjϕj(xk)yk,k=0,1,2,,n

  • 当数据点给定之后( x k , y k x_k,y_k xk,yk给定),上式中残差的大小只取决于 c i c_i ci的取值,所以残差的大小就是衡量参数及拟合程度好坏的标志

    这里插一嘴,我之前看过一些机器学习上的预测类的问题,上面好像说过有一种叫“过拟合”的问题,就是对于已有的数据点并不是符合的越好这个模型就越好,需要看训练出来的模型对未知数据的拟合程度,如果对于未知数据的拟合程度非常的好,那才是真正的好。

    和这里的要求好像不太一样,也不知道这其中的道理是什么。

    还是说这里的最小二乘法只关注已有的数据,而机器学习中的预测问题的关注点在未知数据上。

  • 残差的平方和最小
    r = ∑ k = 0 n e k 2 r = \sum_{k=0}^ne_k^2 r=k=0nek2

  • 按照残差和最小的原则求解,得到 c i c_i ci的过程称为 最小二乘法

  • 一般用残差和的根表示
    ϵ = r = ∑ k = 0 n ( y ^ k − y k ) 2 \epsilon = \sqrt{r} = \sqrt{\sum_{k=0}^n(\hat{y}_k - y_k)^2} ϵ=r =k=0n(y^kyk)2

    求解残差平方和最小

    ∂ r ∂ c k = ∂ ∂ c k [ ∑ k = 0 n e k 2 ] = 0 ( k = 0 , 1 , 2 , ⋯   , n ) \frac{\partial r}{\partial c_k} = \frac{\partial }{\partial c_k}\Big[\sum_{k=0}^ne_k^2\Big] = 0\quad(k=0,1,2,\cdots,n) ckr=ck[k=0nek2]=0(k=0,1,2,,n)

    展开上式
    ∂ r ∂ c k = ∂ ∂ c k [ ∑ i = 0 ( y ^ k − y k ) 2 ] = ∂ ∂ c k [ ∑ i = 0 n ( ∑ j = 0 m c j ϕ j ( x i ) − y i ) 2 ] \frac{\partial r}{\partial c_k} = \frac{\partial}{\partial c_k}\Big[\sum_{i=0}(\hat{y}_k - y_k)^2\Big]=\frac{\partial}{\partial c_k}\Big[\sum_{i=0}^n\Big(\sum_{j=0}^mc_j\phi_j(x_i) - y_i\Big)^2\Big] ckr=ck[i=0(y^kyk)2]=ck[i=0n(j=0mcjϕj(xi)yi)2]
    取出求导项
    ∑ j = 0 m [ ∑ i = 0 n ϕ k ( x i ) ϕ j ( x i ) ] c j = ∑ i = 0 n ϕ k ( x i ) y i \sum_{j=0}^m\Big[ \sum_{i=0}^{n}\phi_k(x_i)\phi_j(x_i) \Big]c_j = \sum_{i=0}^{n}\phi_k(x_i)y_i j=0m[i=0nϕk(xi)ϕj(xi)]cj=i=0nϕk(xi)yi

    这里的具体求导,不是很难,可以自己推导,只要把求偏导符号提到两个求和之间就行了。然后按照复合函数求导。

上式可以看成一个矩阵( n × m n\times m n×m)乘以一个向量( m × 1 m\times 1 m×1)等于另一个向量( n × 1 n\times 1 n×1),即正则方程

P C = Q PC=Q PC=Q

下面开始具体化 P P P矩阵

P = [ p 00 p 01 ⋯ p 0 m p 10 p 11 ⋯ p 1 m ⋮ ⋮ ⋯ ⋮ p n 0 p n 1 ⋯ p n m ] P = \left[ \begin{matrix} p_{00}&p_{01}&\cdots&p_{0m}\\ p_{10}&p_{11}&\cdots&p_{1m}\\ \vdots&\vdots&\cdots&\vdots\\ p_{n0}&p_{n1}&\cdots&p_{nm} \end{matrix} \right] P=p00p10pn0p01p11pn1p0mp1mpnm

P中的每一个元素的值为

p i j = ∑ k = 0 ϕ i ( x k ) ϕ j ( x k ) ( k = 0 , 1 , 2 , ⋯   , n ) p_{ij} = \sum_{k=0}\phi_i(x_k)\phi_j(x_k)\quad(k=0,1,2,\cdots,n) pij=k=0ϕi(xk)ϕj(xk)(k=0,1,2,,n)

Q的具体形式

Q = [ q 0 q 1 q 2 ⋮ q n ] Q = \left[ \begin{matrix} q_0\\ q_1\\ q_2\\ \vdots\\ q_n \end{matrix} \right] Q=q0q1q2qn

Q中的每个元素的值为

q k = ∑ i = 0 n ϕ k ( x i ) y i q_k = \sum_{i=0}^n\phi_k(x_i)y_i qk=i=0nϕk(xi)yi

C的具体形式

C = [ c 0 c 1 ⋮ c m ] C = \left[ \begin{matrix} c_0\\ c_1\\ \vdots\\ c_m \end{matrix} \right] C=c0c1cm

将正则方程写开

[ p 00 p 01 ⋯ p 0 m p 10 p 11 ⋯ p 1 m p 20 p 21 ⋯ p 2 m ⋮ ⋮ ⋯ ⋮ p n 0 p n 1 ⋯ p n m ] [ c 0 c 1 ⋮ c m ] = [ q 0 q 1 q 2 ⋮ q n ] % P矩阵 \left[ \begin{matrix} p_{00}&p_{01}&\cdots&p_{0m}\\ p_{10}&p_{11}&\cdots&p_{1m}\\ p_{20}&p_{21}&\cdots&p_{2m}\\ \vdots&\vdots&\cdots&\vdots\\ p_{n0}&p_{n1}&\cdots&p_{nm}\\ \end{matrix} \right] % C向量 \left[ \begin{matrix} c_0\\ c_1\\ \vdots\\ c_m \end{matrix} \right]= % Q向量 \left[ \begin{matrix} q_0\\ q_1\\ q_2\\ \vdots\\ q_n \end{matrix} \right] p00p10p20pn0p01p11p21pn1p0mp1mp2mpnmc0c1cm=q0q1q2qn

注意上面可以看到C向量的“高度”和P矩阵的“高度”并不相同,我是故意这样写的,之所以这么做是为了明显的强调C向量的行数和P矩阵的列数相同,而Q矩阵的行数和P矩阵的行数相同。

主要是现在课件上没有将这个正则方程写开,导致我在学这一段的时候,写代码的时候,总是因为维度问题出bug,这里写的明明白白的,将P矩阵的行数和列数分析的明明白白的,老师上课总是讲不细致,下课还是得自己看课件,自己推导。

仔细分析P矩阵的行和列都代表着什么。P矩阵的列数为n代表着n元方程组(和Q的行数相同),P矩阵的列数为m,m代表着基函数的个数(和C矩阵的行数相同)。

----------------------------------------------------------------------------------------我是分割线--------------------------------------------------------------------------------

下面只是一些,不重要的推导,我也不知道我在推导什么,直接跳过吧。

按照矩阵相乘的基本规律( a i j = ∑ j = 0 n b i j ∑ k = 0 n c k i a_{ij} = \sum_{j=0}^nb_{ij}\sum_{k=0}^nc_{ki} aij=j=0nbijk=0ncki),P矩阵的每一行都要分别乘上对应位置的 C C C然后求和当做Q矩阵的对应的元素,即
∑ i = 0 , k = i m p 0 i c k 0 = q 00 \sum_{i=0,k=i}^mp_{0i}c_{k0} = q_{00} i=0,k=imp0ick0=q00
同时考虑到 C , Q C,Q C,Q都是向量,所以上式可以写成下式
∑ i = 0 m p 0 i c i = q 0 \sum_{i=0}^mp_{0i}c_{i} = q_{0} i=0mp0ici=q0
上式可以类推出
∑ j = 0 m p i j c j = q i ( i = 0 , 1 , 2 ⋯   , n ) \sum_{j=0}^mp_{ij}c_{j} = q_{i}\quad(i=0,1,2\cdots,n) j=0mpijcj=qi(i=0,1,2,n)
----------------------------------------------------------------------------------------我是分割线--------------------------------------------------------------------------------

一次线性拟合的P矩阵和Q向量

在这里插入图片描述

个人感觉最好还是把上图的 i i i改成 k k k,因为 i , j i,j i,j一般代表的是行数和列数。还有就是上面的i是从1开始的,但是我习惯从0开始所以上面的P矩阵需要改写
P = [ n + 1 ∑ i = 0 n ∑ i = 0 n x i ∑ i = 1 n x 2 ] P = \left[ \begin{matrix} n+1&\sum_{i=0}^{n}\\ \sum_{i=0}^nx_i&\sum_{i=1}^nx^2 \end{matrix} \right] P=[n+1i=0nxii=0ni=1nx2]
虽然,官方是说, p 00 = n p_{00} = n p00=n但是写代码写多了,数组的下标是从0开始的这是我们公认的,再加上有点强迫症,所以这里必须改成从0开始,我才舒服。。。。

例题

1 直接计算

在这里插入图片描述

中间计算矩阵元的的操作按照上面的推导就能得到了,这里就不贴出来了,太大了。

在这里插入图片描述

这里需要用计算器计算。

2 需要进行变量替换

在这里插入图片描述在这里插入图片描述

三点平滑法

原理

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

计算

三点平滑最简单了,感觉就是求平均值,就是在两端不一样。

五点平滑和三点平滑差不多,就是三点平滑求得是三个点的平均值,五点平滑求得是五个点的平均值

首先给定一组数据

x x x x 0 x_0 x0 x 1 x_1 x1 x 2 x_2 x2 x 3 x_3 x3
y y y y 0 y_0 y0 y 1 y_1 y1 y 2 y_2 y2 y 3 y_3 y3

下面对这组数据进行三点平滑

x x x x 0 x_0 x0 x 1 x_1 x1 x 2 x_2 x2 x 3 x_3 x3
y y y 1 3 [ y 0 + y 1 + y 2 ] − 1 2 [ y 2 − y 0 ] \frac{1}{3}\Big[y_0+y_1+y_2\Big] - \frac{1}{2}\Big[y_2-y_0\Big] 31[y0+y1+y2]21[y2y0] 1 3 [ y 0 + y 1 + y 2 ] \frac{1}{3}\Big[y_0+y_1+y_2\Big] 31[y0+y1+y2] 1 3 [ y 1 + y 2 + y 3 ] \frac{1}{3}\Big[y_1+y_2+y_3\Big] 31[y1+y2+y3] 1 3 [ y 1 + y 2 + y 3 ] + 1 2 [ y 3 − y 1 ] \frac{1}{3}\Big[y_1+y_2+y_3\Big] +\frac{1}{2}\Big[y_3-y_1\Big] 31[y1+y2+y3]+21[y3y1]

仔细观察上面表格,不难得到三点平滑的规律,把整个数据表分成三份

在这里插入图片描述

五点平滑法

五点平滑,应该不考。


在这里插入图片描述

局部最优化

在这里插入图片描述

黄金分割法

在这里插入图片描述

所谓黄金分割法就是将区间按照黄金比例分割,两边的那部分长度为区间长度乘以w,剩下的就是中间的,也是寻找最小值所在的区间,然后选取出新的区间,重复上述操作,直到达到收敛标准

这里的收敛标准着实有点让人搞不清,。。。。。

连续抛物线法

在这里插入图片描述

三点法

在这里插入图片描述

所谓的三点法就是将给定的区间等分成三份,计算中间两个点的函数值,再加上两边的两个数据点,现在一共是四个点,然后找出这四个点中最小的那个点,选取包含最小值的区间为新的区间,对新的区间三等分,再次进行上述操作,直到前后两次的误差在允许的范围之内

个人感觉最后的这种判断标准就是在扯淡,如果这个曲线是一个比较平滑的下降呢?每次找到的点之间的差值都比较小,那么直接就停止了吗?

单行法

在这里插入图片描述

广义搜索法

在这里插入图片描述

最速下降法

在这里插入图片描述在这里插入图片描述

注意:上面的课件错了,红字部分是"代入 f ( x , y ) f(x,y) f(x,y)"

其实所谓得最速下降法,就是数学上的梯度。

我记得高数老师曾经讲过,一个人如果想要下山最快得话,就是沿着梯度得方向下山最快。 梯度是一个矢量。让x方向沿着梯度的x分量,y方向沿着梯度的y分量确实是下降速度最快的方向。

这里面主要就是涉及到了计算步长的问题。从上面图片中已经清楚的说明了将 f ( x 2 , y 2 ) f(x_2,y_2) f(x2,y2)最小化,就能得到步长 α \alpha α了。

这里还是想要插一嘴,我们现在就是再求最小化问题,说明我们不会求最小化,这里为什么还要让我们通过 最小化 f ( x 2 , y 2 ) f(x_2,y_2) f(x2,y2) 来求出来 α \alpha α呢?

或者是 f ( x 2 , y 2 ) f(x_2,y_2) f(x2,y2)只是 α \alpha α的函数,求最小化的时候好求吗?

共轭梯度下降法

在这里插入图片描述
在这里插入图片描述

这个共轭梯度下降法,就是在之前最速下降法的基础之上有改进了一点,主要就是在 β \beta β的那一步,老师上课讲的时候,只是说这样改进可以加快收敛速度,说人话就是更快的找到最小值,具体为什么能更快的找到最小值,这个我也不清楚。

但是,记住求解步骤就好了。

例题

在这里插入图片描述

全局最优化

在这里插入图片描述
在这里插入图片描述

全局最优化比之局部最优化的一个最大的区别就是,全局最小只有一个,但是存在很多个全局极小值。全局最优化算法并不会再一股脑的向着当前点下降最快的方向移动了,而是规定了一种新的方式达到下一个新的状态。

众多全局最优化算法之中的核心区别就是产生新状态的方式,或准则不同。

模拟退火算法

模拟退火算法的基本过程

在这里插入图片描述
在这里插入图片描述

在这里插入图片描述
在这里插入图片描述
上面课件中又错了,“若 E j > E i E_j>E_i Ej>Ei,按照一定的概率接受。”课件中写成了小于了。

注意这一页的Metropolis准则是模拟退火算法的核心思想。因为它定义了算法是怎样接受新状态的。

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

模拟退火算法考试应该不会让我们手算一遍的,因为基本上没啥意义。

模拟退火算法的中产生新状态需要随机数,,但是我们用手算的话哪来的随机值?

所以这一部分应该是不会考的。

遗传算法

遗传算法的核心就是对自变量进行编码(编码方式可以有多种多样的,但是对于我们这种非专业的话,就是用二进制编码),然后将编码后的标号看做是DNA链,两个DNA链“纠缠”组合成新的下一个值,然后将新状态中好的那一部分保留下来,将坏的丢掉,继续产生下一轮新的状态。

直到产生的新状态满足我们的要求

这个考试考的可能性很大!!!

遗传算法的定义

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述
在这里插入图片描述

标准遗传算法的步骤

在这里插入图片描述

例题 1

在这里插入图片描述

在这里插入图片描述

其实看了上面的那个例题就明白了,其实就是将一个区间等分,然后再根据当前x的二进制的值转换成10进制。然后回带到原函数中。

例题2

在这里插入图片描述

例题3

在这里插入图片描述

积分

矩形积分法

将被积函数等分成n份,然后将每一块看做成小矩形,函数值(一般是用前面的那个函数值。)看做小矩形的高,区间宽度看做小矩形的宽。

这里注意,矩形法将区间等分成几块就取几个点,因为每一个区间我们只用一个点就能代表这个小矩形的高。

在这里插入图片描述

梯形积分法

梯形积分就是将一个又一个的区间看成小梯形,然后求每个梯形的面积之后将所有梯形的面积加起来。

这里将区间等分成n份,需要取n+1个点,因为每个区间我们需要知道梯形的上底和下底
在这里插入图片描述

做题的话,记住最后的公式就行了。这样快。但是并不是说我们不知道原理,我们知道原理,只不过我们选取了一中更快的方法来做题。

抛物线积分法(辛普森积分法)

推导过程

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

可以看到,辛普森法就是也是需要将区间等分成N份,只不过这次我们将两个小区间看成一个四边形,四边形的左边高,用左边的函数值表示,右边的高用右边的函数值代替,宽就是两个小区间的宽度,然后顶部用一个过这两个小区间中选取出来的三个点的抛物线代替。

这样就是计算这个区间的面积。

最后考试的话,记住每个 C i C_i Ci等于多少就行了。

例题1

在这里插入图片描述

在这里插入图片描述

高斯求积法

在这里插入图片描述

上图给出来 A k A_k Ak的求法,但是还未给出 x k x_k xk的取法。
在这里插入图片描述
在这里插入图片描述

总结一下:

  1. 高斯求积法的通项为:
    ∑ k = 0 n A k f ( x k ) \sum_{k=0}^nA_kf(x_k) k=0nAkf(xk)

  2. 首先看看积分区间是不是[-1,1],如果在[-1,1]之间进入第三步。如果不是的话,先进行换元,将积分区间改到[-1,1]之间。同时注意将被积函数进行替换。

  3. 判断是要求让用几点高斯求积,这里假设是n点

  4. 找到n阶勒让德多项式,然后找到当前阶数的勒让德多项式的所有零点

  5. 所有零点就是要取的 x k x_k xk

  6. 然后计算所有的 A k A_k Ak,这里 A k A_k Ak的表达式记住就好了
    A k = 2 ( 1 − x k 2 ) [ p n + 1 ′ ( x k ) ] 2 A_k = \frac{2}{(1-x_k^2)[p'_{n+1}(x_k)]^2} Ak=(1xk2)[pn+1(xk)]22

  7. 然后回带到公司求积的通项中得到,积分的结果。

例题

在这里插入图片描述

常用高斯积分

1点高斯积分

x 0 = 0 , A 0 = 2 x_0 = 0,A_0 = 2 x0=0,A0=2

2点高斯积分

x 0 = 1 3 , x 1 = − 1 3 x_0 = \frac{1}{\sqrt{3}},x_1 = - \frac{1}{\sqrt{3}} x0=3 1,x1=3 1

A 0 = A 1 = 1 A_0 = A_1 = 1 A0=A1=1

3点高斯积分

x 0 = − 15 5 , x 1 = 0 , x 2 = 15 5 x_0 = -\frac{\sqrt{15}}{5},x_1 = 0,x_2 = \frac{\sqrt{15}}{5} x0=515 ,x1=0,x2=515

A 0 = 5 9 , A 1 = 8 9 , A 2 = 5 9 A_0 = \frac{5}{9},A_1 = \frac{8}{9},A_2 = \frac{5}{9} A0=95,A1=98,A2=95

求解线性方程组

列主消元法

{ a 11 x 1 + a 12 x 2 + ⋯ + a 1 n x n = b 1 a 21 x 1 + a 22 x 2 + ⋯ + a 2 n x n = b 2 a 31 x 1 + a 32 x 2 + ⋯ + a 3 n x n = b 3 a 41 x 1 + a 42 x 2 + ⋯ + a 4 n x n = b 4 a 51 x 1 + a 52 x 2 + ⋯ + a 5 n x n = b 5 ⋮ + ⋮ + ⋯ + ⋮ = ⋮ a n 1 x 1 + a n 2 x 2 + ⋯ + a n n x n = b n \begin{cases} a_{11}x_1 + a_{12}x_2 + \cdots+a_{1n}x_n = b_1\\ a_{21}x_1 + a_{22}x_2 + \cdots+a_{2n}x_n = b_2\\ a_{31}x_1 + a_{32}x_2 + \cdots+a_{3n}x_n = b_3\\ a_{41}x_1 + a_{42}x_2 + \cdots+a_{4n}x_n = b_4\\ a_{51}x_1 + a_{52}x_2 + \cdots+a_{5n}x_n = b_5\\ \quad \vdots\quad+\quad\vdots\quad+ \cdots +\quad \vdots\quad =\quad \vdots\\ a_{n1}x_1 + a_{n2}x_2 + \cdots+a_{nn}x_n = b_n\\ \end{cases} a11x1+a12x2++a1nxn=b1a21x1+a22x2++a2nxn=b2a31x1+a32x2++a3nxn=b3a41x1+a42x2++a4nxn=b4a51x1+a52x2++a5nxn=b5+++=an1x1+an2x2++annxn=bn
增广矩阵
[ a 11 a 12 a 13 ⋯ a 1 n ∣ b 1 a 21 a 22 a 23 ⋯ a 2 n ∣ b 2 a 31 a 32 a 33 ⋯ a 3 n ∣ b 3 ⋮ ⋮ ⋮ ⋯ ⋮ ⋮ ⋮ a n 1 a n 2 a n 3 ⋯ a n n ∣ b n ] \left[ \begin{matrix} a_{11}&a_{12}&a_{13}&\cdots&a_{1n}&|&b_1\\ a_{21}&a_{22}&a_{23}&\cdots&a_{2n}&|&b_2\\ a_{31}&a_{32}&a_{33}&\cdots&a_{3n}&|&b_3\\ \vdots&\vdots&\vdots&\cdots&\vdots&\vdots&\vdots\\ a_{n1}&a_{n2}&a_{n3}&\cdots&a_{nn}&|&b_n\\ \end{matrix} \right] a11a21a31an1a12a22a32an2a13a23a33an3a1na2na3nannb1b2b3bn

接下来就是我们在线性代数上学到的求解线性方程组的方法。对增广矩阵做一系列的初等变换。

追赶法

方程组:
A x = D Ax = D Ax=D
其中
D = [ d 1 d 2 d 3 ⋯ d n ] T x = [ x 1 x 2 x 3 ⋯ x n ] T \begin{aligned} D &= \left[ \begin{matrix} d_1&d_2&d_3&\cdots&d_n \end{matrix} \right]^T \\ x &= \left[ \begin{matrix} x_1 & x_2 & x_3 \cdots & x_n \end{matrix} \right]^T \end{aligned} Dx=[d1d2d3dn]T=[x1x2x3xn]T

这里涉及到了一个矩阵分解,同时系数矩阵还需要时稀疏矩阵

系数矩阵如下
A = [ b 1 c 1 a 2 b 2 c 2   a 3 b 3 c 3     a 4 b 4 c 4     ⋱ ⋱ ⋱       ⋱ ⋱ c n − 1       a n b n ] A = \left[ \begin{matrix} b_1&c_1 & &\\ a_2&b_2 & c_2 \\ \ &a_3&b_3 &c_3\\ \ &\ & a_4 &b_4&c_4\\ \ &\ & &\ddots&\ddots &\ddots\\ \ &\ & &\ &\ddots&\ddots &c_{n-1}\\ \ &\ & &\ & &a_n &b_n\\ \end{matrix} \right] A=b1a2     c1b2a3    c2b3a4c3b4  c4ancn1bn
下面对这个系数矩阵进行分解,这里的分解出来矩阵如下
A = L U A = LU A=LU
其中 L L L是下三角矩阵,
L = [ α 1 γ 2 α 2 γ 3 α 3 γ 4 α 4 ⋱ ⋱ γ n α n ] L = \left[ \begin{matrix} \alpha_1\\ \gamma_2 & \alpha_2\\ &\gamma_3 & \alpha_3\\ &&\gamma_4 & \alpha_4\\ &&&\ddots & \ddots\\ &&&&\gamma_{n}& \alpha_{n}\\ \end{matrix} \right] L=α1γ2α2γ3α3γ4α4γnαn

U U U是上三角单位矩阵
U = [ 1 β 1 1 β 2 1 β 3 ⋱ ⋱ 1 β n − 1 1 ] U = \left[ \begin{matrix} 1&\beta_1\\ &1&\beta_2\\ &&1&\beta_3\\ &&&\ddots&\ddots\\ &&&&1&\beta_{n-1}\\ % 这是注释 % \% 这是转义字符 &&&&&1\\ \end{matrix} \right] U=1β11β21β31βn11
定义 y y y矩阵
y = [ y 1 y 2 y 3 ⋯ y n ] T y = \left[ \begin{matrix} y_1 & y_2 & y_3 & \cdots & y_n \end{matrix} \right]^T y=[y1y2y3yn]T
并且, y y y矩阵满足下面的公式
L U x = D { U x = y L y = D \begin{aligned} LUx = D\\ \begin{cases} Ux &= y \\ Ly &= D\\ \end{cases} \end{aligned} LUx=D{UxLy=y=D

这两个矩阵相乘结果 A A A,所以这两个矩阵中的元素 α , β , γ \alpha,\beta,\gamma α,β,γ满足下面的条件(下面的条件是怎么求出来的,我们不需要知道,只需要将下面的条件记住,就能编程序了,当然,如果是用手算的话,感觉不需要记下面的公式也行。计算物理就是用计算机来算的)
{ α 1 = b 1 α i = b i − a i β i − 1 β 1 = c 1 b 1 β i = c i b i − a i β i − 1 γ i = a i \begin{cases} \alpha_1 = b_1\\ \alpha_i = b_i - a_i \beta_{i-1}\\ \beta_1 = \frac{c_1}{b_1}\\ \beta_i = \frac{c_i}{b_i - a_i \beta_{i-1}}\\ \gamma_i = a_i \end{cases} α1=b1αi=biaiβi1β1=b1c1βi=biaiβi1ciγi=ai
按照上面的条件,通过迭代就能得到所有的 α , β \alpha,\beta α,β

下面再看 L y = D Ly = D Ly=D
[ α 1 γ 2 α 2 γ 3 α 3 γ 4 α 4 ⋱ ⋱ γ n α n ] [ y 1 y 2 y 3 y 4 ⋮ y n ] = [ d 1 d 2 d 3 d 4 ⋮ d n ] % L矩阵 \left[ \begin{matrix} \alpha_1\\ \gamma_2 & \alpha_2\\ &\gamma_3 & \alpha_3\\ &&\gamma_4 & \alpha_4\\ &&&\ddots & \ddots\\ &&&&\gamma_{n}& \alpha_{n}\\ \end{matrix} \right] % y矩阵 \left[ \begin{matrix} y_1 \\ y_2 \\ y_3 \\y_4\\ \vdots \\ y_n \end{matrix} \right]= % D矩阵 \left[ \begin{matrix} d_1\\d_2\\d_3\\d_4\\ \vdots\\d_n \end{matrix} \right] α1γ2α2γ3α3γ4α4γnαny1y2y3y4yn=d1d2d3d4dn
等效于
{ α 1 y 1 = d 1 γ 2 y 1 + α 2 y 2 = d 2 γ 3 y 2 + α 3 y 3 = d 3 ⋯ ⋯ γ n y n − 1 + α n y n = d n \begin{cases} \begin{aligned} \alpha_1 y_1 = d_1\\ \gamma_2 y_1 + \alpha_2 y_2 = d_2\\ \gamma_3 y_2 + \alpha_3 y_3 = d_3\\ \quad \cdots \cdots\quad\quad\\ \gamma_ny_{n-1} + \alpha_n y_n = d_n \end{aligned} \end{cases} α1y1=d1γ2y1+α2y2=d2γ3y2+α3y3=d3γnyn1+αnyn=dn
通过上面的第一行得到
y 1 = d 1 α 1 = d 1 b 1 y_1 = \frac{d_1}{\alpha_1} = \frac{d_1}{b_1} y1=α1d1=b1d1
第二行
y 2 = d 2 − γ 2 y 1 α 2 = d 2 − γ 2 y 1 b 2 − a 2 β 1 \begin{aligned} y_2 &= \frac{d_2 - \gamma_2 y_1}{\alpha_2} \\ &=\frac{d_2 - \gamma_2 y_1}{b_2 - a_2 \beta_1} \end{aligned} y2=α2d2γ2y1=b2a2β1d2γ2y1
从上面的方程组中可以从第二行进行递推出来,
y i = d i − γ i y i − 1 b i − a 2 β i − 1 y_i =\frac{d_i - \gamma_i y_{i-1}}{b_i - a_2 \beta_{i-1}} yi=bia2βi1diγiyi1
根据上面的递推公式,就能得到全部的 y y y,然后再利用下面的公式,求出来所有的 x x x
U x = y Ux = y Ux=y
将上面的公式展开
[ 1 β 1 1 β 2 1 β 3 ⋱ ⋱ 1 β n − 1 1 ] [ x 1 x 2 x 3 ⋮ x n − 1 x n ] = [ y 1 y 2 y 3 ⋮ y n − 1 y n ] L % U矩阵 \left[ \begin{matrix} 1&\beta_1\\ &1&\beta_2\\ &&1&\beta_3\\ &&&\ddots&\ddots\\ &&&&1&\beta_{n-1}\\ &&&&&1\\ \end{matrix} \right] % x矩阵 \left[ \begin{matrix} x_1 \\ x_2 \\ x_3 \\ \vdots \\x_{n-1} \\ x_n \end{matrix} \right]= % y矩阵 \left[ \begin{matrix} y_1 \\ y_2 \\ y_3 \\ \vdots\\y_{n-1} \\ y_n \end{matrix} \right]L 1β11β21β31βn11x1x2x3xn1xn=y1y2y3yn1ynL
所以
{ x n = y n x i = y i − β i x i + 1 \begin{cases} \begin{aligned} x_n &= y_n\\ x_i &= y_i - \beta_i x_{i+1} \end{aligned} \end{cases} {xnxi=yn=yiβixi+1
通过上面的递推公式最终能够得到全部的 x x x

迭代法求解线性方程组

对变量个数较多的线性方程组,常用迭代法进行求解。

  1. 设计一个迭代公式,任选一个初始向量 X 0 → \overrightarrow{X^0} X0

  2. 计算 X 1 → , ⋯   , X k → , ⋯ \overrightarrow{X^1},\cdots,\overrightarrow{X^k},\cdots X1 ,,Xk ,

  3. 若该向量序列收敛,其极限值为原线性方程组的解,即
    lim ⁡ k → ∞ x i ( k ) = x i ∗ \lim_{k\rightarrow\infty}x_i^{(k)} = x_i^* klimxi(k)=xi

  4. X ∗ → = ( X 1 ∗ , X 2 ∗ , ⋯   , X n ∗ ) \overrightarrow{X^*} = \big( X_1^*,X_2^*, \cdots,X_n^*\big) X =(X1,X2,,Xn)

  5. 迭代法的求解过程就相当于是求极限的过程

不明白的地方,上面的 X ∗ X^* X到底是什么意思?

雅克比迭代法

{ a 11 x 1 + a 12 x 2 + ⋯ + a 1 n x n = b 1 a 21 x 1 + a 22 x 2 + ⋯ + a 2 n x n = b 2 ⋮ a n 1 x 1 + a n 2 x 2 + ⋯ + a n n x n = b 1 (1) \begin{cases} \begin{aligned} a_{11} x_1 + a_{12}x_2 + \cdots+a_{1n}x_n = b_{1}\\ a_{21} x_1 + a_{22}x_2 + \cdots+a_{2n}x_n = b_{2}\\ \vdots\quad\quad\quad\quad\quad\quad\quad\quad \\ a_{n1} x_1 + a_{n2}x_2 + \cdots+a_{nn}x_n = b_{1}\\ \end{aligned} \end{cases}\tag{1} a11x1+a12x2++a1nxn=b1a21x1+a22x2++a2nxn=b2an1x1+an2x2++annxn=b1(1)

上面的等效于
A x = B (2) Ax = B\tag{2} Ax=B(2)
a i i ≠ 0 a_{ii}\neq 0 aii=0,则
{ x 1 ( 1 ) = 1 a 11 ( − a 12 x 2 ( 0 ) − a 13 x 3 ( 0 ) − ⋯ − a 1 n x n ( 0 ) + b 1 ) x 2 ( 1 ) = 1 a 22 ( − a 21 x 1 ( 0 ) − a 23 x 3 ( 0 ) − ⋯ − a 2 n x n ( 0 ) + b 2 ) ⋮ x n ( 1 ) = 1 a n n ( − a n 1 x 1 ( 0 ) − a n 2 x 2 ( 0 ) − ⋯ − a 1 , n − 1 x n − 1 ( 0 ) + b n ) (3) \begin{cases} x_1^{(1)} = \frac{1}{a_{11}}(-a_{12}x_2^{(0)} - a_{13}x_3^{(0)}-\cdots-a_{1n}x_n^{(0)} + b_1)\\ x_2^{(1)} = \frac{1}{a_{22}}(-a_{21}x_1^{(0)} - a_{23}x_3^{(0)}-\cdots-a_{2n}x_n^{(0)} + b_2)\\ \quad\quad\quad\quad\quad\quad\quad\quad\vdots\quad\quad\\ x_n^{(1)} = \frac{1}{a_{nn}}(-a_{n1}x_1^{(0)} - a_{n2}x_2^{(0)}-\cdots-a_{1,n-1}x_{n-1}^{(0)} + b_n)\\ \end{cases}\tag{3} x1(1)=a111(a12x2(0)a13x3(0)a1nxn(0)+b1)x2(1)=a221(a21x1(0)a23x3(0)a2nxn(0)+b2)xn(1)=ann1(an1x1(0)an2x2(0)a1,n1xn1(0)+bn)(3)
定义矩阵:
D = [ a 11 a 22 a 33 ⋱ a n n ] , L = [ 0 a 21 0 a 31 a 32 0 ⋮ ⋮ ⋮ ⋱ a n 1 a n 2 a n 3 ⋯ a n n ] , U = [ 0 a 12 a 13 ⋯ a 1 n 0 a 23 ⋯ a 2 n 0 ⋯ a 3 n ⋱ ⋮ 0 ] D = \left[ \begin{matrix} a_{11} \\ & a_{22} \\ & & a_{33} \\ & & & \ddots \\ & & & & a_{nn} \\ \end{matrix} \right] , L = \left[ \begin{matrix} 0\\ a_{21}&0\\ a_{31}&a_{32}&0\\ \vdots&\vdots&\vdots&\ddots\\ a_{n1}&a_{n2}&a_{n3}&\cdots&a_{nn} \end{matrix} \right] , U = \left[ \begin{matrix} 0 & a_{12} & a_{13} & \cdots & a_{1n}\\ & 0 & a_{23} & \cdots & a_{2n}\\ && 0 & \cdots & a_{3n}\\ &&&\ddots &\vdots\\ &&&&0 \end{matrix} \right] D=a11a22a33ann,L=0a21a31an10a32an20an3ann,U=0a120a13a230a1na2na3n0
观察上面的矩阵可以发现
D X ( 1 ) = − ( L + U ) X ( 0 ) + b (4) DX^{(1)} = -(L+U)X^{(0)} + b\tag{4} DX(1)=(L+U)X(0)+b(4)

X ( 1 ) = − D − 1 ( L + U ) X ( 0 ) + D − 1 b = − D − 1 ( L + U ) X ( 0 ) + d (5) \begin{aligned} X^{(1)} &= -D^{-1} (L+U)X^{(0)} + D^{-1}b\\ &=-D^{-1} (L+U)X^{(0)} + d \end{aligned}\tag{5} X(1)=D1(L+U)X(0)+D1b=D1(L+U)X(0)+d(5)
其中 d = D − 1 b d = D^{-1}b d=D1b

式(5)就是迭代公式。

将式(5)中的X标上上标,可以明显的看出其迭代性了
X k + 1 = − D − 1 ( L + U ) X k + d (6) X^{k+1} = -D^{-1}(L+U)X^{k} + d\tag{6} Xk+1=D1(L+U)Xk+d(6)
(这里的 D − 1 D^{-1} D1并不是 D D D矩阵的逆,表示的是矩阵D中的每个元素的倒数。)
我们可以先任意的取一个初值 X 0 X^{0} X0 ,然后将其带入到上面的迭代公式中,就能一次又一次的迭代,直到不满足迭代条件(误差在允许的范围之内),就能停止迭代了,

对于这种方法,我有一个非常不明白的地方,这样求 X X X有什么道理吗?为什么最后求出来的 X X X就是方程的解?还有怎么保证这个迭代是收敛的(循环可以停止,也就是怎么保证这个循环的时候,误差是愈来愈小的?)?

高斯——赛德尔迭代法

高斯——赛德尔迭代法,就是在雅克比迭代法的基础上进一步改进了。

在雅克比迭代法中,我们得到了最后的基本迭代公式,式(6)。

从式(6)中我们可以看出,每次求出的 X k + 1 X^{k+1} Xk+1只和 X k X^{k} Xk有关系,但是注意这里的 X k + 1 , X k X^{k+1},X^k Xk+1,Xk都是一个向量,并且我们求解这个向量的时候其实是需要将向量的的元素一个个的对应着求出来的,也就是说,我们先求出来 x 1 k + 1 x^{k+1}_1 x1k+1然后再求 x 2 k + 1 x^{k+1}_2 x2k+1,以此类推,容易发现,当我们要求 x i k + 1 x_i^{k+1} xik+1的时候,已经求出来的 x 1 , 2 , 3 , ⋯   , i − 1 k + 1 x_{1,2,3,\cdots,i-1}^{k+1} x1,2,3,,i1k+1都没有用上,也就是说浪费了,高斯——赛德尔迭代法就是发现了这个弊端并进行了改进,结果就是求解过程收敛速度加快了。


高斯——赛德尔迭代法,和雅克比迭代法的不同之处就是从式(3)开始的,将式(3)改写如下
{ x 1 ( 1 ) = 1 a 11 ( − a 12 x 2 ( 0 ) − a 13 x 3 ( 0 ) − ⋯ − a 1 n x n ( 0 ) + b 1 ) x 2 ( 1 ) = 1 a 22 ( − a 21 x 1 ( 1 ) − a 23 x 3 ( 0 ) − ⋯ − a 2 n x n ( 0 ) + b 2 ) ⋮ x n ( 1 ) = 1 a n n ( − a n 1 x 1 ( 1 ) − a n 2 x 2 ( 1 ) − ⋯ − a 1 , n − 1 x n − 1 ( 0 ) + b n ) (7) \begin{cases} x_1^{(1)} = \frac{1}{a_{11}}(-a_{12}x_2^{(0)} - a_{13}x_3^{(0)}-\cdots-a_{1n}x_n^{(0)} + b_1)\\ x_2^{(1)} = \frac{1}{a_{22}}(-a_{21}x_1^{(1)} - a_{23}x_3^{(0)}-\cdots-a_{2n}x_n^{(0)} + b_2)\\ \quad\quad\quad\quad\quad\quad\quad\quad\vdots\quad\quad\\ x_n^{(1)} = \frac{1}{a_{nn}}(-a_{n1}x_1^{(1)} - a_{n2}x_2^{(1)}-\cdots-a_{1,n-1}x_{n-1}^{(0)} + b_n)\\ \end{cases}\tag{7} x1(1)=a111(a12x2(0)a13x3(0)a1nxn(0)+b1)x2(1)=a221(a21x1(1)a23x3(0)a2nxn(0)+b2)xn(1)=ann1(an1x1(1)an2x2(1)a1,n1xn1(0)+bn)(7)
从这里的公式可以看出来,上面的第一行还是和式(3)一样的,但是从第二行开始就不一样了,因为从第二行开始就将之前计算出来的 x 1 ( 1 ) x_1^{(1)} x1(1)用上了。

所以高斯——赛德尔迭代法和雅克比迭代法相比,就这一点不同之处。也就相当于是一个算法上的改进吧。

{ x 1 ( 1 ) = 1 a 11 ( − a 12 x 2 ( 0 ) − a 13 x 3 ( 0 ) − ⋯ − a 1 n x n ( 0 ) + b 1 ) x 2 ( 1 ) = 1 a 22 ( − a 21 x 1 ( 1 ) − a 23 x 3 ( 0 ) − ⋯ − a 2 n x n ( 0 ) + b 2 ) ⋮ x n ( 1 ) = 1 a n n ( − a n 1 x 1 ( 1 ) − a n 2 x 2 ( 1 ) − ⋯ − a 1 , n − 1 x n − 1 ( 0 ) + b n ) (7) \begin{cases} x_1^{(1)} = \frac{1}{a_{11}}(-a_{12}x_2^{(0)} - a_{13}x_3^{(0)}-\cdots-a_{1n}x_n^{(0)} + b_1)\\ x_2^{(1)} = \frac{1}{a_{22}}(-a_{21}x_1^{(1)} - a_{23}x_3^{(0)}-\cdots-a_{2n}x_n^{(0)} + b_2)\\ \quad\quad\quad\quad\quad\quad\quad\quad\vdots\quad\quad\\ x_n^{(1)} = \frac{1}{a_{nn}}(-a_{n1}x_1^{(1)} - a_{n2}x_2^{(1)}-\cdots-a_{1,n-1}x_{n-1}^{(0)} + b_n)\\ \end{cases}\tag{7} x1(1)=a111(a12x2(0)a13x3(0)a1nxn(0)+b1)x2(1)=a221(a21x1(1)a23x3(0)a2nxn(0)+b2)xn(1)=ann1(an1x1(1)an2x2(1)a1,n1xn1(0)+bn)(7)
从这里的公式可以看出来,上面的第一行还是和式(3)一样的,但是从第二行开始就不一样了,因为从第二行开始就将之前计算出来的 x 1 ( 1 ) x_1^{(1)} x1(1)用上了。

所以高斯——赛德尔迭代法和雅克比迭代法相比,就这一点不同之处。也就相当于是一个算法上的改进吧。

评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值