研一课程笔记-数值分析4-6

非线性方程迭代解法

n次方程根的个数与次数相同(实根或者复根)

设有非线性方程

f(x)=0(4.1)f(x)=0\tag{4.1}

其中 f(x)f(x) 是一元非线性函数。若常数 ss 使 f(s)=0f(s)=0,则称 ss 是方程(4.1)的根,又称 ss 是函数 f(x)f(x)的零点。若 f(x)f( x) 能分解为

f(x)=(xs)mφ(x)f(x)=(x-s)^m\varphi(x)

其中 mm 是正整数,φ(s)0\varphi(s)\neq0 ,则称 ss 是方程(4.1)的 mm 重根和 f(x)f(x)mm 重零点。当 m=1m=1 时,ss 称为方程(4.1)的单根和 f(x)f(x) 的单零点。

二分法

C[a,b]C[a,b]C(a,b)C(a,b) 分别表示在闭区间 [a,b][a,b] 和开区间 (a,b)(a,b) 上连续的函数集合。

f(x)C[a,b]f(x)\in C[a,b]f(a)f(b)<0f(a)f(b)<0,根据连续函数的介值定理,在区间(a,b)(a,b)内至少有一实数 ss,使 f(s)=0f(s)=0。现假定在(a,b)(a,b)内只有一个实数 ss 使 f(s)=0f(s)=0,并要把 ss 求出来。用对分法求这个 s 的算法如下:

a0=a,b0=ba_0=a,b_0=b

对于k=0,1,,Mk=0,1,\cdotp\cdotp\cdotp,M 执行

(1)计算 xk=ak+bk2x_k=\frac {a_k+b_k}2

(2) 若 bkakεb_k-a_k\leqslant\varepsilonf(xk)η|f(x_k)|\leqslant\eta,则停止计算,取 sxks\approx x_k 否则转(3)。

(3)若 f(ak)f(xk)<0f(a_k)f(x_k)<0,则令 ak+1=ak,bk+1=xka_k+1=a_k,b_{k+1}=x_kf(ak)f(xk)>0f(a_k)f(x_k)>0,则令 ak+1=xk,bk+1=bka_k+1=x_k,b_{k+1}=b_k

(4)若k=Mk=M,则输出MM次迭代不成功的信息;否则继续。

迭代终止条件

一般通过 xk+1xk<ε1\begin{vmatrix}x_{k+1}-x_k\end{vmatrix}<\varepsilon_1 或者 f(x)<ε2\begin{vmatrix}f(x)\end{vmatrix}<\varepsilon_2 来作为迭代的终止条件,一般我们选择前者,因为后者可能会有梯度非常小的情况

简单迭代法及其收敛性

f(x)=0f(x)=0 转化为等价的方程

x=φ(x)(4.2)x=\varphi(x) \tag{4.2}

从而 f(s)=0f( s) = 0 等价于 s=φ(s)s= \varphi ( s)ssφ(x)\varphi(x) 的不动点。

选定ss的初始近似值x0x_0,用递推公式

xk+1=φ(xk)(k=0,1,)(4.3)x_{k+1}=\varphi(x_k)\quad(k=0,1,\cdotp\cdotp\cdotp)\tag{4.3}

产生序列 {xk}\left\{x_k\right\} ,在一定条件下,序列 {xk}\left\{x_k\right\} 收敛于 ss。在 {xk}\left\{x_k\right\} 收敛的情况下,当 kk 足够大时就可取 xkx_k 作为方程(4.1)的近似根。迭代公式(4.3)被称为求解方程(4.1)的简单迭代法,其中 φ(x)\varphi(x)称为迭代函数。

对于任意的式(4.1)可构成多种简单选代法,它们的迭代函数各不相同。为求同一个根,它们所产生的序列{xk}\left\{x_k\right\},有的可能收敛,有的可能不收敛;有的会收敛得快,有的会收敛得慢。这一切都取决于送代函数 φ(x)\varphi(x)在有根区间内的性态。

大范围收敛性定理

设函数 φ(x)C[a,b]\varphi(x)\in C[a,b],在(a,b)a,b)内可导,且满足两个条件:

  1. x[a,b]x\in[a,b]时,φ(x)[a,b];\varphi(x)\in[a,b];
  2. x(a,b)x\in(a,b)时,φ(x)L<1|\varphi^\prime(x)|\leqslant L<1,其中LL 为一常数。

则有如下结论:

  1. 方程(4.2)在区间[a,b]上有唯一的根 ss ;
  2. 对任取的 x0[a,b]x_0\in[a,b],简单迭代法(4. 3)产生的序列 {xk}[a,b]\{x_k\}\subset[a,b] 且收敛于 ss;
  3. 成立误差估计式

sxkLk1Lx1x0(4.4)\mid s-x_k\mid\leqslant\frac{L^k}{1-L}\mid x_1-x_0\mid \tag{4.4}

sxkL1Lxkxk1(4.5)\mid s-x_k\mid\leqslant\frac L{1-L}\mid x_k-x_{k-1}\mid \tag{4.5}

1)令 F(x)=xφ(x)F(x)=x-\varphi(x),则 F(x)C[a,b]F(x)\in C[a,b],并由条件(1)可知

F(a)=aφ(a)0,F(b)=bφ(b)0F(a)=a-\varphi(a)\leqslant0,\quad F(b)=b-\varphi(b)\geqslant0

若上面两个不等式中有一个等号成立,则方程(4.2)有根 s=as=as=bs=b。若两个都是严格不等式,则根据连续函数的介值定理,必存在 s(a,b)s\in(a,b),使 F(s)=sφ(s)=0F(s)=s-\varphi(s)=0,即方程(4.2)有根s(a,b)s\in(a,b)。今设有两个不同的 s1,s2[a,b]s_1,s_2\in[a,b]使 s1=φ(s1),s2=φ(s2)s_1=\varphi(s_1),s_2=\varphi(s_2),则由微分中值定理以及条件(2),有

s1s2=φ(s1)φ(s2)=φ(ξ)s1s2Ls1s2<s1s2\mid s_1-s_2\mid=\mid\varphi(s_1)-\varphi(s_2)\mid=\mid\varphi^{\prime}(\xi)\mid\mid s_1-s_2\mid\leqslant L\mid s_1-s_2\mid<\mid s_1-s_2\mid

其中 ξ\xis1s_1s2s_2 之间,因而 ξ(a,b)\xi\in(a,b)。上式出现的矛盾证实s1=s2s_1=s_2


2)因 x0[a,b]x_0\in[a,b],由条件(1)可知,{xk}[a,b]\{x_k\}\subset[a,b]。又由条件(2),得

xks=φ(xk1)φ(s)=φ(ξk)xk1sLxk1sLkx0s\mid x_k-s\mid=\mid\varphi(x_{k-1})-\varphi(s)\mid=\mid\varphi^{\prime}(\xi_k)\mid\cdot\mid x_{k-1}-s\mid\leqslant L\mid x_{k-1}-s\mid\leqslant\cdotp\cdotp\cdotp\leqslant L^k\mid x_0-s\mid

其中 ξk\xi_kxk1x_k-1ss 之间,因而 ξk(a,b)\xi_k\in(a,b)。因 0L<1\leqslant L<1,故有 limkxk=s\lim_{k\to\infty} x_k=s


3)设 m>km>k,则有

xmxk=i=km1(xi+1xi)x_m-x_k=\sum_{i=k}^{m-1}(x_{i+1}-x_i)

xi+1xi=φ(xi)φ(xi1)Lxixi1Lix1x0\begin{aligned}\mid x_{i+1}-x_i\mid&=\mid\varphi(x_i)-\varphi(x_{i-1})\mid\leqslant\\&L\mid x_i-x_{i-1}\mid\leqslant\ldots\leqslant L^i\mid x_1-x_0\mid\end{aligned}

于是有:

xmxki=km1xi+1xii=km1Lix1x0=Lk1Lmk1Lx1x0\mid x_m-x_k\mid\leqslant\sum_{i=k}^{m-1}\mid x_{i+1}-x_i\mid\leqslant\sum_{i=k}^{m-1}L^i\mid x_1-x_0\mid=L^k\frac{1-L^{m-k}}{1-L}\mid x_1-x_0\mid

xmxk=xmxm1+xm1xkx_m-x_k=x_m-x_{m-1}+x_{m-1}\cdots-x_k

然后使用绝对值不等式

mm\to\infty,由于 0L<1\leqslant L<1,故由上式得到式(4.4)。

又由

xi+1xiLxixi1Lik+1xkxk1\mid x_{i+1}-x_i\mid\leqslant L\mid x_i-x_{i-1}\mid\leqslant\cdots\leqslant L^{i-k+1}\mid x_k-x_{k-1}\mid

xmxki=km1xi+1xii=km1Lik+1xkxk1=L1Lmk1Lxkxk1\mid x_m-x_k\mid\leqslant\sum_{i=k}^{m-1}\mid x_{i+1}-x_i\mid\leqslant\sum_{i=k}^{m-1}L^{i-k+1}\mid x_k-x_{k-1}\mid=L\frac{1-L^{m-k}}{1-L}\mid x_k-x_{k-1}\mid

mm\to\infty,由上式得到式(4.5)。证毕。

局部收敛性定理

注:定理4.2并没有给出邻域 δ\delta 具体是什么,怎么计算,仅仅要求该邻域存在,也即迭代的起始点 x0x_0 足够接近 ss 时 。总能得到收敛于 ss 的迭代序列。

简单迭代法收敛速度

这里 θ\theta 相当于一个仿射变换,此外证明得到的 limkφ(sθek)=φ(s)\lim_{k\to\infty}\mid\varphi^{\prime}(s-\theta e_{k})\mid=\mid\varphi^{\prime}(s)\mid ,也就是说常数 c=φ(s)c=\mid\varphi^{\prime}(s)\mid

证明中一个非常关键的点是连续性

牛顿法

用简单送代法求方程(4.1)的根 ss,十分重要的问题是构造迭代函数 φ(x)\varphi(x)。为了使收敛速度的阶高一些,应尽可能使 φ(x)\varphi(x)x=sx=s 处有更多阶导数等于零。

现在令 φ(x)=x+h(x)f(x)\varphi(x)=x+h(x)f(x)h(x)h(x) 为待定函数,但 h(s)0h(s)\neq0,代入简单迭代法 φ(x)=x\varphi(x)=x,则方程(4.1)与方程

x=x+h(x)f(x)x=x+h(x)f(x)

有共同的根 ss。现用条件 φ(s)=0\varphi^{\prime}(s)=0 确定 h(x)h(x)。由 f(s)=0f(s)=0,则:

φ(s)=1+h(s)f(s)+h(s)f(s)=1+h(s)f(s)=0\varphi^\prime(s)=1+h^\prime(s)f(s)+h(s)f^\prime(s)=1+h(s)f^\prime(s)=0

h(x)h(x) 必须满足 h(s)=1f(s)h(s)=\frac{-1}{f^{\prime}(s)}。显然,取 h(x)=1f(x)h(x)=-\frac1{f^{\prime}(x)} 就具备这个条件,并且也满反 h(s)0h(s)\neq0。于是,φ(x)\varphi(x) 被确定为

φ(x)=xf(x)f(x)\varphi(x)=x-\frac{f(x)}{f^{\prime}(x)}

它满足 φ(s)=0\varphi^\prime(s)=0。由此得出下面的特殊的简单迭代法

xk+1=xkf(xk)f(xk)(k=0,1,)(4.7)x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime}(x_k)}\quad(k=0,1,\cdotp\cdotp\cdotp)\tag{4.7}

式(4.7)所表示的迭代法称为 Newton(牛顿)法。

Newton 法可求方程(4.1)的实数根和复数根。当求实数根时,Newton 法有明显的几何意义。当获得 xkx_k 之后,过曲线 y=f(x)=f(x)上的点(xk,f(xk))(x_k,f(x_k))作该曲线的切线,此切线与 xx 轴相交的交点横坐标就是 Newton 法迭代序列的第k+1k+1个元素xk1x_{k-1}。因此 Newton 法(4.7)又称为切线法。

局部收敛性

定理 4.6 设 ss 是方程(4.1)的根,在包含 s 的某个开区间内 f(x)f^\prime(x) 连续且 f(x)0f^{\prime}(x)\neq0,则存在 δ>0\delta{>}0,当 x0[sδ,s+δ]x_0\in[s-\delta,s+\delta]时,由 Newton 法(4.7)产生的序列 {xk}\{x_k\} 收敛于 ss;若 f(s)0f^\prime(s)\neq0x0sx_0\neq s,则序列 {xk}\left\{x_k\right\} 是平方收敛的。


Newton 法(4.7)的迭代函数 φ(x)\varphi(x)

φ(x)=xf(x)f(x),φ(x)=f(x)f(x)[f(x)]2\varphi(x)=x-\frac{f(x)}{f^{'}(x)},\quad\varphi^{'}(x)=\frac{f(x)f^{''}(x)}{[f^{'}(x)]^2}

f(x)f^\prime(x)连续且 f(x)0f^\prime(x)\neq0,故 φ(x)\varphi^{\prime}(x)在包含 ss 的某个开区间内连续,并且有

φ(s)=s,φ(s)=0\varphi(s)=s,\varphi^\prime(s)=0

根据定理 4.2,必存在 δ>0\delta>0,当 x0[sδ,s+δ]x_0\in[s-\delta,s+\delta]时,由送代法(4.7)产生的序列 {xk}[sδ,s+δ]\{x_k\}\subset[s-\delta,s+\delta],且收敛于 ss

由 Taylor 公式以及 f(s)=0f(s)=0,有

f(s)=f(xk)+f(xk)(sxk)+f(ξk)2!(sxk)2=0f(s)=f(x_k)+f^{\prime}(x_k)(s-x_k)+\frac{f^{\prime\prime}(\xi_k)}{2!}(s-x_k)^2=0

xk+1s=xkf(xk)f(xk)s=f(ξk)2f(xk)(sxk)2x_{k+1}-s=x_k-\frac{f(x_k)}{f^{\prime}(x_k)}-s=\frac{f^{\prime\prime}(\xi_k)}{2f^{\prime}(x_k)}(s-x_k)^2

其中 ξk\xi_kxkx_kss 之间。因 f(s)0f^\prime(s)\neq0 以及 f(x)f^{\prime}(x) 的连续性,故可设当 x[sδ,s+δ]x\in[s-\delta,s+\delta]f(x)0f^{\prime\prime}(x)\neq0。因而当 x0[sδ,s+δ]x_0\in[s-\delta,s+\delta]x0sx_0\neq s,xks(k=1,2,...),x_k\neq s(k=1,2,...)。于是有

limkxk+1sxks2=f(s)2f(s)>0\lim_{k\to\infty}\frac{\mid x_{k+1}-s\mid}{\mid x_k-s\mid^2}=\left|\frac{f^{\prime\prime}(s)}{2f^{\prime}(s)}\right|>0

故序列 {xk}\left\{x_k\right\} 是平方收敛的。

证毕。

定理 4.6 是 Newton 法的局部收敛性定理。如果 f(x)f(x) 只满足此定理的条件,则在一般情况下,初值 x0x_{0} 要比较靠近所要求的根 ss,Newton 法才能收敛,也就是定理结论中的 δ\delta 要比较小。

大范围收敛性(收敛的充分条件

定理 4.7 设函数 f(x)f(x) 在区间 [a,b][a,b]上存在二阶连续导数,且满足条件:

  1. f(a)f(b)<0f( a) f( b) {< }0 ;(根存在)

  2. f(x)f^{\prime\prime}(x) 在区间 [a,b][a,b] 上不变号;(根唯一)

  3. x[a,b]x\in[a,b] 时,f(x)0;f^\prime(x)\neq0;f(x)0f^\prime(x)\neq0,根唯一)

  4. x0[a,b],f(x0)f(x0)>0x_{0}\in [ a, b] , f( x_{0}) f^{\prime \prime }( x_{0}) > 0

则由 Newton 法(4.7)产生的序列{xk}\{ x_k \} 单调收敛于方程(4.1) 在 [a,b][a, b] 内唯一的根 ss,并且至少是平方收敛的。

插值与逼近

代数插值

给定 n+1n+1 个互异的实数 x0,x1,,xnx_0, x_1, \cdots , x_n,实值函数 f(x)f( x) 在包含 x0,x1,,xnx_0,x_1,\cdots,x_n 的某个区间 [a,b][a,b] 内有定义。设函数组

{φk(x)(k=0,1,,n)}\{\varphi_k(x)(k=0,1,\cdots,n)\}

是次数不高于 nn 的多项式组,且在点集 {x0,x1,,xn}\left\{x_0,x_1,\cdots,x_n\right\} 上线性无关。

现在提出如下的问题:在次数不高于 nn 的多项式集合

Dn=Span{φ0,φ1,,φn}\mathscr{D}_n=\mathrm{Span}\{\varphi_0,\varphi_1,\cdots,\varphi_n\}

中寻求多项式

pn(x)=k=0nckφk(x)(5.1)p_n(x)=\sum_{k=0}^nc_k\varphi_k(x)\tag{5.1}

使其满足条件

pn(xi)=f(xi)(i=0,1,,n)(5.2)p_n(x_i)=f(x_i)\quad(i=0,1,\cdots,n)\tag{5.2}

此问题称为一元函数的代数插值问题。x0,x1,,xnx_0,x_1,\cdotp\cdotp\cdotp,x_n 称为插值节点;f(x)f(x) 称为被插值函数; φk(x)(k=0,1,,n)\varphi_k(x)(k=0,1,\cdots,n) 称为插值基函数;条件(5.2)称为插值条件;满足插值条件(5.2)的多项式(5.1)称为 n 次插值多项式。

由于插值基函数组 {φk(x)(k=0,1,,n)}\{\varphi_k(x)(k=0,1,\cdots,n)\} 在点集 {x0,x1,,xn}\left\{x_0,x_1,\cdots,x_n\right\} 上线性无关,所以满足插值条件(5.2)的 nn 次插值多项式 pn(x)p_n(x) 是存在且唯一的。

又由于插值基函数组限定为次数不高于nn的多项式组,所以对于不同的插值基函数组,只要满足同一插值条件(5.2),则所得的nn次插值多项式也是唯一的。

对于 n+1n+1 个插值点, pn(xi)=f(xi)p_n(x_i)=f(x_i)n+1n+1 个方程:p(xi)=j=0najxij=f(xi)p(x_i)=\sum^n_{j=0}a_jx_i^j=f(x_i)

最后可以化为矩阵 XA=FX\cdot A=F (分别是 列向量 xix_i,值 ai,f(xi)a_i,f(x_i) 构成的矩阵)

拉格朗日插值

取插值基函数

lk(x)=j=0nxxjxkxj(k=0,1,,n)(5.3)l_k(x)=\prod_{j=0}^n\frac{x-x_j}{x_k-x_j}\quad(k=0,1,\cdotp\cdotp\cdotp,n) \tag{5.3}

显然,lk(x)(k=0,1,...,n)l_k(x)(k=0,1,...,n) 都是 nn 次多项式,且具有下列性质

lk(xi)={1,i=k0,ikl_k(x_i)=\begin{cases}1,&i=k\\0,&i\neq k\end{cases}

显然,通过这种方式得到的 lk(xi)l_k(x_i) 仅第 ii 个分量为1,所以彼此一定会线性无关

因此,函数组 {lk(x)(k=0,1,,n)}\left\{l_k(x)(k=0,1,\cdotp\cdotp\cdotp,n)\right\} 必在点集 {x0,x1,,xn}\left\{x_0,x_1,\cdotp\cdotp\cdotp,x_n\right\} 上线性无关,并且

pn(x)=k=0nlk(x)f(xk)(5.4)p_n(x)=\sum_{k=0}^nl_k(x)f(x_k)\tag{5.4}

ckl(xk)=0,kjc_kl(x_k)=0,k\neq j

pn(xj)=cil(xi)=cjl(xj)=cj=f(xj)p_n(x_j)=\sum c_i l(x_i)=c_jl(x_j)=c_j=f(x_j)

pn(x)=k=0n(j=0nxxjxkxj)f(xk)(5.5)p_n(x)=\sum_{k=0}^n\Big(\prod_{j=0}^n\frac{x-x_j}{x_k-x_j}\Big)f(x_k)\tag{5.5}

就是满足插值条件(5.2)的 nn 次插值多项式。

牛顿插值

拉格朗日插值在计算过程中若需要改变节点或增加节点时, 全部基函数 lk(x)l_k(x) 都需要重新计算 ,这在实际计算中很不方便

目标是寻找新的基函数,使得每加一个节点时,只附加一项 上去 ,而不必重新计算全部基函数

Newton 插值多项式的构造

基函数如下:

{φk(x)}={1,xx0,(xx0)(xx1),,(xx0)(xx1)(xxn1)}\{\varphi_k(x)\}=\{1,x-x_0,(x-x_0)(x-x_1),\ldots,(x-x_0)(x-x_1)\ldots(x-x_{n-1})\}

(φ0(x0),φ1(x0),,φn(x0)φ0(x1),φ1(x1),,φn(x1)φ0(xn),φ1(xn),,φn(xn))=(10001x1x0001x2x0(x2x0)(x2x1)01xnx0(xnx0)(xnx1)j=0n1(xnxj))\begin{pmatrix}\varphi_0(x_0),\varphi_1(x_0),\ldots,\varphi_n(x_0)\\\varphi_0(x_1),\varphi_1(x_1),\ldots,\varphi_n(x_1)\\\vdots&\vdots&\vdots\\\varphi_0(x_n),\varphi_1(x_n),\ldots,\varphi_n(x_n)\end{pmatrix}=\begin{pmatrix}1&0&0&0\\1&x_1-x_0&0&0\\1&x_2-x_0&(x_2-x_0)(x_2-x_1)&\vdots\\\vdots&\vdots&\vdots&0\\\\1&x_n-x_0&(x_n-x_0)(x_n-x_1)&\prod_{j=0}^{n-1}(x_n-x_j)\end{pmatrix}

所以这组函数线性无关,故插值函数为:

pn(x)=c0+c1(xx0)+c2(xx0)(xx1)++cn(xx0)(xx1)(xxn1)\begin{aligned}p_{n}(x)&=c_0+c_1(x-x_0)+c_2(x-x_0)(x-x_1)\\&+\cdots+c_n(x-x_0)(x-x_1)\cdots(x-x_{n-1})\end{aligned}

pn(xi)=f(xi)(i=0,1,,n)p_n(x_i)=f(x_i)\quad(i=0,1,\cdots,n)

{f(x0)=pn(x0)=c0f(x1)=pn(x1)=c0+c1(x1x0)f(x2)=pn(x2)=c0+c1(x2x0)+c2(x2x0)(x2x1)f(x3)=pn(x3)=c0+c1(x3x0)+c2(x3x0)(x3x1)f(xn)=pn(xn)=c0+c1(xnx0)+c2(xnx0)(xnx1)++cn(xnx0)(xnxn1)\begin{cases}&f(x_0)=p_n(x_0)=c_0\\&f(x_1)=p_n(x_1)=c_0+c_1(x_1-x_0)\\&f(x_2)=p_n(x_2)=c_0+c_1(x_2-x_0)+c_2(x_2-x_0)(x_2-x_1)\\&f(x_3)=p_n(x_3)=c_0+c_1(x_3-x_0)+c_2(x_3-x_0)(x_3-x_1)\\&f(x_n)=p_n(x_n)=c_0+c_1(x_n-x_0)+c_2(x_n-x_0)(x_n-x_1)+\cdots+c_n(x_n-x_0)\cdots(x_n-x_{n-1})\end{cases}

可得 ck=f[x0,x1,,xk]c_k=f[x_0,x_1,\cdots,x_k]

k 阶差商

一阶差商:f[x0,x1]=f(x1)f(x0)x1x0f[x_0,x_1]=\frac{f(x_1)-f(x_0)}{x_1-x_0}

二阶差商:f[x0,x1,x2]=f[x0,x2]f[x0,x1]x2x1f[x_0,x_1,x_2]=\frac{f[x_0,x_2]-f[x_0,x_1]}{x_2-x_1}

k阶插商:

f[x0,x1,,xk]=f[x0,x1,,xk2,xk]f[x0,x1,,xk1]xkxk1f[x_0,x_1,\cdots,x_k]=\frac{f[x_0,x_1,\cdots,x_{k-2},x_k]-f[x_0,x_1,\cdots,x_{k-1}]}{x_k-x_{k-1}}

牛顿插值公式利用差商可简单地表为

N0(x)=f(x0)N1(x)=f(x0)+f[x0,x1](xx0)=N0(x)+f[x0,x1](xx0)N2(x)=N1(x)+f[x0,x1,x2](xx0)(xx1)......Nk+1(x)=Nk(x)+f[x0,,xk+1](xx0)(xx1)(xxk)\begin{aligned} &N_0(x)=f(x_0) \\ &N_1(x)=f(x_0)+f[x_0,x_1](x-x_0) \\ & \quad \quad \quad =N_0(x)+f[x_0,x_1](x-x_0) \\ &N_2(x)=N_1(x)+f[x_0,x_1,x_2](x-x_0)(x-x_1) \\ &...... \\ &N_{k+1}(x)=N_k(x)+\underbrace{f[x_0,\cdots,x_{k+1}](x-x_0)(x-x_1)\cdots(x-x_k)} \end{aligned}

因此, 每增加一个结点, Newton 插值多项式只增加一项(即括号标注的一项)

插值余项

x0,x1,,xnx_0, x_1, \cdots , x_n 是互异实数,对给定的实数xx,实值函数f(x)f( x) 在区间 IxI_x 上有 n+1n+ 1 阶导数,则拉格朗日插值多项式与牛顿插值多项式有相同的插值余项:

Rn(x)=f(x)pn(x)=f(n+1)(ξ)(n+1)!i=0n(xxi).R_n(x)=f(x)-p_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^n(x-x_i).

罗尔定理的推论:若 φ(x)\varphi ( x) 充分光滑,且 ϕ(x1)==ϕ(xn)=0\phi ( x_1) = \cdots = \phi ( x_n) = 0\Rightarrow 存在 nln-lξI~\xi \in \tilde{I} 使得 φ(l)(ξ)=0\varphi ^{( l) }( \xi ) = 0

此外,过 n+1n +1 个点的 nn 次多项式是唯一的:Nn(x)Ln(x)N_n(x)\equiv L_n(x)

Rn(x)R_n(x)n+1n+1 个根 x0x_0~xnx_n

证明:

由于 Rn(xi)=0R_n( x_i) = 0 , i=0,1,...,n\boldsymbol{i}= 0, 1, . . . , n \Rightarrow Rn(x)=u(x)i=0(xxi)=f(x)pn(x)R_{n}(x)=u(x)\prod_{i=0}(x-x_{i})=f(x)-p_{n}(x)

任意固定xxix\neq x_i (i=0,...,n)( i= 0, . . . , n),构造辅助函数

g(t)=Rn(t)u(x)i=0n(txi)=f(t)pn(t)u(x)i=0n(txi)(1)g(t)=R_n(t)-u(x)\prod_{i=0}^n(t-x_i)=f(t)-p_n(t)-u(x)\prod_{i=0}^n(t-x_i)\tag{1}

g(t)g( t)n+2n+ 2 个不同的根 x,x0,,xnx, x_0, \ldots , x_n 反复使用Rolle定理,可得存在 ξI~x\xi\in\tilde{I}_x,使得:

g(n+1)(ξ)=0,ξI~xg^{(n+1)}(\xi)=0,\xi\in\tilde{I}_{x}

u(x)=f(n+1)(ξ)(n+1)!Rn(x)=f(n+1)(ξ)(n+1)!i=0n(xxi).(2)\Rightarrow u(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\Rightarrow R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^n(x-x_i).\tag{2}

显然 ξ\xixx 有关

对 (1) 式求导,由于 pn(t)p_n(t)nn 次多项式,所以 pn(n+1)(t)=0p_n^{(n+1)}(t)=0[u(x)i=0n(txi)](n+1)=u(x)(n+1)![u(x)\prod_{i=0}^n(t-x_i)]^{(n+1)}=u(x)(n+1)! ,所以:

g(n+1)(t)=f(n+1)(t)u(x)(n+1)!g^{(n+1)}(t)=f^{(n+1)}(t)-u(x)(n+1)!

即可得到 (2) 式

正交多项式

权函数

定义 权函数 在区间(a,b)a,b)上,若非负函数ρ(x)\rho ( x) 满足

  1. 对一切整数n0,abxnρ(x)dxn\geq0,\int_a^bx^n\rho(x)dx 存在;

  2. (a,b)( a, b) 上的非负连续函数 f(x)f( x),若 abρ(x)f(x)dx=0\int _{a}^{b}\rho ( x) f( x) dx= 0,则在区间 (a,b)( a, b)f(x)0.f( x) \equiv 0.

那么,称 ρ(x)\rho(x)(a,b)(a,b) 上的权函数。

常见的权函数:ρ(x)1,axb;\rho ( x) \equiv 1, a\leq x\leq b;

内积

给定 f(x),g(x)C[a,b],ρ(x)f( x) , g( x) \in C[ a, b] , \rho ( x)(a,b)(a,b) 上的权函数,称

(f,g)=abρ(x)f(x)g(x)dx(f,g)=\int_a^b\rho(x)f(x)g(x)dx

ff,gg[a,b][a,b] 上的内积.

内积的性质
  1. (f,g)=(g,f);\left(f,\:g\right)=\left(\:g,\:f\right);

  2. (kf,g)=(f,kg)=k(f,g),k(kf,g)=(f,kg)=k(f,g),k 为常数,kR;k\in R;

  3. (f1+f2,g)=(f1,g)+(f2,g);\left(f_{1}+f_{2},\:g\right)=\left(f_{1},g\:\right)+\left(f_{2},g\:\right);

  4. 若在 [a,b][a,b]f(x)≢0f(x)\not\equiv0,则 (f,f)>0.(f,f)>0.

正交与正交函数系

若内积

(f,g)=abρ(x)f(x)g(x)dx=0(f,g)=\int_a^b\rho(x)f(x)g(x)dx=0

则称 ff, g 在 [a,b][a,b] 上带权 ρ(x)\rho(x) 正交.

正交函数系

定义 若函数系 {φ0(x),φ1(x),...,φn(x),...}\{\varphi_0(x),\varphi_1(x),...,\varphi_n(x),...\} 满足

(φi,φj)=abρ(x)φi(x)φj(x)dx={0,ijai>0,i=j(\varphi_i,\varphi_j)=\int_a^b\rho(x)\varphi_i(x)\varphi_j(x)dx=\begin{cases}0, i\neq j\\a_i>0,i=j\end{cases}

则称 {φn(x)}\{\varphi_n(x)\}[a,b][a,b] 上带权 ρ(x)\rho(x) 的正交函数系.

  • [a,b][a,b] 上带权 ρ(x)\rho(x) 的正交函数系必是线性无关的函数系,而不论 ρ(x)\rho(x) 是什么函数.

    • 这里函数系线性无关与前面插值里的线性无关不是一个概念
  • φk(x)\varphi_k(x) 是最高项系数不为零的 kk 次多项式,k=0,1,2,;k=0,1,2,\ldots;{φk(x)}\{\boldsymbol{\varphi}_k(x)\}[a,b][a,b] 上带权 ρ(x)\rho(x)的正交多项式系

    \Leftrightarrow(充分必要条件)

    对任何次数不高于 k1k-1 的多项式 q(x)q(x),总有

    (q,φk)=aρ(x)q(x)φk(x)dx=0,k=1,2,.(q,\varphi_k)=\int_a\rho(x)q(x)\varphi_k(x)dx=0,k=1,2,\cdots.

性质

  1. (结构)设 {φk(x)}\{\varphi_k(x)\}[a,b][a,b] 上带权 ρ(x)\rho(x) 的正交多项式系,则{ckφk(x)}\left\{c_k\varphi_k(x)\right\}也是[a,b][a,b]上带权 ρ(x)\rho(x)的正交多项式系。其中,ckc_k是非零常数.

  2. (唯一性){φk(x)}\varphi_k(x)\}是区间[a,b][a,b]上带权 ρ(x)\rho(x)的正交多项式系,若各个φk(x)\varphi_k(x)的最高次项系数是1,则这样的{φk(x)}\{\varphi_k(x)\}是唯一的.
    证 设 {φk(x)},{gk(x)}\{\varphi_k(x)\},\{g_k(x)\} 为性质2中要求正交多项式系,
    k=0k=0时,φ0(x)1g0(x).\varphi_0(x)\equiv1\equiv g_0(x).
    k1k\geq 1 时,gkφkg_k- \varphi _k 是次数不高于 k1k- 1的多项式,由定理5.6知,

    (gkφk,φk)=(gkφk,gk)=0(gkφk,gkφk)=0,gkφk0,gk(x)φk(x),axb\begin{aligned}&(\:g_k-\varphi_k,\:\varphi_k\:)=(\:g_k-\varphi_k,\:g_k\:)=0\\&\Leftrightarrow(g_k-\varphi_k,g_k-\varphi_k)=0,\\&\Leftrightarrow g_k-\varphi_k\equiv0,\quad g_k(x)\equiv\varphi_k(x)_,a\leq x\leq b\end{aligned}

  3. (正交多项式的根)若 φk(x){\varphi_k(x)}[a,b][a,b] 上带权 ρ(x)\rho(x) 的 正交多项式系,则 k1k\geq1时,kk 次正交多项式 φk(x)\varphi_k(x)(a,b)(a,b) 内有 kk 个互异的实根.

    事实上,可以证明 kk 次正交多项式有kk 个单根.

  4. 设{φk(x)}\varphi_k(x)\}是区间[a,b][a,b]上带权 ρ(x)\rho(x)的正交多项式系,则 k1k\geq1时,有如下的递推关系式:

ϕk+1(x)=ak+1ak(xβk)ϕk(x)ak+1ak1ak2λk1ϕk1(x)\phi_{k+1}(x)=\frac{a_{k+1}}{a_k}(x-\beta_k)\phi_k(x)-\frac{a_{k+1}a_{k-1}}{a_k^2}\lambda_{k-1}\phi_{k-1}(x)

其中,aka_kφk(x)\varphi_k(x) 的最高次项系数,且

βk=(xφk,φk)(φk,φk),λk1=(φk,φk)(φk1,φk1).\beta_{k}=\frac{(x\varphi_{k},\varphi_{k})}{(\varphi_{k},\varphi_{k})},\quad\lambda_{k-1}=\frac{(\varphi_{k},\varphi_{k})}{(\varphi_{k-1},\varphi_{k-1})}.

函数的最佳平方逼近

函数的最佳逼近问题:对函数类 AA 中给定的函数 f(x)f( x),需要在另一类较简单的便于计算的函数类 BB (BA)( B\in A) 中,找一个函数P(x)P( x),使 P(x)P(x)f(x)f (x) 之差在某种度量意义下达到最小。

函数f(x)f(x) 称为被逼近函数;P(x)P( x) 称为逼近函数,两者之差称为逼近误差。

函数类 AA 通常是区间 [a,b][a,b] 上的连续函数,记作 C[a,b]C[a,b];函数类 BB 通常是代数多项式,有理多项式 ,三角多项式,分段多项式等容易计算的函数。

常用度量标准

  1. 一致逼近(均匀逼近)以 maxaxbf(x)p(x)\max _a\leq x\leq b| f( x) - p( x) | 作为度量误差 f(x)P(x)f(x)-P( x) 的“大小” 标准

  2. 平方逼近(均方逼近)以 以 ab[f(x)p(x)]2dx\int_a^b[f(x)-p(x)]^2dx 作为度量误差 f(x)P(x)f(x)-P( x) 的“大小” 标准

最佳平方逼近的概念与解法

Hn=Span{φ0(x),φ1(x),...,φn(x)}H_n= \operatorname { Span} \{ \varphi _0( x) , \varphi _1( x) , . . . , \varphi _n( x) \}φi(x)\varphi_i(x) 张成的线性空间

定理5.7 设f(x)C[a,b],p(x)Hnf(x)\in C[a,b],p^*(x)\in H_n,在 HnH_n 中,p(x)p^* ( x)f(x)f( x) 最佳平方逼近的函数(fp,φj)=0,j=0,1,...,n\Leftrightarrow(f-p^*,\varphi_j)=0,j=0,1,...,n 其中,{φ0(x),φ1(x),...,φn(x)}\{ \varphi _0( x) , \varphi _1( x) , . . . , \varphi _n( x) \} 为子空间 HnH_n 的一组基

()( \Rightarrow )p(x)=ciφi(x)p^*(x)=\sum c_i^*\varphi _i(x) 代入 (fp,φj)=0(f-p^*,\varphi_j)=0(fciφi,φj)=0(f-\sum c_i^*\varphi _i,\varphi_j)=0

(f,φj)=(ciφi,φj)=ci(φi,φj)(f,\varphi_j)=(\sum c_i^*\varphi _i,\varphi_j)=\sum c_i^*(\varphi_i,\varphi_j)

即:

[(φ0,φ0)(φ0,φ1)(φ0,φn)(φ1,φ0)(φ1,φ1)(φ1,φn)(φn,φ0)(φn,φ1)(φn,φn)][c0c1cn]=[(f,φ0)(f,φ1)(f,φn)]\begin{bmatrix}(\varphi_0,\varphi_0)&(\varphi_0,\varphi_1)&\cdots&(\varphi_0,\varphi_n)\\(\varphi_1,\varphi_0)&(\varphi_1,\varphi_1)&\cdots&(\varphi_1,\varphi_n)\\\vdots&\vdots&\cdots&\vdots\\(\varphi_n,\varphi_0)&(\varphi_n,\varphi_1)&\cdots&(\varphi_n,\varphi_n)\end{bmatrix}\begin{bmatrix}c_0^*\\c_1^*\\\vdots\\c_n^*\end{bmatrix}=\begin{bmatrix}(f,\varphi_0)\\(f,\varphi_1)\\\vdots\\(f,\varphi_n)\end{bmatrix}

从而求解系数 cic_i^*

()( \Leftarrow )(fp,φj)=0,j=0,1,...,n( f- p^*, \varphi _j) = 0, j= 0, 1, . . . , n 成立,对任意的 p(x)Hnp( x) \in H_n,有

(fp,fp)=(fp+pp,fp+pp)=(fp,fp)+2(fp,pp)+(pp,pp)\begin{aligned}(f-p,f-p)&=(f-p^*+p^*-p,f-p^*+p^*-p)\\&=(f-p^*,f-p^*)+2(f-p^*,p^*-p)+(p^*-p,p^*-p)\end{aligned}

由于 (fp,pp)=j=0n(cjcj)(fp,φj)=0\left ( f- p^* , p^* - p\right ) = \sum _{j= 0}^n\left ( c_j^* - c_j\right ) \left ( f- p^* , \varphi _j\right ) = 0(pp,pp)0( p^* - p, p^* - p) \geq 0,故 (fp,fp)(fp,fp).(f-p,f-p)\geq(f-p^*,f-p^*).

pp^*HnH_n 中是对 ff 的最佳平方逼近函数

即设 f(x)C[a,b],p(x)Hnf(x)\in C[a,b],p^*(x)\in H_n,在 HnH_n 中,p(x)p^* ( x) 是对 f(x)f( x) 最佳平方逼近的函数\Leftrightarrow对任意的 p(x)Hnp( x) \in H_n 均有

(fp,p)=0(f-p^*,p)=0

定理5.8 设 f(x)C[a,b]f( x) \in C[ a, b] 在子空间 HnH_n中,对 f(x)f( x) 最佳平方逼近的函数是唯一的

最佳平方逼近的求解

$$ (f,g)=\int_a^b\rho(x)f(x)g(x)dx $$ 一般来说权函数默认 $\rho(x)=1$

正交函数系在最佳平方逼近中的应用

对于正交函数系,(φk,φj)=ρ(x)φk(x)φj(x)dx=0,kj(\varphi_{k},\varphi_{j})=\int\rho(x)\varphi_{k}(x)\varphi_{j}(x)dx=0,k\neq j

[(φ0,φ0)(φ0,φ1)(φ0,φn)(φ1,φ0)(φ1,φ1)(φ1,φn)(φn,φ0)(φn,φ1)(φn,φn)][c0c1cn]=[(f,φ0)(f,φ1)(f,φn)]\begin{bmatrix}(\varphi_0,\varphi_0)&(\varphi_0,\varphi_1)&\cdots&(\varphi_0,\varphi_n)\\(\varphi_1,\varphi_0)&(\varphi_1,\varphi_1)&\cdots&(\varphi_1,\varphi_n)\\\vdots&\vdots&\cdots&\vdots\\(\varphi_n,\varphi_0)&(\varphi_n,\varphi_1)&\cdots&(\varphi_n,\varphi_n)\end{bmatrix}\begin{bmatrix}c_0\\c_1\\\vdots\\c_n\end{bmatrix}=\begin{bmatrix}(f,\varphi_0)\\(f,\varphi_1)\\\vdots\\(f,\varphi_n)\end{bmatrix}

[(φ0,φ0)(φ1,φ1)(φn,φn)][c0c1cn]=[(f,φ0)(f,φ1)(f,φn)]\Rightarrow\begin{bmatrix}(\varphi_0,\varphi_0)\\&(\varphi_1,\varphi_1)\\&&\ddots\\&&&(\varphi_n,\varphi_n)\end{bmatrix}\begin{bmatrix}c_0\\c_1\\\vdots\\c_n\end{bmatrix}=\begin{bmatrix}(f,\varphi_0)\\(f,\varphi_1)\\\vdots\\(f,\varphi_n)\end{bmatrix}

ck=(f,φk)(φk,φk),k=0,1,2,...,n\Rightarrow c_k=\frac{(f,\varphi_k)}{(\varphi_k,\varphi_k)},\quad k=0,1,2,...,n

曲线拟合

与插值的一个非常大的区别就是曲线拟合的函数 y(x)y^*(x) 不要求经过所有数据点。

最小二乘问题

最小二乘问题的矩阵表示

最小二乘解的特征

最小二乘法的精度

正交多项式用于曲线拟合

数值积分

目的是求一些牛顿-莱布尼茨积分公式无法求解的积分函数

求积公式及其代数精度

一求积公式的一般形式:

abf(x)dxk=0nλkf(xk)(6.1)\int_a^bf(x)dx\approx\sum_{k=0}^n\lambda_kf(x_k)\tag{6.1}

式中 xk(k=0,1,,n)x_k( k= 0, 1, \cdots , n) 称为求积节点,它们满足 :

ax0<x1<<xnba\leq x_0<x_1<\ldots<x_n\leq b

λk(k=0,1,,n)\lambda _{k}( k= 0, 1, \cdots , n) 称为求积系数,与被积函数无关 。

当求积节点与求积系数确定了以后,求积公式就确定了,因此任务是确定求积节点和求积系数。

截断误差或余项:

Rn=abf(x)dxk=0nλkf(xk)R_n=\int_a^bf(x)dx-\sum_{k=0}^n\lambda_kf(x_k)

代数精度

对求积公式(6.1),当 f(x)f(x) 为任何次数不高于 mm 的多项式时都成立,而当 f(x)f(x) 为某个 m+1m+1 次多项式时不成立,则称它具有 mm 次代

数精度。

如何求解代数精度:对求积公式(6.1),当 f(x)f(x)1,x,x2,,xn1,x,x^2,\cdots,x^n 时,求积公式均成立,则对于任何次数不高于 mm 的多项式求积公式也必然成立。

差值型求积公式

利用插值多项式 pn(x)p_n(x) 逼近代替 f(x)f(x) 计算积分

[a,b]\left [ a, b\right ] 上取 ax0<x1<<xmba\leq x_0< x_1< \ldots < x_m\leq b 。构 造 ffnn 次Lagrange插值基函数

pn(x)=k=0nlk(x)f(xk),lk(x)=i=0nxxixkxi,k=0,1,,np_n(x)=\sum_{k=0}^nl_k(x)f(x_k),l_k(x)=\prod_{i=0}^n\frac{x-x_i}{x_k-x_i},k=0,1,\cdots,n

abf(x)dxk=0n(ablk(x)dx)f(xk)(6.1.1)\int_a^bf(x)dx\approx\sum_{k=0}^n\left(\int_a^bl_k(x)\mathrm{d}x\right)f(x_k)\tag{6.1.1}

λk(n)=ablk(x)dx=abi=0iknxxixkxidx,k=0,1,,n(6.2)\lambda_k^{(n)}=\int_a^bl_k(x)\mathrm{d}x=\int_a^b\prod_{i=0\atop i\neq k}^n\frac{x-x_i}{x_k-x_i}\mathrm{d}x,k=0,1,\cdots,n\tag{6.2}

其中,λk(n)\lambda_k^{(n)} 由节点决定与 f(x)f(x) 无关。

设函数 f(x)f(x) 在区间 [a,b][a,b] 上足够光滑,则由定理5.1得

f(x)=k=0nlk(x)f(xk)+f(n+1)(ξ)(n+1)!j=0n(xxj)f(x)=\sum_{k=0}^nl_k\left(x\right)f(x_k)+\frac{f^{(n+1)}\left(\xi\right)}{\left(n+1\right)!}\prod_{j=0}^n\left(x-x_j\right)

其中 x[a,b],ξ=ξ(x)(a,b)x\in[a,b],\xi=\xi(x)\in(a,b) 因此,可得

Rn=abf(x)dxk=0nλk(n)f(xk)=ab[f(x)Ln(x)]dx=k=0n(ablk(x)dx)f(xk)+abf(n+1)(ξ)(n+1)!j=0n(xxj)dxk=0n(ablk(x)dx)f(xk)=abf(n+1)(ξ)(n+1)!j=0n(xxj)dxR_n=\int_a^bf(x)\mathrm{d}x-\sum_{k=0}^n\lambda_k^{(n)}f(x_k)=\int_a^b[f(x)-L_n(x)]\mathrm{d}x\\=\sum_{k=0}^n\left(\int_a^bl_k\left(x\right)dx\right)f(x_k)+\int_a^b\frac{f^{(n+1)}\left(\xi\right)}{\left(n+1\right)!}\prod_{j=0}^n\left(x-x_j\right)dx-\sum_{k=0}^n\left(\int_a^bl_k^{\prime}(x)\mathrm{d}x\right)f(x_k)\\=\int_a^b \frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{j=0}^n\left(x-x_j\right)dx

定理与相关推论

  • 定理6.1 n+1n+1 个节点的求积公式(6.1.1),(6.2) 至少有 nn 次代数精度

  • 推论 n+1n+1个节点的插值型求积公式中的求积系数 λk(n)\lambda_k^{(n)} 满足

k=0nλk(n)=ba\sum_{k=0}^n\lambda_k^{(n)}=b-a

其中,a,ba,b 为积分下上限

  • 定理6.2 n+1n+1个节点的求积公式 (6.1)至少有 nn 次代数精度则它一定是插值型求积公式。

Newton-Cotes 求积公式

等距节点插值积分

将积分区间 [a,b][a,b] 划分为 nn 等份,步长为 h=banh=\frac{b-a}n ,节点等距为 xk=a+khx_k=a+kh 构造出的插值求积公式

abf(x)dxk=0n(ablk(x)dx)f(xk)\int_a^bf(x)dx\approx\sum_{k=0}^n\left(\int_a^bl_k(x)\mathrm{d}x\right)f(x_k)

λk(n)=ablk(x)dx=abi=0iknxxixkxidx,k=0,1,,n\lambda_k^{(n)}=\int_a^bl_k(x)\mathrm{d}x=\int_a^b\prod_{i=0\atop i\neq k}^n\frac{x-x_i}{x_k-x_i}\mathrm{d}x,k=0,1,\cdots,n

Newton-Cotes求积公式,式中 λk(n)\lambda_{k}^{(n)} 称为Newton-Cotes求积系数

  • 定理6.3 当 nn 为偶数时,n+1n+1个节点的 nn 阶Newton-Cotes 求积公式的代数精度至少是n+1n+1 次。

收敛性

考察是否对任何在[a,b][ a, b] 上 可 积 的 函 数 f(x)f( x) 都有

limn[k=0nλk(n)f(a+kban)]=abf(x)dx.\lim_{n\to\infty}\biggl[\sum_{k=0}^n\lambda_k^{(n)}f\biggl(a+k\frac{b-a}n\biggr)\biggr]=\int_a^bf(x)\mathrm{d}x.

其中,λk(n)\lambda _k^{( n) } 是 Newton-Cotes 求积系数。答案是否定的

几个常用公式

常微分方程初值问题的数值解法

一般概念

一常微分方程初值问题

{y=f(t,y),t0tT(7.1)y(t0)=y0(7.2)\begin{cases}y'=f\left(t,y\right),t_0\leq t\leq T&(7.1)\\y\left(t_0\right)=y_0&(7.2)\end{cases}

函数 f(x)f( x) 在区域

D0={(t,y)t0tT,y<}D_0=\left\{\left(t,y\right)|t_0\leq t\leq T,\mid y\mid<\infty\right\}

内连续目对变量 yy 满足Lipschitz(李普希兹)条件

f(t,u1)f(t,u2)Lu1u2\left|f\left(t,u_1\right)-f\left(t,u_2\right)\right|\leq L\left|u_1-u_2\right|

其中(t,u1),(t,u2)D0( t, u_1) , ( t, u_2) \in D_0 这样初值问题 (7.1),(7.2)( 7. 1) , ( 7. 2)的解存在且唯一。

额外假设:γ(t),t[t0,T]\gamma(t),t\in[t_0,T] 足够光滑

注意这里只需要对 yy 满足李普希兹条件

解的存在性

  • 定理1 设函数 f(t,y)f(t,y) 在凸集 DR2D\subset R^2 中有定义,若存在常数 LL(t,y)D\forall(t,y)\in D ,有 f(t,y)yL\left|\frac\partial f(t,y){\partial y}\right|\leq L ,则 ffDD 内关于 γ\gamma 满足Lipschtiz条件。

  • 定理2 设函数 f(t,y)f(t,y)D={(t,y)t0tT,<y<}D=\left\{(t,y)|t_0\leq t\leq T,-\infty<y<\infty\right\} 内连续且在 DD 内关于 yy 满足 Lipschitz 条件,则初值问题(7.1),(7.2)在区间 [t0,T][t_0,T] 内存在唯一解。

注:满足定理2条件的初值问题也称良态问题。

数值解法

单步法与多步法

给定步长 h>0h> 0 ,取节点 tn=t0+nh,n=0,1,,Mt_n= t_0+ nh, n= 0, 1, \cdots , M

其中 tMTt_M\leq T 通过数值方法,求出各个节点处满足 (7.1),(7.2)( 7. 1) , ( 7. 2) 的近似值 yny(tn),n=1,,My_n\approx y( t_n) , n= 1, \cdots , M 如下

F(tn,yn,yn+1,,yn+k,h)=0,n=0,1,,Mk(7.3)F( t_n, y_n, y_{n+ 1}, \cdots , y_{n+ k}, h) = 0, n= 0, 1, \cdots , M- k\tag{7.3}

上式称为关于 yn,n=0,1,,My_n, n= 0, 1, \cdots , M 的差分方程

  • k=1,F(tn,yn,yn+1,h)=0,n=0,1,,M1(7.4)k= 1, F( t_n, y_n, y_{n+ 1}, h) = 0, n= 0, 1, \cdots , M- 1 (7.4) 称之为单步法

  • k2k\geq2 则数值解法(7.3)统称为多步法,具体地称为 kk 步法

显式法与隐式法

若差分方程(7.3)能表示为如下显函数:

yn+k=G(tn,yn,yn+1,,yn+k1,h),n=0,1,,Mk(7.5)y_{n+k}=G(t_n,y_n,y_{n+1},\cdots,y_{n+k-1},h),n=0,1,\cdots,M-k\tag{7.5}

称数值解法(7.5)为显式法;否则称为隐式法(即 G()=0G(\cdots)=0 的形式)

整体截断误差

y(t)y(t) 是初值问题(7.1),(7.2)的解,y1,y2,,yMy_1,y_2,\cdots,y_M 是由数值解法(7.3)解出的初值问题(7.1),(7.2)的数值解,则称误差

εn=y(tn)yn\varepsilon_n=y\left(t_n\right)-y_n

为数值解法(7.3)在节点 tnt_n 处的整体截断误差。

显示单步法

Euler折线法

设良态微分方程初值问题

{y=f(t,y),t0tTy(t0)=y0\begin{cases}y'=f(t,y),t_0\leq t\leq T\\y(t_0)=y_0\end{cases}

解存在且唯一。取等距结点:

tn=t0+nh,n=0,1,,M;h=(Tt0)/Mt_n=t_0+nh,n=0,1,\cdots,M;h=(T-t_0)/M

由泰勒公式,有:

y(tn+1)=y(tn)+hy(tn)+h22y(ξn)=y(tn)+hf(tn,y(tn))+h22y(ξn)y(t_{n+1})=y(t_n)+hy^{\prime}(t_n)+\frac{h^2}2y^{\prime\prime}(\xi_n)\\=y(t_n)+hf(t_n,y(t_n))+\frac{h^2}2y^{\prime\prime}(\xi_n)

假定 yny(tn)y_n\approx y(t_n) ,并舍去 h2h^2 项,可构造出Eular折线法

{y0=y(t0)yn+1=yn+hf(tn,yn)n=0,1,,M\begin{cases}y_0=y(t_0)\\y_{n+1}=y_n+hf(t_n,y_n)\end{cases}\quad n=0,1,\cdots,M

几何意义

由于 f(t0,y0)f(t_0,y_0) 已知,则 y(t)y(t)(t0,y0)(t_0,y_0) 处必有切线方程,其斜率为\frac{dy}{dt}\Bigg|_{(t_0,y_0)}=f(t_0,y_0)

dydt(t0,y0)=f(t0,y0)\frac{dy}{dt}\Bigg|_{(t_0,y_0)}=f(t_0,y_0)

由点斜式写出切线方程

y=y0+(tt0)dydt(x0,y0)=y0+(tt0)f(t0,y0)y=y_0+\left(t-t_0\right)\frac{dy}{dt}\Bigg|_{(x_0,y_0)}=y_0+\left(t-t_0\right)f(t_0,y_0)

步长为 hh,则 t1t0=ht_1-t_0=h ,可由切线算出 y1=y0+hf(t0,y0)y_1=y_0+hf(t_0,y_0) 按此逐步计算 y(tn)y(t_n) ,在 tn+1t_n+1 处的值 yn+1=yn+hf(tn,yn)y_n+1=y_n+hf(t_n,y_n)

注意:y(tn)y(t_n) 为函数真值,yny_n 为迭代值

一般形式

显示单步法的一般形式

yn+1=yn+hφ(tn,yn,h),n=0,1,,M1(7.6)y_{n+1}=y_n+h\varphi\left(t_n,y_n,h\right),n=0,1,\cdots,M-1\tag{7.6}

其中,函数 φ()\varphi(\cdot) 与函数 ff 有关,称为增量函数

整体截断误差

tn+1t_{n+1} 处的整体截断误差为:

εn+1=y(tn+1)yn+1=Rn+1+εn+h[φ(tn,y(tn),h)φ(tn,yn,h)]=y(tn+1)y(tn)hφ(tn,y(tn),h)+y(tn)yn+h[φ(tn,y(tn),h)φ(tn,yn,h)](7.7)\varepsilon_{n+1}=y\left(t_{n+1}\right)-y_{n+1}\\ =R_{n+1}+\varepsilon_{n}+h\left[\varphi\left(t_n,y\left(t_n\right),h\right)-\varphi\left(t_n,y_n,h\right)\right]\\ =y(t_{n+1})-y(t_n)-h\varphi(t_n,y(t_n),h)+y\left(t_n\right)-y_n+h\left[\varphi\left(t_n,y\left(t_n\right),h\right)-\varphi\left(t_n,y_n,h\right)\right]\tag{7.7}

局部截断误差

γ(t)\gamma(t) 是初值问题(7.1),(7.2)的解,则称:

Rn+1=y(tn+1)y(tn)hφ(tn,y(tn),h)(7.8)R_{n+1}=y\left(t_{n+1}\right)-y\left(t_n\right)-h\varphi\left(t_n,y\left(t_n\right),h\right)\tag{7.8}

为单步法(7.6)的局部截断误差

定义 若单步法(7.6)的局部截断误差(7.8)与 hp+1h^{p+1}pp 为正整数)同阶,即 Rn+1=O(hp+1)R_n+1=O(h^{p+1}) 。则称单步法(7.6)是 pp 阶方法。

定理

定理 1 设增量函数 φ(t,y,h)\varphi(t,y,h) 在区域

D={(t,y,h)t0tT,y<,0hh0}D=\{(t,y,h)|t_0\leq t\leq T,|y|<\infty,0\leq h\leq h_0\}

内对变量 yy 满足Lipschitz条件,即存在常数 KK 使对 DD 内任何两点 (t,u1,h),(t,u2,h)(t,u_1,h),(t,u_2,h) ,不等式

φ(t,u1,h)φ(t,u2,h)Ku1u2\begin{vmatrix}\varphi(t,u_1,h)-\varphi(t,u_2,h)\end{vmatrix}\leq K\begin{vmatrix}u_1-u_2\end{vmatrix}

成立,那么,若单步法(7.6)的局部截断 Rn+1R_{n+1}hp+1(p1)h^{p+1}(p\geq1) 同阶,即

Rn+1=O(hp+1)R_{n+1}=O(h^{p+1})

则单步法(7.6)的整体截断误差 εn+1\varepsilon_{n+1}hph^p 同阶,即有

εn+1=O(hp)\varepsilon_{n+1}=O(h^p)

证明

由式(7.7),有

εn+1Rn+1+y(tn)yn+hφ(tn,y(tn),h)φ(tn,yn,h)Rn+1+εn+hKεn=Rn+1+εn(1+hK)Rn+1+Rn(1+hK)+εn1(1+hK)2k=0nRn+1k(1+hK)k+ε0(1+hK)n+1\left|\varepsilon_{n+1}\right|\leq \left|R_{n+1}\right|+\left|y(t_n)-y_n\right|+h\left|\varphi(t_n,y(t_n),h)-\varphi(t_n,y_n,h)\right|\\ \leq |R_{n+1}|+|\varepsilon_{n}|+hK|\varepsilon_{n}|\\ =\left|R_{n+1}\right|+\left|\varepsilon_n\right|(1+hK)\\\leq\left|R_{n+1}\right|+\left|R_n\right|(1+hK)+\left|\varepsilon_{n-1}\right|(1+hK)^2\\ \cdots\leq\sum_{k=0}^n\left|R_{n+1-k}\right|(1+hK)^k+\left|\varepsilon_0\right|(1+hK)^{n+1}

R=max1nMRnR=\max_{1\leq n\leq M}|R_n|,并注意到 ε0=0\boldsymbol{\varepsilon}_0=0 ,就有

εn+1Rk=0n(1+hK)k=RhK[(1+hK)n+11]RhK(e(n+1)hK1)1hKO(hp+1)(e(tn+1t0)K1)=O(hp)\left|\varepsilon_{n+1}\right|\leq R\sum_{k=0}^n(1+hK)^k=\frac R{hK}[(1+hK)^{n+1}-1]\leq\\\frac R{hK}(\mathrm{e}^{(n+1)hK}-1)\leq\frac1{hK}O(h^{p+1})(\mathrm{e}^{(t_{n+1}-t_0)K}-1)=O(h^p)

这里用到了 1+xex1+x\leq e^x

由此可知,εn+1\varepsilon_{n+1} 可表示为 εn+1=c(tn+1)hp+O(hp+1)\varepsilon_n+1=c(t_{n+1})h^p+O(h^{p+1})

Runge-Kutta方法

称为显式Runge-Kutta(龙格-库塔)方法,简称R-K方法, 其中正整数 NN 称为R-K方法的级,所有ci,ai,bijc_i,a_i,b_{ij} 都是待定常数。

根据定义(7.8),NN 级R-K方法(7.9)的局部截断误差为

Rn+1=y(tn+1)y(tn)hi=1Nciki(7.10)R_{n+1}=y(t_{n+1})-y(t_n)-h\sum_{i=1}^Nc_ik_i\tag{7.10}

其中 k1,k2,,kNk_1,k_2,\cdots,k_N 中的 yny_n 都换成 y(tn)y(t_n)。如果系数 ci,ai,bijc_i,a_i,b_{ij} 能使局部截断误差(7.10)与 hp+1h^p+1 同阶,则相应的 NN 级R-K方法就是 pp 阶方法。

希望:NN 确定时选择一组系数 ci,ai,bijc_i,a_i,b_{ij},使阶数 pp 达到最高

改进的Euler法

c1=c2=12,a2=1c_1=c_2=\frac12,a_2=1,式(7.14)成为

{yn+1=yn+h2(k1+k2)k1=f(tn,yn)k2=f(tn+h,yn+hk1)(7.19)\begin{cases}y_{n+1}=y_n+\frac h2(k_1+k_2)\\k_1=f(t_n,y_n)\\k_2=f(t_n+h,y_n+hk_1)\end{cases}\tag{7.19}

方法(7.19)又称为改进的Euler法。

中点公式

c1=0,c2=1,a2=12c_1=0,c_2=1,a_2=\frac12 ,则式(7.14)成为

{yn+1=yn+hk,k1=f(tn,yn)k2=f(tn+12h,yn+12hk1)(7.20)\begin{cases}y_{n+1}=y_n+hk,\\k_1=f(t_n,y_n)\\k_2=f(t_n+\frac12h,y_n+\frac12hk_1)\end{cases}\tag{7.20}

式(7.20)又称为中点公式。

Heun(休恩)公式

c1=14,c2=34,a2=23c_1=\frac14,c_2=\frac34,a_2=\frac23,则式(7.14)成为

{yn+1=yn+h4(k1+3k2)k1=f(tn,yn)k2=f(tn+23h,yn+23hk1)(7.21)\begin{aligned}&\begin{cases}y_{n+1}=y_n+\dfrac{h}{4}(k_1+3k_2)\\k_1=f(t_n,y_n)\\k_2=f(t_n+\dfrac{2}{3}h,y_n+\dfrac{2}{3}hk_1)\end{cases}\end{aligned}\tag{7.21}

式(7.21)又称为Heun(休恩)公式。

相容性、收敛性和绝对稳定性

$$ e_n=y_n-\tilde{y}_n\\ \begin{align} e_{n+1}&=y_{n+1}-\tilde{y}_{n+1}\\ &=(1+h\lambda)(y_n-\tilde{y}_n)\\ &=(1+h\lambda)e_n \end{align} $$ 当且仅当 $|1+h\lambda|<1$ 时有 $\lim_{n\to\infty} e_n=0$。所以,Euler法(7.13)的绝对稳定区域为

1+hλ<1|1+h\lambda|<1

λ\lambda 为实数时,Euler法的绝对稳定区间 2<hλ<0-2<h\lambda<0

步长选择