数值分析重难点总结「二」
六、插值法
在工程应用中,常用一个简单的函数 y=y(x)y=y(x)y=y(x) 来近似代替函数 y=f(x)y=f(x)y=f(x),使得
y(xi)=f(xi),i=0,1,2,⋯ ,n
y(x_i)=f(x_i),\quad i=0,1,2,\cdots,n
y(xi)=f(xi),i=0,1,2,⋯,n
称函数 y=y(x)y=y(x)y=y(x) 为函数 y=f(x)y=f(x)y=f(x) 在点 x0,x1,⋯ ,xnx_0,x_1,\cdots,x_nx0,x1,⋯,xn 处的插值函数。
1. Lagrange 插值
引进记号 Pn+1(x)=∏i=0n(x−xi)=(x−x0)(x−x1)⋯(x−xn)P_{n+1}(x)=\prod\limits_{i=0}^n(x-x_i)=(x-x_0)(x-x_1)\cdots(x-x_n)Pn+1(x)=i=0∏n(x−xi)=(x−x0)(x−x1)⋯(x−xn),则有 Pn+1′(xj)=∏i=0,i≠jn(xj−xi)=(xj−x0)⋯(xj−xj−1)(xj−xj+1)⋯(xj−xn)P'_{n+1}(x_j)=\prod\limits_{i=0,i\ne j}^n(x_j-x_i)=(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_n)Pn+1′(xj)=i=0,i=j∏n(xj−xi)=(xj−x0)⋯(xj−xj−1)(xj−xj+1)⋯(xj−xn)
令 lj(x)=(x−x0)⋯(x−xj−1)(x−xj+1)⋯(x−xn)(xj−x0)⋯(xj−xj−1)(xj−xj+1)⋯(xj−xn)l_j(x)=\dfrac{(x-x_0)\cdots(x-x_{j-1})(x-x_{j+1})\cdots(x-x_n)}{(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_n)}lj(x)=(xj−x0)⋯(xj−xj−1)(xj−xj+1)⋯(xj−xn)(x−x0)⋯(x−xj−1)(x−xj+1)⋯(x−xn)
=∏i=0,i≠jnx−xixj−xi=Pn+1(x)(x−xj)Pn+1′(xj),j=0,1,2,⋯ ,n\qquad\quad\ =\prod\limits_{i=0,i\ne j}^n\dfrac{x-x_i}{x_j-x_i}=\dfrac{P_{n+1}(x)}{(x-x_j)P'_{n+1}(x_j)},\quad j=0,1,2,\cdots,n =i=0,i=j∏nxj−xix−xi=(x−xj)Pn+1′(xj)Pn+1(x),j=0,1,2,⋯,n
则 lj(x)l_j(x)lj(x) 是 nnn 次多项式,且 lj(xi)={1, j=i,0, j≠i,l_j(x_i) =\begin{cases} 1,\ j=i,\\0,\ j\ne i,\end{cases}lj(xi)={1, j=i,0, j=i, 称 l0(x),l1(x),⋯ ,ln(x)l_0(x),l_1(x),\cdots,l_n(x)l0(x),l1(x),⋯,ln(x) 为 Lagrange 插值基函数。令
y(x)=∑j=0nlj(x)f(xj) y(x) =\sum\limits_{j=0}^nl_j(x)f(x_j) y(x)=j=0∑nlj(x)f(xj)
则 y(xi)=f(xi)y(x_i)=f(x_i)y(xi)=f(xi),称之为 Lagrange 插值多项式,记为 Ln(x)L_n(x)Ln(x).
当 n=1n=1n=1 时,L1(x)L_1(x)L1(x) 为线性插值 L1(x)=x−x1x0−x1f(x0)+x−x0x1−x0f(x1).L_1(x)=\dfrac{x-x_1}{x_0-x_1}f(x_0)+\dfrac{x-x_0}{x_1-x_0}f(x_1).L1(x)=x0−x1x−x1f(x0)+x1−x0x−x0f(x1).
当 n=2n=2n=2 时,L2(x)L_2(x)L2(x) 为二次插值,即抛物线插值
L2(x)=(x−x1)(x−x2)(x0−x1)(x0−x2)f(x0)+(x−x0)(x−x2)(x1−x0)(x1−x2)f(x1)+(x−x0)(x−x1)(x2−x0)(x2−x1)f(x2) L_2(x)=\dfrac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}f(x_0)+ \dfrac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}f(x_1)+ \dfrac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}f(x_2) L2(x)=(x0−x1)(x0−x2)(x−x1)(x−x2)f(x0)+(x1−x0)(x1−x2)(x−x0)(x−x2)f(x1)+(x2−x0)(x2−x1)(x−x0)(x−x1)f(x2)
Lagrange 插值的插值余项为:
E(x)=f(x)−y(x)=f(n+1)(ξ)(n+1)!Pn+1(x),∀x∈[a,b] E(x) =f(x)-y(x) =\dfrac{f^{(n+1)}(\xi)}{(n+1)!}P_{n+1}(x),\quad\forall x \in[a,b] E(x)=f(x)−y(x)=(n+1)!f(n+1)(ξ)Pn+1(x),∀x∈[a,b]
其中 Pn+1(x)=∏i=0n(x−xi)P_{n+1}(x)=\prod\limits_{i=0}^n(x-x_i)Pn+1(x)=i=0∏n(x−xi),ξ∈[a,b]\xi\in[a,b]ξ∈[a,b] 且与 xxx 有关。
2. Newton 插值
一阶差商:f[x0,x1]=f(x0)−f(x1)x0−x1f[x_0,x_1]=\dfrac{f(x_0)-f(x_1)}{x_0-x_1}f[x0,x1]=x0−x1f(x0)−f(x1)
二阶差商:f[x0,x1,x2]=f[x0,x1]−f[x1,x2]x0−x2f[x_0,x_1,x_2]=\dfrac{f[x_0,x_1]-f[x_1,x_2]}{x_0-x_2}f[x0,x1,x2]=x0−x2f[x0,x1]−f[x1,x2]
kkk 阶差商: f[x0,x1,⋯ ,xk−1,xk]=f[x0,x1,…,xk−1]−f[x1,x2,…,xk]x0−xkf[x_0,x_1,\cdots,x_{k-1},x_k]=\dfrac{f [x_0,x_1,\dots,x_{k-1}]-f [x_1,x_2,\dots,x_k]}{x_0-x_k}f[x0,x1,⋯,xk−1,xk]=x0−xkf[x0,x1,…,xk−1]−f[x1,x2,…,xk]
差商有以下性质:
1) kkk 阶差商可表示为 f(x0),f(x1),⋯ ,f(xn)f(x_0),f(x_1),\cdots,f(x_n)f(x0),f(x1),⋯,f(xn) 的线性组合,即:f[x0,x1,⋯ ,xk]=∑i=0kf(xi)Pk+1′(xi).f[x_0,x_1,\cdots,x_k]=\sum\limits_{i=0}^k\dfrac{f(x_i)}{P'_{k+1}(x_i)}.f[x0,x1,⋯,xk]=i=0∑kPk+1′(xi)f(xi).
2) 差商有对称性:f[x0,x1,⋯ ,xk]=f[x1,x0,⋯ ,xk]=f[x1,x2,…,xk,x0].f[x_0,x_1,\cdots,x_k]=f[x_1,x_0,\cdots,x_k]=f[x_1,x_2,\dots,x_k,x_0].f[x0,x1,⋯,xk]=f[x1,x0,⋯,xk]=f[x1,x2,…,xk,x0].
3) f[x0,x1,⋯ ,xk]=f(k)(ξ)k!f[x_0,x_1,\cdots,x_k]=\dfrac{f^{(k)}(\xi)}{k!}f[x0,x1,⋯,xk]=k!f(k)(ξ),其中 ξ\xiξ 与节点 x0,x1,⋯ ,xkx_0,x_1,\cdots,x_kx0,x1,⋯,xk 有关。
Newton 插值多项式为:
Nn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+⋯+f[x0,x1,⋯ ,xn](x−x0)(x−x1)⋯(x−xn−1)
\begin{aligned}
N_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,\cdots,x_n](x-x_0)(x-x_1)\cdots(x-x_{n-1})
\end{aligned}
Nn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+⋯+f[x0,x1,⋯,xn](x−x0)(x−x1)⋯(x−xn−1)
计算差商时,可以作如下差商表:
| xix_ixi | f(xi)f(x_i)f(xi) | 一阶差商 | 二阶差商 | 三阶差商 |
|---|---|---|---|---|
| x0x_0x0 | f(x0)f(x_0)f(x0) | |||
| x1x_1x1 | f(x1)f(x_1)f(x1) | f[x0,x1]f[x_0,x_1]f[x0,x1] | ||
| x2x_2x2 | f(x2)f(x_2)f(x2) | f[x1,x2]f[x_1,x_2]f[x1,x2] | f[x0,x1,x2]f[x_0,x_1,x_2]f[x0,x1,x2] | |
| x3x_3x3 | f(x3)f(x_3)f(x3) | f[x2,x3]f[x_2,x_3]f[x2,x3] | f[x1,x2,x3]f[x_1,x_2,x_3]f[x1,x2,x3] | f[x0,x1,x2,x3]f[x_0,x_1,x_2,x_3]f[x0,x1,x2,x3] |
| ⋮\vdots⋮ | ⋮\vdots⋮ | ⋮\vdots⋮ | ⋮\vdots⋮ | ⋮\vdots⋮ |
其余项
E(x)=f(x)−Nn(x)=f[x,x0,x1,⋯ ,xn]Pn+1(x) E(x) =f(x)-N_n(x)=f[x,x_0,x_1,\cdots,x_n]P_{n+1}(x) E(x)=f(x)−Nn(x)=f[x,x0,x1,⋯,xn]Pn+1(x)
根据插值多项式的唯一性,虽然 Lagrange 插值多项式与 Newton 插值多项式的构造方式不同,但必有 Nn(x)≡Ln(x)N_n(x)\equiv L_n(x)Nn(x)≡Ln(x),因此
f[x0,x1,⋯ ,xn]=∑i=0nf(xi)Pn+1′(xi),f[x,x0,x1,⋯ ,xn]=f(n+1)(ξ)(n+1)! f[x_0,x_1,\cdots,x_n]=\sum\limits_{i=0}^n\dfrac{f(x_i)}{P'_{n+1}(x_i)},\quad f[x,x_0,x_1,\cdots,x_n]=\dfrac{f^{(n+1)}(\xi)}{(n+1)!} f[x0,x1,⋯,xn]=i=0∑nPn+1′(xi)f(xi),f[x,x0,x1,⋯,xn]=(n+1)!f(n+1)(ξ)
3. Hermite 插值
给定的函数关系中含有导数的插值即称为 Hermite 插值。
待定系数法:
求满足 P(x0)=f(x0),P(x1)=f(x1)P(x_0)=f(x_0),P(x_1)=f(x_1)P(x0)=f(x0),P(x1)=f(x1) 及 $ P’(x_0)=f’(x_0)$ 的次数不超过2次的插值多项式 P(x)P(x)P(x).
也即求 Hermite 两点两次插值,令 H2(x)=h0(x)f(x0)+h1(x)f(x1)+h‾0(x)f′(x0)H_2(x)=h_0(x)f(x_0)+h_1(x)f(x_1)+\overline{h}_0(x)f'(x_0)H2(x)=h0(x)f(x0)+h1(x)f(x1)+h0(x)f′(x0).
可知
{H2(x0)=f(x0)⇒h0(x0)=1, h1(x0)=0, h‾0(x0)=0H2(x1)=f(x1)⇒h0(x1)=0, h1(x1)=1, h‾0(x1)=0H2′(x0)=f′(x0)⇒h0′(x0)=0, h1′(x0)=0, h‾0′(x0)=1
\begin{cases}
H_2(x_0)=f(x_0)&\Rarr h_0(x_0)=1,\ h_1(x_0)=0,\ \overline{h}_0(x_0)=0\\
H_2(x_1)=f(x_1)&\Rarr h_0(x_1)=0,\ h_1(x_1)=1,\ \overline{h}_0(x_1)=0\\
H'_2(x_0)=f'(x_0)&\Rarr h'_0(x_0)=0,\ h'_1(x_0)=0,\ \overline{h}'_0(x_0)=1\\
\end{cases}
⎩⎨⎧H2(x0)=f(x0)H2(x1)=f(x1)H2′(x0)=f′(x0)⇒h0(x0)=1, h1(x0)=0, h0(x0)=0⇒h0(x1)=0, h1(x1)=1, h0(x1)=0⇒h0′(x0)=0, h1′(x0)=0, h0′(x0)=1
求解,得
h0(x)=1−(x−x0)2(x1−x0)2, h1(x)=(x−x0)2(x1−x0)2, h‾0(x)=(x−x0)(x−x1)(x0−x1)h_0(x)=1-\dfrac{(x-x_0)^2}{(x_1-x_0)^2},\ h_1(x)=\dfrac{(x-x_0)^2}{(x_1-x_0)^2},\ \overline{h}_0(x)=\dfrac{(x-x_0)(x-x_1)}{(x_0-x_1)}h0(x)=1−(x1−x0)2(x−x0)2, h1(x)=(x1−x0)2(x−x0)2, h0(x)=(x0−x1)(x−x0)(x−x1)
即 P(x)=H2(x)=[1−(x−x0)2(x1−x0)2]f(x0)+(x−x0)2(x1−x0)2f(x1)+(x−x0)(x−x1)(x0−x1)f′(x0)P(x)=H_2(x)=[1-\dfrac{(x-x_0)^2}{(x_1-x_0)^2}]f(x_0)+\dfrac{(x-x_0)^2}{(x_1-x_0)^2}f(x_1)+\dfrac{(x-x_0)(x-x_1)}{(x_0-x_1)}f'(x_0)P(x)=H2(x)=[1−(x1−x0)2(x−x0)2]f(x0)+(x1−x0)2(x−x0)2f(x1)+(x0−x1)(x−x0)(x−x1)f′(x0)
两点三次插值:
H3(x)=(1+2x−x0x1−x0)(x−x1x0−x1)2f(x0)+(1+2x−x1x0−x1)(x−x0x1−x0)2f(x1)+(x−x0)(x−x1x0−x1)2f′(x0)+(x−x1)(x−x0x1−x0)2f′(x1)
\begin{aligned}
H_3(x)=&(1+2\dfrac{x-x_0}{x_1-x_0})(\dfrac{x-x_1}{x_0-x_1})^2f(x_0)+
(1+2\dfrac{x-x_1}{x_0-x_1})(\dfrac{x-x_0}{x_1-x_0})^2f(x_1)+\\[8pt]
&(x-x_0)(\dfrac{x-x_1}{x_0-x_1})^2f'(x_0)+(x-x_1)(\dfrac{x-x_0}{x_1-x_0})^2f'(x_1)
\end{aligned}
H3(x)=(1+2x1−x0x−x0)(x0−x1x−x1)2f(x0)+(1+2x0−x1x−x1)(x1−x0x−x0)2f(x1)+(x−x0)(x0−x1x−x1)2f′(x0)+(x−x1)(x1−x0x−x0)2f′(x1)
*带重节点的差商表:
首先明确一个定义,之前提到差商的性质时,有一条性质为:f[x0,x1,…,xn]=f(n)(ξ)n!f[x_0,x_1,\dots,x_n]=\dfrac{f^{(n)}(\xi)}{n!}f[x0,x1,…,xn]=n!f(n)(ξ).
因此,利用这条性质可以得到:f[xi,xi,…,xi]=f(n)(xi)n!f[x_i,x_i,\dots,x_i]=\dfrac{f^{(n)}(x_i)}{n!}f[xi,xi,…,xi]=n!f(n)(xi).
例如,给定了 f′(x1)=y1′f^{'} ( x_1 ) =y_1^{'}f′(x1)=y1′ ,则作出的重差商表为:
| xix_ixi | f(xi)f(x_i)f(xi) | 一阶差商 | 二阶差商 | 三阶差商 |
|---|---|---|---|---|
| x0x_0x0 | f(x0)f(x_0)f(x0) | |||
| x1x_1x1 | f(x1)f(x_1)f(x1) | f[x0,x1]f[x_0,x_1]f[x0,x1] | ||
| x1x_1x1 | f(x1)f(x_1)f(x1) | f[x1,x1]f[x_1,x_1]f[x1,x1] | f[x0,x1,x1]f[x_0,x_1,x_1]f[x0,x1,x1] | |
| x2x_2x2 | f(x2)f(x_2)f(x2) | f[x1,x2]f[x_1,x_2]f[x1,x2] | f[x1,x1,x2]f[x_1,x_1,x_2]f[x1,x1,x2] | f[x0,x1,x1,x2]f[x_0,x_1,x_1,x_2]f[x0,x1,x1,x2] |
| ⋮\vdots⋮ | ⋮\vdots⋮ | ⋮\vdots⋮ | ⋮\vdots⋮ | ⋮\vdots⋮ |
利用此法要比书上所说的待定系数求得插值多项式的方法简单很多。
七、函数逼近与曲线拟合
1. 函数的内积与范数
定义 7.2 设 f(x),g(x)∈C[a,b]f(x),g(x)\in C[a,b]f(x),g(x)∈C[a,b],称
(f,g)=∫abρ(x)f(x)g(x)dx
(f,g)=\int_a^b\rho(x)f(x)g(x){\rm d}x
(f,g)=∫abρ(x)f(x)g(x)dx
为函数 f(x),g(x)f(x),g(x)f(x),g(x) 在 [a,b][a,b][a,b] 上的内积,其中 ρ(x)\rho(x)ρ(x) 为有权函数,计算时一般不带权,其取值恒为 1.
定义 7.3 设 f(x)f(x)f(x) 在 C[a,b]C[a,b]C[a,b] 上的范数定义为
∣∣f∣∣1=∫ab∣f(x)∣dx∣∣f∣∣2=[∫abf2(x)dx]12∣∣f∣∣∞=maxa∈[a,b]∣f(x)∣
\begin{aligned}
||f||_1=&\int_a^b|f(x)|{\rm d}x\\
||f||_2=&\left[\int_a^bf^2(x){\rm d}x\right]^{\frac{1}{2}}\\
||f||_\infin=&\max_{a\in[a,b]}|f(x)|\\
\end{aligned}
∣∣f∣∣1=∣∣f∣∣2=∣∣f∣∣∞=∫ab∣f(x)∣dx[∫abf2(x)dx]21a∈[a,b]max∣f(x)∣
2. 最佳平方逼近
函数的最佳平方逼近,列出法方程,也称为正则方程:
[(φ0,φ0)(φ0,φ1)⋯(φ0,φm)(φ1,φ0)(φ1,φ1)⋯(φ1,φm)⋮⋮⋱⋮(φm,φ0)(φm,φ1)⋯(φm,φm)][a0a1⋮am]=[(f,φ0)(f,φ1)⋮(f,φm)]
\begin{bmatrix} (\varphi_0,\varphi_0) & (\varphi_0,\varphi_1) & \cdots & (\varphi_0,\varphi_m)\\
(\varphi_1,\varphi_0) & (\varphi_1,\varphi_1) & \cdots & (\varphi_1,\varphi_m)\\
\vdots & \vdots & \ddots & \vdots\\ (\varphi_m,\varphi_0) & (\varphi_m,\varphi_1) & \cdots & (\varphi_m,\varphi_m)\\ \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\\vdots\\ a_m \end{bmatrix} =
\begin{bmatrix} (f,\varphi_0)\\ (f,\varphi_1)\\ \vdots\\ (f,\varphi_m) \end{bmatrix}
(φ0,φ0)(φ1,φ0)⋮(φm,φ0)(φ0,φ1)(φ1,φ1)⋮(φm,φ1)⋯⋯⋱⋯(φ0,φm)(φ1,φm)⋮(φm,φm)a0a1⋮am=(f,φ0)(f,φ1)⋮(f,φm)
方程组存在一组唯一解 a0∗,a1∗,⋯ ,am∗a_0^*,a_1^*,\cdots,a_m^*a0∗,a1∗,⋯,am∗,平方误差为
∣∣E(x)∣∣22=∣∣f∣∣22−∑j=0m(f,φj)aj∗ ||E(x)||_2^2=||f||_2^2-\sum\limits_{j=0}^m(f,\varphi_j)a_j^* ∣∣E(x)∣∣22=∣∣f∣∣22−j=0∑m(f,φj)aj∗
若取基函数空间 Φm=span{1,x,x2,⋯ ,xn}, ρ(x)≡1, f(x)∈C[0,1]\Phi_m=\text{span}\{1,x,x^2,\cdots,x^n\},\ \rho(x)\equiv 1,\ f(x)\in C[0,1]Φm=span{1,x,x2,⋯,xn}, ρ(x)≡1, f(x)∈C[0,1] ,则有:
[112⋯1n1213⋯1n+1⋮⋮⋱⋮1n1n+1⋯12n−1][a0a1⋮an]=[(f,φ0)(f,φ1)⋮(f,φn)] \begin{bmatrix} 1 & \dfrac12 & \cdots & \dfrac1n\\[8pt] \dfrac12 & \dfrac13 & \cdots & \dfrac1{n+1}\\[8pt] \vdots & \vdots & \ddots & \vdots\\[8pt] \dfrac1n & \dfrac1{n+1} & \cdots & \dfrac1{2n-1}\\ \end{bmatrix} \begin{bmatrix} a_0\\[8pt] a_1\\[8pt] \vdots\\[8pt] a_n \end{bmatrix}= \begin{bmatrix} (f,\varphi_0)\\[8pt] (f,\varphi_1)\\[8pt] \vdots\\[8pt] (f,\varphi_n) \end{bmatrix} 121⋮n12131⋮n+11⋯⋯⋱⋯n1n+11⋮2n−11a0a1⋮an=(f,φ0)(f,φ1)⋮(f,φn)
3. 数据拟合与最小二乘法
列出法方程,也称为正则方程的矩阵形式为:
[(φ0,φ0)(φ0,φ1)⋯(φ0,φm)(φ1,φ0)(φ1,φ1)⋯(φ1,φm)⋮⋮⋱⋮(φm,φ0)(φm,φ1)⋯(φm,φm)][a0a1⋮am]=[(f,φ0)(f,φ1)⋮(f,φm)] \begin{bmatrix} (\varphi_0,\varphi_0) & (\varphi_0,\varphi_1) & \cdots & (\varphi_0,\varphi_m)\\ (\varphi_1,\varphi_0) & (\varphi_1,\varphi_1) & \cdots & (\varphi_1,\varphi_m)\\ \vdots & \vdots & \ddots & \vdots\\ (\varphi_m,\varphi_0) & (\varphi_m,\varphi_1) & \cdots & (\varphi_m,\varphi_m)\\ \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\\vdots\\ a_m \end{bmatrix} = \begin{bmatrix} (f,\varphi_0)\\ (f,\varphi_1)\\ \vdots\\ (f,\varphi_m) \end{bmatrix} (φ0,φ0)(φ1,φ0)⋮(φm,φ0)(φ0,φ1)(φ1,φ1)⋮(φm,φ1)⋯⋯⋱⋯(φ0,φm)(φ1,φm)⋮(φm,φm)a0a1⋮am=(f,φ0)(f,φ1)⋮(f,φm)
方程组存在一组唯一解 a0∗,a1∗,⋯ ,am∗a_0^*,a_1^*,\cdots,a_m^*a0∗,a1∗,⋯,am∗,数据拟合的平方误差为
∣∣E(x)∣∣22=∣∣f∣∣22−∑j=0maj∗(f,φj) ||E(x)||_2^2=||f||_2^2-\sum\limits_{j=0}^ma_j^*(f,\varphi_j) ∣∣E(x)∣∣22=∣∣f∣∣22−j=0∑maj∗(f,φj)
取 Φm=span{1,x,x2,⋯ ,xm}, ωi(x)≡1,m<n\Phi_m=\text{span}\{1,x,x^2,\cdots,x^m\},\ \omega_i(x)\equiv 1,m<nΦm=span{1,x,x2,⋯,xm}, ωi(x)≡1,m<n,法方程为:
[n+1∑i=0nxi⋯∑i=0nxim∑i=0Nxi∑i=0nxi2⋯∑i=0nxim+1⋮⋮⋱⋮∑i=0nxim∑i=0nxim+1⋯∑i=0nxi2m][a0a1⋮an]=[∑i=0nfi∑i=0nfixi⋮∑i=0nfixim] \begin{bmatrix} n+1 & \sum\limits_{i=0}^nx_i & \cdots & \sum\limits_{i=0}^nx_i^m\\[8pt] \sum\limits_{i=0}^Nx_i & \sum\limits_{i=0}^nx_i^2 & \cdots & \sum\limits_{i=0}^nx_i^{m+1}\\[8pt] \vdots & \vdots & \ddots & \vdots\\[8pt] \sum\limits_{i=0}^nx_i^m & \sum\limits_{i=0}^nx_i^{m+1} & \cdots & \sum\limits_{i=0}^nx_i^{2m}\\[8pt] \end{bmatrix} \begin{bmatrix} a_0\\[9pt] a_1\\[9pt] \vdots\\[9pt] a_n \end{bmatrix} = \begin{bmatrix} \sum\limits_{i=0}^nf_i\\[8pt] \sum\limits_{i=0}^nf_ix_i\\[8pt] \vdots\\[8pt] \sum\limits_{i=0}^nf_ix_i^m \end{bmatrix} n+1i=0∑Nxi⋮i=0∑nximi=0∑nxii=0∑nxi2⋮i=0∑nxim+1⋯⋯⋱⋯i=0∑nximi=0∑nxim+1⋮i=0∑nxi2ma0a1⋮an=i=0∑nfii=0∑nfixi⋮i=0∑nfixim
八、数值积分与数值微分
1. 代数精度
定义 8.1 对于 ∫abf(x)dx≈∑k=0nAkf(xk)\int_a^bf(x){\rm d}x\approx\sum\limits_{k=0}^nA_kf(x_k)∫abf(x)dx≈k=0∑nAkf(xk) 的数值求积公式近似表示求积结果,xkx_kxk 为求积节点,AkA_kAk 为求积系数。
若该求积公式对于任意 f(x)=Pr(x)f(x)=P_r(x)f(x)=Pr(x) 精确成立,即 ∫abPr(x)dx=∑k=0nAkPr(xk)\int_a^bP_r(x)dx=\sum\limits_{k=0}^nA_kP_r(x_k)∫abPr(x)dx=k=0∑nAkPr(xk) ,但是至少有一个 f(x)=Pr+1(x)f(x)=P_{r+1}(x)f(x)=Pr+1(x) 不精确成立,即 ∫abPr+1(x)dx≠∑k=0nAkPr+1(xk)\int_a^bP_{r+1}(x)dx\ne\sum\limits_{k=0}^nA_kP_{r+1}(x_k)∫abPr+1(x)dx=k=0∑nAkPr+1(xk) ,则称该求积公式具有 rrr 次代数精度,其中 Pr(x)P_r(x)Pr(x) 表示 rrr 次多项式。
定理 8.1 对任意给定的 n+1n+1n+1 个相异节点 a⩽x0<x1<⋯<xn⩽ba\leqslant x_0<x_1<\cdots<x_n\leqslant ba⩽x0<x1<⋯<xn⩽b,总存在相应的求积系数 A0,A1,⋯ ,AnA_0,A_1,\cdots,A_nA0,A1,⋯,An 使求积公式至少具有 n 次代数精度。
2. 求解求积公式的节点或系数
列线性方程组
{A0+A1+⋯+An=∫ab1dx=b−aA0x0+A1x1+⋯+Anxn=∫abxdx=12(b2−a2)A0x02+A1x12+⋯+Anxn2=∫abx2dx=13(b3−a3)⋯⋯⋯A0x0n+A1x1n+⋯+Anxnn=∫abxndx=1n+1(bn+1−an+1) \begin{cases} A_0&+A_1&+\cdots&+A_n&=\int_a^b1{\rm d}x&=b-a\\[2ex] A_0x_0&+A_1x_1&+\cdots&+A_nx_n&=\int_a^bx{\rm d}x&=\dfrac{1}{2}(b^2-a^2)\\[2ex] A_0x_0^2&+A_1x_1^2&+\cdots&+A_nx_n^2&=\int_a^bx^2{\rm d}x&=\dfrac{1}{3}(b^3-a^3)\\[2ex] & &\cdots&\cdots&\cdots\\ A_0x_0^n&+A_1x_1^n&+\cdots&+A_nx_n^n&=\int_a^bx^n{\rm d}x&=\dfrac{1}{n+1}(b^{n+1}-a^{n+1}) \end{cases} ⎩⎨⎧A0A0x0A0x02A0x0n+A1+A1x1+A1x12+A1x1n+⋯+⋯+⋯⋯+⋯+An+Anxn+Anxn2⋯+Anxnn=∫ab1dx=∫abxdx=∫abx2dx⋯=∫abxndx=b−a=21(b2−a2)=31(b3−a3)=n+11(bn+1−an+1)
3. Newton-Cotes求积公式
将区间 [a,b] n[a,b]\ n[a,b] n 等分,步长 h=b−anh=\dfrac{b-a}{n}h=nb−a.
易知,n=1n=1n=1 时,求积公式为梯形公式:
∫abf(x)dx≈b−a2[f(a)+f(b)] \int_a^bf(x){\rm d}x\approx\frac{b-a}{2}\left[f(a)+f(b)\right] ∫abf(x)dx≈2b−a[f(a)+f(b)]
余项为:
E[f]=−h312f′′(η),η∈(a,b) E[f]=-\frac{h^3}{12}f''(\eta),\quad\eta\in(a,b) E[f]=−12h3f′′(η),η∈(a,b)
n=2n=2n=2 时,求积公式为 Simpson 公式
∫abf(x)dx≈b−a6[f(a)+4f(a+b2)+f(b)] \int_a^bf(x){\rm d}x\approx \frac{b-a}6\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b) \right] ∫abf(x)dx≈6b−a[f(a)+4f(2a+b)+f(b)]
余项为:
E[f]=−h590f(4)(η),η∈(a,b) E[f]=-\frac{h^5}{90}f^{(4)}(\eta),\quad\eta\in(a,b) E[f]=−90h5f(4)(η),η∈(a,b)
n=4n=4n=4 时,求积公式为 Cotes 公式
∫abf(x)dx=b−a90[7f(x0)+32f(x1)+12f(x2)+32f(x3)+7f(x4)]−8h7945f(6)(η) \int_a^bf(x){\rm d}x=\frac{b-a}{90}\left[7f(x_0)+32f(x_1)+12f(x_2)+32f(x_3)+7f(x_4)\right]- \frac{8h^7}{945}f^{(6)}(\eta) ∫abf(x)dx=90b−a[7f(x0)+32f(x1)+12f(x2)+32f(x3)+7f(x4)]−9458h7f(6)(η)
其中,xi=a+ib−a4(i=0,1,2,3,4).x_i=a+i\dfrac{b-a}{4}(i=0,1,2,3,4).xi=a+i4b−a(i=0,1,2,3,4).
4. 复化求积公式
将积分区间 [a,b][a,b][a,b] 划分为 nnn 等分,步长 h=b−anh=\dfrac{b-a}nh=nb−a ,分点为 xi=a+ih,i=0,1,2,⋯ ,nx_i=a+ih,\quad i=0,1,2,\cdots,nxi=a+ih,i=0,1,2,⋯,n . 所谓复化求积法,是指先用低阶的 Newton-Cotes 公式求得每个子区间 [xi,xi+1][x_i,x_{i+1}][xi,xi+1] 上的积分值 IiI_iIi ,然后再求和,用 ∑i=0n−1Ii\sum\limits_{i=0}^{n-1}I_ii=0∑n−1Ii 作为所求积分 III 的近似值.则有:
复化梯形公式:
Tn=∑i=0n−1h2[f(xi)+f(xi+1)]=h2[f(a)+2∑i=1n−1f(xi)+f(b)] T_n=\sum\limits_{i=0}^{n-1}\frac h2\left[f(x_i)+f(x_{i+1})\right]=\frac h2\left[f(a)+2\sum\limits_{i=1}^{n-1}f(x_i)+f(b)\right] Tn=i=0∑n−12h[f(xi)+f(xi+1)]=2h[f(a)+2i=1∑n−1f(xi)+f(b)]
余项为:
E[f]=I−Tn=−b−a12h2f′′(η),η∈(a,b) E[f]=I-T_n=-\frac{b-a}{12}h^2f''(\eta),\quad\eta\in(a,b) E[f]=I−Tn=−12b−ah2f′′(η),η∈(a,b)
复化Simpson公式:
Sn=∑i=0n−1h6[f(xi)+4f(xi+12)+f(xi+1)]=h6[f(a)+4∑i=0n−1f(xi+12)+2∑i=1nf(xi+1)+f(b)] \begin{aligned} S_n&=\sum_{i=0}^{n-1}\frac h6\left[f(x_i)+4f\left(x_{i+\frac12}\right)+f(x_{i+1})\right]\\&=\frac h6\left[f(a)+4\sum_{i=0}^{n-1}f\left(x_{i+\frac12}\right)+2\sum_{i=1}^{n}f\left(x_{i+1}\right)+f(b)\right] \end{aligned} Sn=i=0∑n−16h[f(xi)+4f(xi+21)+f(xi+1)]=6h[f(a)+4i=0∑n−1f(xi+21)+2i=1∑nf(xi+1)+f(b)]
余项为:
E[f]=I−Sn=−b−a2880h4f(4)(η),η∈(a,b) E[f]=I-S_n=-\frac{b-a}{2880}h^4f^{(4)}(\eta),\quad\eta\in(a,b) E[f]=I−Sn=−2880b−ah4f(4)(η),η∈(a,b)
5. Romberg 求积公式
实际计算中,一般并不需要使用 Simpson 公式与 Cotes 公式,而是可以由递推技巧,仅需计算梯形公式就能一步一步推导得到 Romberg 公式结果,其推导公式为
{Sn=43T2n−13TnCn=1615S2n−115SnRn=6463C2n−163Cn⋯ \begin{cases} S_n=\dfrac{4}{3}T_{2n}-\dfrac{1}{3}T_{n}\\[2ex] C_n=\dfrac{16}{15}S_{2n}-\dfrac{1}{15}S_{n}\\[2ex] R_n=\dfrac{64}{63}C_{2n}-\dfrac{1}{63}C_{n}\\ \cdots \end{cases} ⎩⎨⎧Sn=34T2n−31TnCn=1516S2n−151SnRn=6364C2n−631Cn⋯
其中,TnT_nTn 为复化梯形公式的求积结果,n=1,2,⋯ ,nn=1,2,\cdots,nn=1,2,⋯,n. 故要计算 R1R_1R1,仅需计算 T1,T2,T4,T8T_1,T_2,T_4,T_8T1,T2,T4,T8 即可。一般可以列出如下的表格进行计算:
| TnT_nTn | SnS_nSn | CnC_nCn | RnR_nRn |
|---|---|---|---|
| T1T_1T1 | |||
| T2T_2T2 | S1S_1S1 | ||
| T4T_4T4 | S2S_2S2 | C1C_1C1 | |
| T8T_8T8 | S4S_4S4 | C2C_2C2 | R1R_1R1 |
6. Gauss 型求积公式
定义 8.3 若求积公式
∫abf(x)dx=∑j=0nAjf(xj)+E(f) \int_a^bf(x){\rm d}x=\sum\limits_{j=0}^nA_jf(x_j)+E(f) ∫abf(x)dx=j=0∑nAjf(xj)+E(f)
具有 2n+12n+12n+1 次的代数精度,则称此求积公式为 Gauss 型求积公式,相应的求积节点 xjx_jxj 称为 Gauss 点。
Gauss_Legendre 求积公式:
| 节点数 | xjx_jxj | AjA_jAj |
|---|---|---|
| 111 | 000 | 222 |
| 222 | ±0.5773502692\pm0.5773502692±0.5773502692 | 111 |
| 333 | ±0.77459666920\pm0.7745966692\\0±0.77459666920 | 0.55555555560.88888888890.5555555556\\0.88888888890.55555555560.8888888889 |
7. 数值微分
插值型数值微分:
两点公式 (n=1)(n=1)(n=1) :
{f′(x0)=1h[f(x1)−f(x0)]−h2f′′(ξ1)f′(x1)=1h[f(x1)−f(x0)]+h2f′′(ξ2) \begin{cases} f'(x_0)=\dfrac{1}{h}[f(x_1)-f(x_0)]-\dfrac{h}{2}f''(\xi_1)\\[2ex] f'(x_1)=\dfrac{1}{h}[f(x_1)-f(x_0)]+\dfrac{h}{2}f''(\xi_2) \end{cases} ⎩⎨⎧f′(x0)=h1[f(x1)−f(x0)]−2hf′′(ξ1)f′(x1)=h1[f(x1)−f(x0)]+2hf′′(ξ2)
三点公式 (n=2)(n=2)(n=2) :
{f′(x0)=12h[−3f(x0)+4f(x1)−f(x2)]+h23f′′′(ξ1)f′(x1)=12h[f(x2)−f(x0)]−h26f′′′(ξ2)(中心公式)f′(x2)=12h[f(x0)−4f(x1)+3f(x2)]+h23f′′′(ξ3) \begin{cases} f'(x_0)=\dfrac{1}{2h}[-3f(x_0)+4f(x_1)-f(x_2)]+\dfrac{h^2}{3}f'''(\xi_1)\\[2ex] f'(x_1)=\dfrac{1}{2h}[f(x_2)-f(x_0)]-\dfrac{h^2}{6}f'''(\xi_2)\quad\text{(中心公式)}\\[2ex] f'(x_2)=\dfrac{1}{2h}[f(x_0)-4f(x_1)+3f(x_2)]+\dfrac{h^2}{3}f'''(\xi_3) \end{cases} ⎩⎨⎧f′(x0)=2h1[−3f(x0)+4f(x1)−f(x2)]+3h2f′′′(ξ1)f′(x1)=2h1[f(x2)−f(x0)]−6h2f′′′(ξ2)(中心公式)f′(x2)=2h1[f(x0)−4f(x1)+3f(x2)]+3h2f′′′(ξ3)
用差商代替微商:
1) 向前差商:f′(x)=f(x+h)−f(x)h−h2f′′(ξ)f'(x)=\dfrac{f(x+h)-f(x)}{h}-\dfrac{h}{2}f''(\xi)f′(x)=hf(x+h)−f(x)−2hf′′(ξ)
2) 向后差商:f′(x)=f(x)−f(x−h)h+h2f′′(ξ)f'(x)=\dfrac{f(x)-f(x-h)}{h}+\dfrac{h}{2}f''(\xi)f′(x)=hf(x)−f(x−h)+2hf′′(ξ)
3) 中心差商:f′(x)=f(x+h)−f(x−h)2h−h26f′′′(ξ)f'(x)=\dfrac{f(x+h)-f(x-h)}{2h}-\dfrac{h^2}{6}f'''(\xi)f′(x)=2hf(x+h)−f(x−h)−6h2f′′′(ξ)
4) 二阶中心差商:f′′(x)=f(x+h)−2f(x)+f(x−h)h2−h212f(4)(ξ)f''(x)=\dfrac{f(x+h)-2f(x)+f(x-h)}{h^2}-\dfrac{h^2}{12}f^{(4)}(\xi)f′′(x)=h2f(x+h)−2f(x)+f(x−h)−12h2f(4)(ξ)
九、常微分方程的数值解法
主要考虑微分方程的初值问题
{dydx=f(x,y),a⩽x⩽by(a)=y0 \begin{cases} \dfrac{{\rm d}y}{{\rm d}x}=f(x,y),\quad a\leqslant x\leqslant b\\[2ex] y(a)=y_0 \end{cases} ⎩⎨⎧dxdy=f(x,y),a⩽x⩽by(a)=y0
1. Euler 方法
y(xi+1)≈y(xi)+hf(xi,y(xi)) y(x_{i+1})\approx y(x_i)+hf(x_i,y(x_i)) y(xi+1)≈y(xi)+hf(xi,y(xi))
2. Euler 预测-校正法
{y‾i+1=yi+hf(xi,yi)yi+1=yi+h2[f(xi,yi)+f(xi+1,y‾i+1)] \begin{cases} \overline{y}_{i+1}=y_i+hf(x_i,y_i)\\[2ex] y_{i+1}=y_i+\dfrac{h}{2}[f(x_i,y_i)+f(x_{i+1},\overline{y}_{i+1})] \end{cases} ⎩⎨⎧yi+1=yi+hf(xi,yi)yi+1=yi+2h[f(xi,yi)+f(xi+1,yi+1)]
3. 四阶经典 Runge-Kutta 公式
yi+1=yi+16(K1+2K2+2K3+K4){K1=hf(xi,yi)K2=hf(xi+h2,yi+K12)K3=hf(xi+h2,yi+K22)K4=hf(xi+h,yi+K3) y_{i+1}=y_i+\dfrac{1}{6}(K_1+2K_2+2K_3+K_4)\\[2ex] \begin{cases} K_1=hf(x_i,y_i)\\[2ex] K_2=hf(x_i+\dfrac{h}{2},y_i+\dfrac{K_1}{2})\\[2ex] K_3=hf(x_i+\dfrac{h}{2},y_i+\dfrac{K_2}{2})\\[2ex] K_4=hf(x_i+h,y_i+K_3)\\[2ex] \end{cases} yi+1=yi+61(K1+2K2+2K3+K4)⎩⎨⎧K1=hf(xi,yi)K2=hf(xi+2h,yi+2K1)K3=hf(xi+2h,yi+2K2)K4=hf(xi+h,yi+K3)
10万+

被折叠的 条评论
为什么被折叠?



