第四章插值和曲線擬合

1、第四章插值和曲線擬合第四章插值和曲線擬合在實際問題和科學實驗中所遇到的函數y=f(x),往往沒有解析表達式 , 只能根據試驗觀察或其它方法提供一系列點的函數值; 有時盡管可以寫出表達式，但是比較復雜, 直接使用它感到不方便。我們經常需要利用已知的數據去尋求某個簡單的函數(x)來逼近f(x),即用(x)作為f(x)的近似表達式。本章的插值法和曲線擬合就是兩種用來求求 f(x) 的近似函數的近似函數(x) 的重要方法。第一節(jié)插值法的基本理論第一節(jié)插值法的基本理論一、 插值問題 設函數 y = f(x) 給出了一組函數值 yi = f(xi) , i = 0, 1, , n ，或者給出了如下的一張表

2、 x0 , x1 , x2 , , xn y0 , y1 , y2 , , yn構造一個簡單的函數(x) 作為f(x)的近似表達式，以滿足 (xi) = yi , i = 0 , 1 , , n我們稱這樣的問題為插值問題插值問題。 其中(xi) = yi 稱為插值原則插值原則;(x)稱為f(x)的插值函數；f(x)稱為被插值函數; x0 , x1 , x2 , , xn稱為插值基點(或節(jié)點)。 根據插值原則插值原則求其余點x的函數值(x)(x)稱為插值,x稱為插值點；根據插值原則插值原則求f(x)近似函數(x)(x)的方法稱為插值法。插值法的幾何意義插值法的幾何意義 插值法的幾何意義就是通過n

3、+1個點: (xi,yi) (i=0,1,2,n) 作一條近似曲線y= (x) 代替y=f(x)。如下圖所示。xxnx2x1x0Xn-1y0(x0,y0)(x1,y1)(x2,y2)(xn-1,yn-1)(xn,yn)y= (x)y=f(x)插值函數插值函數(x)的類型的類型 在插值問題中，插值函數(x)的類型可有不同的選擇,如代代數多項式、三角數多項式、三角多項式、有理函數多項式、有理函數等，但是最簡單而常用的是代數多項式，這時就稱為代數多項式代數多項式插值插值。在本章，我們主要討論代數多項式插值。 代數多項式插值的任務就是根據 n+1個點 x0 , x1 , x2 , , xn y0 ,

4、y1 , y2 , , yn構造一個次數不超過 n 的多項式 Pn(x) = a0 + a1x + a2x2 + + anxn使?jié)M足插值原則插值原則 Pn(xi) = yi , i = 0 , 1 , , n 。 Pn(x)稱為稱為 f(x) 的的 n次插值次插值多項式多項式。 二、插值多項式的誤差二、插值多項式的誤差 函數 f(x)用n次插值多項式Pn(x)近似代替時，截斷誤差記為Rn(x)=f(x)-Pn(x)稱 Rn(x)為n次插值多項式Pn(x)的余項。 定理 設函數f(x)在包含基點x0 , x1 , x2 , , xn 的區(qū)間a,b上具有n+1階導數, Pn(x)為滿足Pn(xi)

5、 = yi的n次插值多項式,則對任一點xa,b,總存在相應的點 ,使 其中 wn+1(x) = (x-x0)(x-x1) (x-xn),(ba)()!1()()(1)1(xwnfxRnnn第二節(jié)第二節(jié) 拉格朗日插值拉格朗日插值 為了得到n次的拉格朗日插值多項式，我們從最簡單的一次、二次插值開始。 一、一次插值（線性插值）一、一次插值（線性插值） 已知 x0 x1 求 P1(x) y0 y1 因 P1(x0)=y0 , P1(x1)=y1 所以 P1(x)=y0+(y1-y0)/(x1-x0)*(x-x0) （線性插值多項式） 上式可改寫為： P1(x)=y0L0(x)+y1L1(x ) （拉格

6、朗日線性插值多項式） L0(x)=(x-x1)/(x0-x1),L1(x)=(x-x0)/(x1-x0) L0(x)、L1(x)特點： L0(x)= 1 , x = x0 L1(x)= 1 , x = x1 0 , x = x1 , 0 , x = x0線性插值舉例線性插值舉例例 已知 1001/2 =10，1211/2 =11 求 1151/2解 P1(x) = y0+(y1-y0)/(x1-x0)*(x-x0) P1(115) = 10+(11-10)/(121-100)*(115-100) 或 P1(x)=y0L0(x)+y1L1(x ) P1(115) = 10*(115-121)/(

7、100-121) +11*(115-100)/(121-100) 二、二次插值二、二次插值(拋物線插值拋物線插值) 二次插值問題：已知f(x)在三個互異點x0,x1,x2的函數值y0,y1,y2 ,要構造次數不超過二次的多項式P2(x)=a0+a1x+a2x2,使?jié)M足 P2(xi)=yi , i = 0, 1, 2 設 P2(x) = L0(x)y0+L1(x)y1+L2(x)y2 , 則當x = x0 時, P2(x0) = y0 L0(x) = 1, L1(x) = 0, L2(x) = 0 當x = x1時，P2(x1) = y1 L0(x) = 0, L1(x) = 1, L2(x)

8、= 0當x = x2時，P2(x2) = y2 L0(x) = 0, L1(x) = 0, L2(x) = 1由上知 L0(x) = 1, x = x0 0, x = x1, x2令 L0(x)=A0(x-x1)(x-x2) 則 A0=1/(x0-x1)(x0-x2)所以 L0(x)= (x-x1)(x-x2)/(x0-x1)(x0-x2)同理可得 L1(x)=(x-x0)(x-x2)/(x1-x0)(x1-x2) L2(x)=(x-x0)(x-x1)/(x2-x0)(x2-x1)綜上可得 P2(x)=y0(x-x1)(x-x2)/(x0-x1)(x0-x2) +y1(x-x0)(x-x2)/

9、(x1-x0)(x1-x2) +y2 (x-x0)(x-x1)/(x2-x0)(x2-x1)該式稱為拉格朗日二次插值多項式。二次插值舉例二次插值舉例 例 已知函數y=f(x)的觀測數據如下表所示，試求其拉格朗日插值多項式，并計算f(1.5)的近似值。 4 -1 2 y 2 1 0 x 解 P2(x) = y0(x-x1)(x-x2)/(x0-x1)(x0-x2) +y1(x-x0)(x-x2) /(x1-x0)(x1-x2)+y2 (x-x0)(x-x1)/(x2-x0)(x2-x1) = 2*(x-1)(x-2)/(0-1)(0-2) +(-1)*(x-0)(x-2)/(1-0)(1- 2)

10、+4*(x-0)(x-1)/(2-0)(2-1) = 4x2-7x+2 f(1.5) P2 (1.5)=4*1.52-7*1.5+2 = 0.5 三、三、n次拉格朗日插值次拉格朗日插值 仿照P2 (x)的構造方法，可得出 Pn(x)=L0(x)y0+L1(x)y1+Ln(x)yn其中 L0(x)=(x-x1)(x-x2)(x-xn)/ (x0-x1)(x0-x2)(x0-xn) Lk(x)= (x-x0)(x-xk-1)(x-xk+1) (x-xn) /(xk-x0)(xk-xk-1)(xk-xk+1) (xk-xn) ( k = 0, 1, , n ) 這就是n次拉格朗日插值多項式。 也可寫

11、為 niinikkkikniiinyxxxxyxLxP0, 00)()(n次拉格朗日插值舉例次拉格朗日插值舉例例例 已知函數表 x 1.1275 1.1503 1.1735 1.972 y 0.1191 0.13954 0.15932 0.17903應用朗格拉日插值公式計算f(1.1300)的近似值。 解 P3(x) = L0(x)y0+L1(x)y1+L2(x)y2 +L3(x)y3 = f(1.1300) P3(1.1300) = 0.1214 n次拉格朗日插值計算機實現次拉格朗日插值計算機實現 按按n次拉格朗日插值公式實現次拉格朗日插值公式實現分段插值 七、八次以上的高次插值在實際中很少

12、采用。因為理論研究和實例都表明，插值基點增加并不能保證 Pn(x)在非基點處逼近f(x)的精度得到提高, 某些情況下甚至誤差反而變大。所以總是對每個插值點 x選擇其附附近的幾個插值基點近的幾個插值基點作低次內插(將 x 放在插值基點之間),或者采用分段低次插值(一次、二次插值 )。 為什么要選擇 x 附近的幾個插值基點？ 根據)()!1()()(1)1(xwnfxRnnn其中 wn+1(x) = (x-x0)(x-x1) (x-xn)第三節(jié)第三節(jié) 牛頓插值牛頓插值 拉格朗日插值多項式結構對稱, 使用方便, 但公式不具備遞推性,當需要增加基點時必須全部重新計算。因此，我們希望構造具有如下形式的插

13、值多項式 Pn(x) = a0 + a1(x-x0) + a2(x-x0)(x-x1) + + an(x-x0)(x-x1) (x-xn-1)這種形式的優(yōu)點是便于改變基點數，每增加一個基點只需增加相應的一項即可 (具有遞推性) 。為了確定出a0 、a1、an , 我們就需要討論牛頓差商插值多項式。下面首先介紹差商的概念。一、差商及差商表一、差商及差商表 1. 差商定義差商定義 在區(qū)間a,b上，函數f(x)關于兩點xi , xj的一階差商定義為 f xi, xj = f(xj)-f(xi)/(xj-xi)f(x)關于三點xi, xj, xk的二階差商定義為 f xi, xj, xk=(fxj,

14、xk - fxi, xj)/(xk-xi)f(x)關于k+1個點xi-k, xi-k+1, , xi 的k階差商定義為 f xi-k,xi-k+1, xi = (f xi-k+1,xi f xi-k, , xi-1)/(xi-xi-k) f(x)關于一個點 xi 的零階差商定義為函數本身，即 f xi = f(xi) 不論幾階差商，差商均有對稱性（任意改變基點的次序后其值不變）。即 f x0,x1, xk = f 其中 是 x0,x1, xk 的任一種排列。（證略）kjjjxxx,.,10kjjjxxx,.,10 2. 差商表差商表對于給定的基點及其函數值，我們可按表計算各階差商，這樣的表就叫

15、差商表。如下:xi四階差商一階差商二階差商三階差商x0 x4x3 x2x1f(x4)f(x3)f(x2)f(x1)f(x0)f x3, x4f x2, x3f x1, x2f x0, x1f x1, x2, x3f x2, x3, x4f x0, x1, x2f x0, x1, x2, x3f x1, x2, x3, x4 f x0 , x1, x2 , x3 , x4.零階差商二、牛頓差商插值多項式二、牛頓差商插值多項式 由差商定義和差商性質有f(x) = f(x0)+f x0,x(x-x0) ( fx0,x=f(x)-f(x0)/(x-x0) )f x0, x= f x0,x1+f x0,

16、x1,x(x-x1)f x0, x1,x = f x0,x1,x2+ f x0,x1,x2,x(x-x2) f x0,x1,xn-1,x = f x0,x1,xn+ fx0,x1,xn,x(x-xn)f(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) + f x0,x1,xn,x(x-x0)(x-x1)(x-xn) = Pn(x) + Rn(x)Pn(x) = f(x0) + f x0,x1(x-x0) + f x0,x1,x2(x-x0)(x-x1) + + f x0

17、,x1,xn(x-x0)(x-x1)(x-xn-1) Pn(x)由于滿足 Pn(xi)=f(xi) 稱作 n 次牛頓(差商)插值多項式。 Rn(x) = f(x)-Pn(x) = w(x)*f(n+1)()/(n+1)! （w(x)=(x-x0)(x-x1)(x-xn)）稱為n次牛頓插值多項式的余項。牛頓差商插值多項式的兩個特殊形式牛頓差商插值多項式的兩個特殊形式 牛頓(差商)插值多項式為 當n=1時 P1(x) = f(x0) + f x0,x1(x-x0)即為線性插值。 當n=2時 P2(x)=f(x0)+f x0,x1(x-x0)+f x0,x1,x2(x-x0)(x-x1)即為拋物線插

18、值。 可見, 增加一個基點, 只是增加了f x0,x1,x2(x-x0)(x-x1)這一項。 注意注意, 可以證明牛頓插值多項式與拉格朗日插值多項式是等價可以證明牛頓插值多項式與拉格朗日插值多項式是等價的的, 只不過形式不一樣而已。所以只不過形式不一樣而已。所以, 兩者的截斷誤差是一樣的。兩者的截斷誤差是一樣的。 Pn(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)牛頓差商插值多項式舉例牛頓差商插值多項式舉例 例 已知函數y=f(x)的觀測數據如下表所示，試用全部基點構

19、造牛頓差商插值多項式，并用二次插值求f(3)的近似值。 解 差商表為P4(x) = 1+2(x-0)+0(x-0)(x-2)+(-1)(x-0)(x-2)(x-4)+(x-0)(x-2)(x-4)(x-5) = x4-12x3+44x2-46x+1 13 -4 9 5 1 f(x) 6 5 4 2 0 x xi 零階 一階 二階 三階 四階 0 1 2 5 2 4 9 2 0 5 -4 -13 -5 -1 6 13 17 15 5 1牛頓差商插值多項式舉例（續(xù)）牛頓差商插值多項式舉例（續(xù)） 用二次插值求 f(3) 時, 取 x0 = 2, x1 = 4, x2 = 5 f(3) P2(3) =

20、 f(2)+f2,4(3-2)+f2,4,5(3-2)(3-4) = 5+2(3-2)-5(3-2)(3-4) = 5+2+5 = 12 xi 零階 一階 二階 三階 四階 0 1 2 5 2 4 9 2 0 5 -4 -13 -5 -1 6 13 17 15 5 1 差商表為P2(x)=f(x0)+f x0,x1(x-x0)+f x0,x1,x2(x-x0)(x-x1)牛頓差商插值多項式的計算機實現牛頓差商插值多項式的計算機實現 按牛頓差商插值多項式公式實現。分兩大步： 1. 求各階差商 2. 按秦九韶方法求多項式的值。第四節(jié)第四節(jié) 曲線擬合曲線擬合 插值法和曲線擬合都是用來求列表函數f(x

21、)(只知道一些點的函數稱為列表函數)的近似函數(x) 。插值法插值法求出的近似曲線求出的近似曲線 y =(x) 要完全通過所有要完全通過所有 n+1 個已知個已知點點 ( 即要滿足插值原則即要滿足插值原則 ) ; 而曲線擬合求出的近似曲線而曲線擬合求出的近似曲線 y =(x) 不要求完全通過所有不要求完全通過所有 n+1個已知點個已知點 , 只要求求只要求求得的近似曲線得的近似曲線 y =(x) 能反映數據的基本趨勢即可。能反映數據的基本趨勢即可。曲線擬合求得的近似曲線 y =(x) 比插值法求得的近似曲線 y = (x)更能反映客觀實際。因為列表函數中的點往往是通過實驗、測量或計算得來的,

22、而實驗、測量或計算得來的數據經常帶有誤差, 如果要求所得出的曲線 y=(x)通過所有n+1個已知點 (xi,yi), 就會使曲線 y =(x) 保留著這些誤差，而這是我們所不希望的。一、曲線擬合問題一、曲線擬合問題 設函數設函數 y = f(x) 在在 n+1 個互異點的觀測數據為個互異點的觀測數據為 x0 , x1 , , xn y0 , y1 , , yn 構造函數構造函數(x)在包含全部基點的區(qū)間上在包含全部基點的區(qū)間上“最好最好”地逼近地逼近(或靠或靠近近)f(x), 這就是曲線擬合問題。這就是曲線擬合問題。如下圖所示如下圖所示, ,就是使曲線就是使曲線y=(x)盡量靠近已知點盡量靠近

23、已知點(xi,yi) (i=0,1,2,n)。 xy0y=(x)(x0,y0)(x1,y1)(xn,yn)(xn-1,yn-1)(x2,y2)nn-1210二、二、 最小二乘法最小二乘法 假設 y =(x)其中(x)=a0+a1x+a2x2+amxm）為給定的一組數據(xi,yi)(i=0,1,2,n)的擬合曲線，則將這n+1個點代入(x)得以下式子 a0+a1x0+a2x02+amx0m y0a0+a1x1+a2x12+amx1m y1 a0+a1xn+a2xn2+amxnm yn最小二乘法最小二乘法(續(xù)續(xù))若將a0+a1x0+a2x02+amx0m y0a0+a1x1+a2x12+amx1

24、m y1 a0+a1xn+a2xn2+amxnm = yn 中的“”換為“=”該式變?yōu)榉匠探Ma0+a1x0+a2x02+amx0m = y0 a0+a1xn+a2xn2+amxnm yn a0+a1x1+a2x12+amx1m = y1 方程組中有n+1個方程、m+1個未知數。 若 n+1=m+1, 方程組有唯一解; 若n+1m+1, 方程組無解(無精確解), 稱為超定方程組超定方程組。最小二乘法最小二乘法(續(xù)續(xù)) 超定方程組無精確解但可求其近似解。 那么解近似到什么程度才算近似解？這有不同的標準。 若所求得的近似解使得誤差平方和( )達到最小，我們稱這組近似解為最優(yōu)近似解最優(yōu)近似解。根據誤差

25、平方根據誤差平方和達到最小這一標準求最優(yōu)近似解的方法就稱為最小二乘和達到最小這一標準求最優(yōu)近似解的方法就稱為最小二乘法法。 下面具體解釋一下什么是超定方程組的最小二乘法。 a0+a1x1+a2x12+amx1m = y1 a0+a1xn+a2xn2+amxnm = yn a0+a1x0+a2x02+amx0m = y0 njmijijiyxa020)(超定方程組的最小二乘法超定方程組的最小二乘法 設超定方程組 根據高數知識， 達到最小必須滿足條件 a0+a1x1+a2x12+amx1m = y1a0+a1x0+a2x02+amx0m = y0 a0+a1xn+a2xn2+amxnm = yn

26、( i = 0, 1, 2, , m )，誤差平方和為),.,()(10020mnjmijijiaaaQyxa ),.,1 ,0(0miaQi),.,(10maaaQ),.,1 ,0(0miaQi方程組近似解的方法稱為解超定方程組的最小二乘法最小二乘法。 (即y=(x)=a0+a1x+a2x2+amxm 在xi處的值與yi的誤差平方和)解此方程組便可求出最優(yōu)近似解。稱為法方程組。根據法方程組求超定的近似解為ai根據超定方程組的最小二乘法知, 上式系數 ai (i=0,1,2,m)可由法方程組即即求得。niiminimimniminiminiminiiinimimniiniiniiniinimimniiniiyxxaxaxaxayxxaxaxaxayxaxaxaan002022011000010320210000022010.) 1( 三、代數多項式擬合三、代數多項式擬合 若擬合函數形式(x)=a0+