数值分析重难点总结「二」

数值分析重难点总结「二」

六、插值法

在工程应用中,常用一个简单的函数 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=0n(xxi)=(xx0)(xx1)(xxn),则有 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=jn(xjxi)=(xjx0)(xjxj1)(xjxj+1)(xjxn)

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)=(xjx0)(xjxj1)(xjxj+1)(xjxn)(xx0)(xxj1)(xxj+1)(xxn)

 =∏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=jnxjxixxi=(xxj)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=0nlj(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)=x0x1xx1f(x0)+x1x0xx0f(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)=(x0x1)(x0x2)(xx1)(xx2)f(x0)+(x1x0)(x1x2)(xx0)(xx2)f(x1)+(x2x0)(x2x1)(xx0)(xx1)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=0n(xxi)ξ∈[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]=x0x1f(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]=x0x2f[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,,xk1,xk]=x0xkf[x0,x1,,xk1]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=0kPk+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](xx0)+f[x0,x1,x2](xx0)(xx1)++f[x0,x1,,xn](xx0)(xx1)(xxn1)

计算差商时,可以作如下差商表:

xix_ixif(xi)f(x_i)f(xi)一阶差商二阶差商三阶差商
x0x_0x0f(x0)f(x_0)f(x0)
x1x_1x1f(x1)f(x_1)f(x1)f[x0,x1]f[x_0,x_1]f[x0,x1]
x2x_2x2f(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_3x3f(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=0nPn+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)=0h0(x1)=0, h1(x1)=1, h0(x1)=0h0(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(x1x0)2(xx0)2, h1(x)=(x1x0)2(xx0)2, h0(x)=(x0x1)(xx0)(xx1)
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(x1x0)2(xx0)2]f(x0)+(x1x0)2(xx0)2f(x1)+(x0x1)(xx0)(xx1)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+2x1x0xx0)(x0x1xx1)2f(x0)+(1+2x0x1xx1)(x1x0xx0)2f(x1)+(xx0)(x0x1xx1)2f(x0)+(xx1)(x1x0xx0)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_ixif(xi)f(x_i)f(xi)一阶差商二阶差商三阶差商
x0x_0x0f(x0)f(x_0)f(x0)
x1x_1x1f(x1)f(x_1)f(x1)f[x0,x1]f[x_0,x_1]f[x0,x1]
x1x_1x1f(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_2x2f(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.2f(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.3f(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∣∣∞=max⁡a∈[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} ∣∣f1=∣∣f2=∣∣f=abf(x)dx[abf2(x)dx]21a[a,b]maxf(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)a0a1am=(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=∣∣f22j=0m(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} 121n12131n+11n1n+112n11a0a1an=(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)a0a1am=(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=∣∣f22j=0maj(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=0Nxii=0nximi=0nxii=0nxi2i=0nxim+1i=0nximi=0nxim+1i=0nxi2ma0a1an=i=0nfii=0nfixii=0nfixim

八、数值积分与数值微分

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)dxk=0nAkf(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=0nAkPr(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=0nAkPr+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 bax0<x1<<xnb,总存在相应的求积系数 A0,A1,⋯ ,AnA_0,A_1,\cdots,A_nA0A1,,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=ba=21(b2a2)=31(b3a3)=n+11(bn+1an+1)

3. Newton-Cotes求积公式

将区间 [a,b] n[a,b]\ n[a,b] n 等分,步长 h=b−anh=\dfrac{b-a}{n}h=nba.

易知,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)dx2ba[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)dx6ba[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=90ba[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+i4ba(i=0,1,2,3,4).

4. 复化求积公式

将积分区间 [a,b][a,b][a,b] 划分为 nnn 等分,步长 h=b−anh=\dfrac{b-a}nh=nba ,分点为 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=0n1Ii 作为所求积分 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=0n12h[f(xi)+f(xi+1)]=2h[f(a)+2i=1n1f(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]=ITn=12bah2f′′(η),η(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=0n16h[f(xi)+4f(xi+21)+f(xi+1)]=6h[f(a)+4i=0n1f(xi+21)+2i=1nf(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]=ISn=2880bah4f(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=34T2n31TnCn=1516S2n151SnRn=6364C2n631Cn

其中,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_nTnSnS_nSnCnC_nCnRnR_nRn
T1T_1T1
T2T_2T2S1S_1S1
T4T_4T4S2S_2S2C1C_1C1
T8T_8T8S4S_4S4C2C_2C2R1R_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=0nAjf(xj)+E(f)

具有 2n+12n+12n+1 次的代数精度,则称此求积公式为 Gauss 型求积公式,相应的求积节点 xjx_jxj 称为 Gauss 点。

Gauss_Legendre 求积公式:

节点数xjx_jxjAjA_jAj
111000222
222±0.5773502692\pm0.5773502692±0.5773502692111
333±0.77459666920\pm0.7745966692\\0±0.774596669200.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(xh)+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(xh)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(xh)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),axby(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个

红包金额最低5元

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

抵扣说明:

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

余额充值