980 likes | 1.17k Views
数值计算方法. 第 5 章 函数的插值与 最佳平方逼近. 实践中常有这样的问题: (1) 由实验得到某一函数 f ( x ) 在一系列点 x 0 , x 1 , … , x n 处的值 f 0 , f 1 , … , f n ,其函数的解析表达式是未知的 (2) 或者 f ( x ) 虽有解析式,但计算复杂,不便于使用 需要构造一个简单函数 y ( x ) 近似地代替 f ( x ) —— 这就是函数逼近问题. 5.0 基本概念 1. 逼近函数与被逼近函数
E N D
数值计算方法 第5章 函数的插值与 最佳平方逼近
实践中常有这样的问题: • (1) 由实验得到某一函数f (x)在一系列点x0,x1,…,xn处的值f0,f1,…,fn,其函数的解析表达式是未知的 • (2) 或者f (x)虽有解析式,但计算复杂,不便于使用 • 需要构造一个简单函数y(x)近似地代替f (x) —— 这就是函数逼近问题
5.0 基本概念 • 1. 逼近函数与被逼近函数 • 函数逼近问题中的函数f (x)称为被逼近函数,y(x)称为逼近函数,其中所谓简单函数指可用四则运算进行计算的函数(如:有理、多项式、分段多项式) • 2. 逼近的度量 • (1) 以 为度量的逼近称为一致逼近 • (2) 以 为度量的逼近称为平方逼近
3. 插值与拟合 • 设已知被逼近函数f (x)在离散点xi [a,b]上的值 • f (xi) = fi, • (1) 要求y(x)满足 • (甚至 )的问题称为函数插值。 • (2) 要求y(x)满足 为最小的问题称为数据拟合(曲线拟合)
4. 简单函数类 • 设φ0,φ1,…,φn线性无关,令 • Φ = span{φ0,φ1,…,φn}为简单函数类, • 其中φ0,φ1,…,φn称为Φ的基函数。 • 逼近问题即用y(x) = a0φ0(x) + a1φ1(x) +…+ anφn(x)来做逼近,问题归结为求其中的待定系数a0,a1,…,an。
5.1 多项式插值 • 即:求多项式pn(x)满足插值条件: • pn(xi) = f(xi) = fii = 0,1,2,…,n (5.1-1) • 其中点xi [a,b] i = 0,1,2,…,n,称为插值节点,区间[a,b]称为插值区间,pn(x)称为插值多项式
定理5.1-1 存在唯一pn(x) Pn[x]满足插值条件(5.1-1) • pn(xi) = f(xi) = fii = 0,1,2,…,n • 证明:取Pn[x]的一组基{1,x,x2,…,xn },则pn(x) Pn[x]表为 • pn(x) = a0 + a1x +…+ anxn (5.1-2) • 由(5.1-1)知 • pn(xi) = a0 + a1xi +…+ anxin = fi (i = 0,1,2,…,n) (5.1-3) • (5.1-3)的系数行列式为范德蒙行列式:
(5.1-3)的系数行列式为范德蒙行列式: • 因为x0,x1,…,xn互异,所以Vn ≠ 0 • 即(5.1-3)存在唯一解,从而存在唯一的pn(x) Pn[x] 满足插值条件(5.1-1)。
例1 给定数据 • 求次数不小于3的插值多项式p3(x) • 解:设p3(x) = a0 + a1x+ a2x2 +a3x3 • 依题意有 • 解之得:a0 = 10,a1 = 5,a2 = – 5,a3= 2 • 即有p3(x) = 10 + 5x– 5x2 +2x3 • 注:(1) 范德蒙矩阵的条件数很大 —— 误差大计算量大 • (2) 选择适当基函数使插值多项式具有特殊形式
1. Lagrange插值 • 因为 • 所以先考虑特殊的插值问题。求次数不大于n的多项式li(x)满足 • (5.1-4)
由定理5.1-1知,li(x)唯一存在,且有n个零点:x0,...,xi-1,xi+1,...,xn由定理5.1-1知,li(x)唯一存在,且有n个零点:x0,...,xi-1,xi+1,...,xn • 所以li(x) = bi(x – x0)... (x – xi-1)(x – xi+1)... (x – xn) • 又由li(xi) = 1,得 • 即 • (5.1-5)
注:(1) 易知{l0,...,ln}为Pn[x]的一组基,称为 • 以x0,...,...,xn为节点的Lagrange插值基函数。 • (2) 令 • (5.1-8) • 则易知(5.1-8)所示的pn(x)为次数不大于n的多项式,且满足插值条件(5.1-1) • (j = 0,1,...,n) • 称pn(x)为Lagrange插值多项式。
(3) n = 1时称为线性插值: • , p1(x) = l0(x)f0 + l1(x)f1 • n = 2时称为抛物插值: • p2(x) = l0(x)f0 + l1(x)f1 + l2(x)f2
例2 已知离散数据如下: • (1) 求以x2=0,x3=1为节点的线性插值多项式,并预测x=0.3时f 的近似值。 • (2) 求以x1=-1,x2=0,x3=1为节点的二次插值多项式,并预测x=0.3时f 的近似值。 • 解:(1) 线性插值多项式是: • 由于p1(0.3)=0.3+1=1.3,所以当x=0.3时f 的近似值为1.3。
例2 已知离散数据如下: • (2) 求以x1=-1,x2=0,x3=1为节点的二次插值多项式,并预测x=0.3时f 的近似值。 • 解:(2) 二次插值多项式是: • 由于p2(0.3)=1.2475,所以当x=0.3时f 的近似值为1.2475。
Matlab代码: • syms x • x0 = [-1, 0, 1]; y0 = [0.5, 1, 2]; • n = length(x0); s = 0; • for k = 1:n • p = 1.0; • for j = 1:n • if j~=k • p = p*(x-x0(j))/(x0(k)-x0(j)); • end • end • s = s + p*y0(k); • end • s = vpa(s,8) • s = expand(s); s = simple(s) • s = subs(s,'0.3',x); s = vpa(s,8)
2. Newton插值 • L插值的缺点:每增加一个新节点,其插值基函数li(x)要重新计算,能否充分利用已有结果呢? • 为此作基函数: • 即将pn(x)表示为: • (5.1-10) • 这样,当增加一个新节点时,只需增加一个新项
利用插值条件(5.1-1)可得 • c0 = f0 • 称为一阶差商 • 称为二阶差商
… • 称为i阶差商 (i = 1,2,…,n) • 从而有 • pn(x) = f0 +f [x0,x1](x – x0) + … • +f [x0,x1,…,xn](x – x0)(x – x1)…(x – xn-1) (5.1-11) • 称为Newton插值多项式
注:差商的基本性质如下: • (1) i阶差商f [x0,x1,…,xi]可表为f (x0),f (x1),…,f (xi)的线性组合 • (2) 差商具有对称性,即差商与它所含节点的排列顺序无关 • (3) 若f为m次多项式,则i阶差商:
(4) 差商表 • 例3(p369) • 求(1) 次数不大于3的插值多项式p3(x)通过前四个数据点 • (2) 次数不大于4的插值多项式p4(x)通过所给五个数据点,并计算f (4.0)的近似值。
例3(p369) • 求(1) 次数不大于3的插值多项式p3(x)通过前四个数据点 • (2) 次数不大于4的插值多项式p4(x)通过所给五个数据点,并计算f (4.0)的近似值。 • 解:计算差商表如下:
解:计算差商表如下: • (1) • p3(x) = 14.2 + 2.118(x – 1.0) + 2.855(x – 1.0)(x – 2.7) • – 0.535(x – 1.0)(x – 2.7)(x – 3.2)
(1) p3(x) = 14.2 + 2.118(x – 1.0) + 2.855(x – 1.0)(x – 2.7) • – 0.535(x – 1.0)(x – 2.7)(x – 3.2) • (2) p4(x)=14.2+2.118(x–1.0)+2.855(x–1.0)(x–2.7) • –0.535(x–1.0)(x – 2.7)(x – 3.2) • + 0.266(x – 1.0)(x – 2.7)(x – 3.2)(x – 4.8) • = {[[0.266(x–4.8)–0.535](x–3.2)+2.855](x–2.7)+2.118}(x–1.0) +14.2 • f(4.0) • p4(4.0)={[[0.266(4.0–4.8)–0.535](4.0–3.2) + 2.855](4.0 – 2.7) + 2.118}(4.0 – 1.0) +14.2 • = 29.355
3. 插值余项 • 称R(x) = f(x) – pn(x)为插值余项 • 定理5.1-2 设f(x)在插值区间[a,b]上存在n + 1阶导数, • 则对于任意x [a,b],存在= (x) (a,b),使 • (5.1-12)
证明:显然R(xi) = 0,(i = 0,1,2,…,n), • 所以R(xi)至少有n + 1个零点 • 故可设 • 作辅助函数 • 则(t)在[a,b]上具有n + 1阶导数,且有n + 2个互异的零点:x0,x1,...,xn,x • 由罗尔定理,存在(a,b)使 • f (n+1)() – K(x)(n + 1)! = 0 •
注: • (1) 通常余项估计: • 其中 • (5.1-14) • (2) 对于指定x,应取靠近x的n + 1个节点作插值,以使 • | n + 1(x) |较小
例4(p371) • 设f(x) = ex,给出f (x)在[0,1]上n + 1个等距节点 • 处的函数值 • (1) 若按所给的函数值表用线性插值求 • f (x) = ex,(0 x 1) • 的近似值,使它们的绝对误差限不大于 ,问n应取多大 • (2) 每个表值f(xi)应取几位有效数字
解:(1) 由于x [0,1],故必有一个i,使得xi-1xxi, • 若p1(x)是f(x)在[xi-1,xi]上的线性插值多项式,且用p1(x)的值作为f (x) = ex的近似值,则由(5.1-15)及 • 有 • 由此解得n> 824。造表时一个合适的选取是n = 1000。此时步长为0.001。 • (2) 因0 x 1,有1 ex e,即个位数非零,由有小数位的概念知,每个表值至少应有7位有效数字。
4. 余项估计方法 • (5.1-12) • ——通常f (x)不知道,故f (n + 1)()也不知道 • 设[a,b]上n + 2个互异的零点:x0,x1,...,xn,xn+1。 • 取x0,x1,...,xn做n次插值多项式pn(1)(x) • 取x1,x2,...,xn+1做n次插值多项式pn(2)(x) • 由定理5.1-2有 • 其中1,2 (a,b)
若f (n + 1)(x)在(a,b)连续且变化较小,则有 • f (n + 1)(1) f (n + 1)(2) • 从而 • • (5.1-16) • 注:(1) 第一式为事后估计: • pn(2)(x)与f (x)之差可通过pn(2)(x) – pn(1)(x)来估计 • (2) 第二式表明: • 作f 的近似比pn(2)(x)作f 的近似更好。
5. 差商与导数的关系 • 将x看成与x0,x1,...,xn互异的节点,则按Newton插值公式得: • 上式令t = x,则因为f(x) = pn(x) (插值条件)
• 由余项公式(5.1-12)知 • (5.1-18)
6. 重节点差商公式 • 在(5.1-18)中令x0,x1,...,xn都趋于x,则也趋于x,从而有 • 注:要求f(n+1)(x)在(a,b)内变化很小,即要求f(x)的n+1阶差商变化很小。
7. Hermite插值 • 设f(xi) = fi,f '(xi) = f i',...,f (mi-1)(xi) = f i(mi-1), • i = 0,1,2,...,n, • 求次数尽可能低的多项式H(x),满足: • H (j)(xi) = f i(j)j = 0,1,...,mi – 1, i = 0,1,...,n (5.1-20) • 可以证明:存在唯一的满足(5.1-20)的多项式Hm(x), • ( ) • Hm(x)称为Hermite插值多项式,其余项: • (5.1-21)
8. Newton形式的Hermite插值 • 设f(xi) = fi,f '(xi) = f i',...,f (mi-1)(xi) = f i(mi-1), • i = 0,1,2,...,n • 显然,xi为f的mi重根,不妨设f 的零点:z0,z1,...,zm • 则 • 注:其中重节点差商按(5.1-19)计算。
例5 求次数最低的Hermite插值多项式H(x),使 • H(0) = - 1,H '(0) = - 2,H(1) = 0,H '(1) = 10,H ''(1) = 40,并求其余项表达式。 • 解:因为x = 0有一阶导数条件,所以为二重节点(余项的二重零点)x = 1有二阶导数条件,所以为三重节点
差商表如下: • 所以插值多项式为 • H(x) = – 1 – 2x + 3x2 + 6x2(x – 1) + 5x2(x – 1)2 • 余项为
9. 算法 • (1) 计算差商f [x0,x1],f [x0,x1,x2] … • (2) 计算多项式值(秦九韶法)
也可以使用二维数组存放:A[n][n] • 存储方式:x[i]存放节点,A[i, j]存放j阶差商, • i = 0,1,2,…,n,j = 0,1,2,…,m, • 其中零阶差商即函数值f(x),j阶差商有n – j + 1个
5.2 样条插值 • 1. 分段低次插值 • 为使插值更准确,插值节点之间间距应较小,即增加节点,此时插值多项式的次数增高——高次插值。 高次插值具有的震荡特性,使插值效果反而不好——龙格现象(p378例) 为了避免龙格现象,通常采用分段低次多项式插值,但此时在节点处一般不光滑。
2. 三次样条插值多项式 • 若S(x)在[a,b]上满足: • (1) S(x),S'(x),S''(x)连续; • (2) S(xi) = fi,i = 0,1,2,…,n ——插值条件 • (3) S在[xi,xi+1]上为三次多项式 • 则称S为[a,b]上的三次样条插值多项式
3. 样条插值公式 • 因为S在[xi,xi+1]上为三次多项式,所以其二阶导数必为一次式,记Mi表示S''(xi)(i = 0,1,2,…,n)则在[xi,xi+1]上有: • (5.2-1) • 其中hi = xi+1 – xi • 对(5.2-1)两次积分得: • 其中φ(x)为一次多项式,设为
由插值条件S(xi) = fi,S(xi+1) = fi+1得 • 从而在[xi,xi+1]上有 • i = 0,1,2,…,n – 1 (5.2-2) • 即S(x)有n + 1个未知的参数Mi (i = 0,1,2,…,n)需确定
而在(xi,xi+1)中有 • (5.2-3) • 由S ',S ''连续,在(5.2-3)中令 • xxi有 • xxi+1有 • 从而
因为S '(xi+0) = S '(xi – 0)所以 • (i = 0,1,2,…,n – 1) (5.2-4) • 记 , • 而 • (5.2-4)改写为:Mi – 1+ 2M i +M i + 1 =6f [x i – 1,x i,x i +1] • (i = 0,1,2,…,n – 1) (5.2-5) • 注:(5.2-5)有n – 1个方程,n + 1个未知数Mi。
注:(5.2-5)有n – 1个方程,n + 1个未知数Mi。 • 附加边界条件: • 第一类边界条件:S '(x0) = m0,S '(xn) = mn • 第二类边界条件:S ''(x0) = M0,S ''(xn) = Mn • 第三类边界条件:以b – a为周期 • S(x0 + 0) = S(xn – 0) • S '(x0 + 0) = S '(xn – 0) • S ''(x0 + 0) = S ''(xn – 0)
(1) 已知S ''(x0) = M0,S ''(xn) = Mn • 则(5.2-5)改写为: • (5.2-6)
(2) 已知S '(x0) = m0,S '(xn) = mn • 则因为 • 所以 • (5.2-8)
(3) 已知S以b – a为周期,即f 0 = fn,m0 = mn,M0 = Mn • (5.2-9) • 其中,x-1 = xn-1,