插值與數(shù)據(jù)擬合

1、插值與數(shù)據(jù)擬合第1頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 插值與數(shù)據(jù)擬合就是通過一些已知數(shù)據(jù)去確定某類函數(shù)的參數(shù)或?qū)ふ夷硞€(gè)近似函數(shù)，使所得的函數(shù)與已知數(shù)據(jù)具有較高的精度，并且能夠使用數(shù)學(xué)分析的工具分析數(shù)據(jù)所反映的對(duì)象的性質(zhì) 幾種常用的方法： 1、一般插值法 2、樣條插值法 3、最小二乘曲線 4、曲面的擬合數(shù)據(jù)擬合與插值建模第2頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 上大學(xué)二年級(jí)的小華正在做作業(yè)，“爸爸，計(jì)算這道題要用到sin ，可是我的計(jì)算器壞了，怎么辦。”當(dāng)工程師的老張從厚厚的一摞舊書底下抽出一本數(shù)學(xué)用表來，“給你，這是我念大學(xué)時(shí)用的，那時(shí)候啊，計(jì)算器

2、聽都沒聽說過?！毙∪A拿著表翻了一會(huì)兒，無奈的說：“表上每 才有一個(gè)函數(shù)值，這里只sin 和sin “表中沒有的，都可以用插值方法計(jì)算”“插值！我們的數(shù)學(xué)實(shí)驗(yàn)課就要學(xué)了，不過，今天我要先自己想個(gè)辦法，用這個(gè)算出sin ”這本四位數(shù)學(xué)用表給出sin 0.576，sin =0.5783。小華認(rèn)為在sin 到sin 這樣小的范圍內(nèi)，正弦可以近似為線性函數(shù)，于是很容易地得到Sin =0.576+(0.5783-0.5760)0.6=0.5774第3頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 聰明的小華用的這個(gè)辦法是一種插值方法分段線性插值。實(shí)際上，插值可以理解為，要根據(jù)一個(gè)用表格表示的函

3、數(shù)，計(jì)算表中沒有的函數(shù)值，表中有的。表中有的，如（sin ，0.5760）(sin ,0.5783)稱為節(jié)點(diǎn)；要計(jì)算的，如sin ，稱為插值點(diǎn)，結(jié)果（0.5774）即為插值。小華作的線性函數(shù)為插值函數(shù)，插值函數(shù)所表示的直線當(dāng)然要通過節(jié)點(diǎn)。第4頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 插值最初來源于天體計(jì)算由若干觀測值（即節(jié)點(diǎn)）計(jì)算任意時(shí)刻星球的位置（即插值點(diǎn)和插值）的需要?，F(xiàn)在，雖然人們已很少需要用它從函數(shù)表計(jì)算函數(shù)值了，但是插值仍然在諸如機(jī)械加工等工程技術(shù)和數(shù)據(jù)處理等科學(xué)研究中有著許多直接的應(yīng)用，另一方面，插值又是數(shù)值微分、數(shù)值積分、常微分方程數(shù)值等數(shù)值計(jì)算的基礎(chǔ)。第5頁，

4、共78頁，2022年，5月20日，15點(diǎn)20分，星期三 幾天后，小華在物理實(shí)驗(yàn)里又碰到一個(gè)看起來非常類似的問題：有一只對(duì)溫度敏感的電阻，已經(jīng)測得了一組溫度T和電阻R數(shù)據(jù)。 現(xiàn)在想知道 時(shí)的電阻多大。溫度t(0C) 20.5 32.7 51.0 73.0 95.7電阻R() 765 826 873 942 1032第6頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 小華征求老師的意見，老師給了他兩點(diǎn)提示：一是在直角坐標(biāo)系中把5個(gè)點(diǎn)（T，R）畫一下，看看電阻R和溫度T之間大致有什么關(guān)系；二是測量數(shù)據(jù)總有相當(dāng)大的誤差，這與用函數(shù)表作插值計(jì)算應(yīng)該有不同之處吧（雖然函數(shù)表也存在舍入誤差，但很

5、小，可以認(rèn)為表中數(shù)值是精確的） 通過圖形小華看到，R與T大致呈直線關(guān)系，于是用手畫了一條靠近5個(gè)點(diǎn)的直線，又想起中學(xué)物理學(xué)過，金屬材料的電阻率與溫度成正比，從而確定R與T的關(guān)系應(yīng)該是 R=at+b 其中a，b為待定常數(shù)。第7頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 正是由于測量誤差的存在，由R= at+b表示的直線不可能通過全部5個(gè)點(diǎn)，所以，與插值曲線要通過全部節(jié)點(diǎn)不同，小華打算作一條盡量靠近所有的點(diǎn)的直線，求出a，b待定常數(shù)，由此計(jì)算t= 的R就十分簡單了。第8頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 根據(jù)一組（二組）數(shù)據(jù)，即平面上的若干點(diǎn)，確定一個(gè)一元函

6、數(shù)，即曲線，使這些節(jié)點(diǎn)與曲線總體來說盡量接近，這就是曲線擬合。 函數(shù)值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函 數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法是完全不同的。第9頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 數(shù) 據(jù) 擬 合1. 擬合的基本原理；2. 最小二乘擬合；3. 用Matlab作最小二乘擬合；4.如何用擬合解決實(shí)際問題。第10頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 t (h) 0.25 0.5 1 1.5 2 3 4 6 8c (g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01 對(duì)

7、某人用快速靜脈注射方式一次性注射某種藥物300mg后，經(jīng)過時(shí)間t采集血樣，測得血藥濃度c如下表：求血藥濃度隨時(shí)間的變化規(guī)律c(t).半對(duì)數(shù)坐標(biāo)系(semilogy)下的圖形Log10c(t)=a t + b引例1：血藥濃度的變化規(guī)律第11頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三曲 線 擬 合 問 題 的 提 法 已知一組（二維）數(shù)據(jù)，即平面上 n個(gè)點(diǎn)（xi,yi) i=1,n, 尋求一個(gè)函數(shù)（曲線）y=f(x), 使 f(x) 在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。 +xyy=f(x)(xi,yi)ii 為點(diǎn)（xi,yi) 與曲線 y=f(x) 的距離第12頁，

8、共78頁，2022年，5月20日，15點(diǎn)20分，星期三最小二乘擬合 第一步：先選定一類函數(shù)f(x，a1，a2， ，am) 其準(zhǔn)則為（最小二乘準(zhǔn)則）：使n個(gè)點(diǎn)（xi,yi) 與曲線 y=f(x ，a1，a2， ，am) 的距離i 的平方和最小 。其中 a1,a2, am 為待定常數(shù)。f可以為一些簡單的“基函數(shù)”（如冪函數(shù)，三角函數(shù)等等）的線性組合：第二步：確定參數(shù)a1,a2, am，第13頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三記問題歸結(jié)為，求 a1,a2, am 使 J(a1,a2, am) 最小。 這樣的擬合稱為最小二乘擬合。 除了最小二乘準(zhǔn)則（即各點(diǎn)誤差的平方和最小），你

9、認(rèn)為還可以用怎樣的擬合準(zhǔn)則？ 比較起來，最小二乘準(zhǔn)則有什么優(yōu)點(diǎn)？思考第14頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三最小二乘擬合函數(shù) f(x，a1， am)的選取 +f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1exp(a2x)+f=a1exp(a2x)1. 通過機(jī)理分析建立數(shù)學(xué)模型來確定 f；2. 將數(shù)據(jù) (xi,yi) i=1, n 作圖，通過直觀判斷確定 f：第15頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三2. 作一般的最小二乘曲線擬合，可利用已有程序curvefit,其調(diào)用格式為： a=curvefit(f, a0, x,

10、y) 1. 作多項(xiàng)式f(x)=a1xm+ +amx+am+1函數(shù)擬合,可利用已有程序polyfit,其調(diào)用格式為:a=polyfit(x,y,m)用MATLAB作最小二乘擬合數(shù)據(jù)點(diǎn)擬合多項(xiàng)式次數(shù)系數(shù)注：f為擬合函數(shù)y=f(a,x)的函數(shù)M文件，f(a,x)為擬合函數(shù)。數(shù)據(jù)點(diǎn)待定常數(shù)a的初值函數(shù)M文件第16頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三用MATLAB作多項(xiàng)式最小二乘擬合2.用命令polyfit(x,y,m)得到 a1=3.3940, a2=702.49181. 選取函數(shù) R= a1t+a2溫度t(0C) 20.5 32.7 51.0 73.0 95.7電阻R() 76

11、5 826 873 942 1032例. 由數(shù)據(jù)擬合R=f(t)第17頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三用MATLAB作最小二乘曲線擬合例：用函數(shù)f(x)=a1*exp(-a2*x)+a3*exp(-a4*x)擬合下列數(shù)據(jù)點(diǎn)：xdata=0:.1:2 ydata=5.8955 3.5639 2.5173 1.9790 1.8990 1.3938 1.1359 1.0096 1.0343 0.8435 0.6856 0.6100 0.5392 0.3946 0.3903 0.5474 0.3459 0.1370 0.2211 0.1704 0.2636用命令curvefi

12、t(f, a0, x, y) 第18頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三擬合的應(yīng)用參數(shù)辨識(shí) 數(shù)學(xué)建模的方法：機(jī)理分析和測試分析。 機(jī)理分析是根據(jù)對(duì)客觀事物特性的認(rèn)識(shí)，找出反映內(nèi)部機(jī)理的數(shù)量規(guī)律，建立的模型常有明確的物理意義。 測試分析將研究的對(duì)象看作一個(gè)“黑箱”，通過對(duì)實(shí)驗(yàn)數(shù)據(jù)的統(tǒng)計(jì)分析，找出與數(shù)據(jù)擬合得最好的模型。 機(jī)理分析模型結(jié)構(gòu)實(shí)驗(yàn)數(shù)據(jù)未知參數(shù)第19頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三范例：薄膜滲透率的測定 一、問題： 某種醫(yī)用薄膜，具有從高濃度的溶液向低濃度的溶液擴(kuò)散的功能，在試制時(shí)需測定薄膜被物質(zhì)分子穿透的能力。 測定方法：用面積為S的薄膜

13、將容器分成體積分別為 的兩部份，在兩部分中分別注滿該物質(zhì)的兩種不同濃度的溶液。此時(shí)該物質(zhì)分子就會(huì)從高濃度溶液穿過薄膜向低濃度溶液中擴(kuò)散。平均每單位時(shí)間通過單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比，比例系數(shù)K表征了薄膜被該物質(zhì)分子穿透的能力，稱為滲透率。定時(shí)測量容器中薄膜某一側(cè)的溶液濃度，以此確定K。VAVBS第20頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三二、問題分析 考察時(shí)段t，t+t薄膜兩側(cè)容器中該物質(zhì)質(zhì)量的變化。 設(shè) ，對(duì)容器的B部分溶液濃度的測試結(jié)果如下表：（濃度單位 ） 1） 在容器的一側(cè)，物質(zhì)質(zhì)量的增加是由于另一側(cè)的物質(zhì)向該側(cè)滲透的結(jié)果，因此物質(zhì)質(zhì)量的增量

14、應(yīng)等于另一側(cè)的該物質(zhì)向這側(cè)的滲透量。第21頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 以容器A側(cè)為例，在時(shí)段t，t+t物質(zhì)質(zhì)量的增量為：分別表示在時(shí)刻t膜兩側(cè)溶液設(shè)的濃度，濃度單位: 由于平均每單位時(shí)間通過單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比，比例系數(shù)為K。 因此，在時(shí)段t，t+t，從B側(cè)滲透至A側(cè)的該物質(zhì)的質(zhì)量為：第22頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三于是有：兩邊除以t，并令t0取極限再稍加整理即得：分別表示在初始時(shí)刻兩側(cè)溶液的濃度其中（1）2) 注意到整個(gè)容器的溶液中含有該物質(zhì)的質(zhì)量不變,與初始時(shí)刻該物質(zhì)的含量相同，因此 第23頁，共

15、78頁，2022年，5月20日，15點(diǎn)20分，星期三從而：加上初值條件：代入式（1）得：便可得出CB(t)的變化規(guī)律，從而根據(jù)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行擬合，估計(jì)出參數(shù)K, 。第24頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三三、數(shù)學(xué)模型假設(shè)：1）薄膜兩側(cè)的溶液始終是均勻的；2）平均每單位時(shí)間通過單位面積薄膜的物質(zhì)分子量與膜兩側(cè)溶液的濃度差成正比。3）薄膜是雙向同性的即物質(zhì)從膜的任何一側(cè)向另一側(cè)滲透的性能是相同的?；诩僭O(shè)和前面的分析，B側(cè)的濃度CB(t)應(yīng)滿足如下微分方程和初始條件：第25頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三四、求解方法：1. 函數(shù)擬合法前面得到的模型是

16、一個(gè)帶初值的一階線性微分方程，解之得：問題歸結(jié)為利用CB在時(shí)刻tj的測量數(shù)據(jù)Cj(j=1,2,.,N)來辨識(shí) K 和 。第26頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三引入從而 用函數(shù)CB(t)來擬合所給的實(shí)驗(yàn)數(shù)據(jù)，從而估計(jì)出其中的參數(shù)a,b,K。將代入上式有：第27頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三用MATLAB軟件進(jìn)行計(jì)算.1）編寫函數(shù)M-文件 nongdu.mfunction f = nongdu(x,tdata)f = x(1)+x(2)*exp(-0.02*x(3)*tdata);其中 x(1) = a;x(2) = b;x(3) = k;2)

17、 在工作空間中執(zhí)行以下命令(test1.m) tdata = linspace(100,1000,10); cdata =4.54 4.99 5.35 5.65 5.90 6.10 . 6.26 6.39 6.50 6.59; x0 = 0.2,0.05,0.05; x = curvefit(nongdu,x0,tdata,cdata)3) 輸出結(jié)果: x = 0.007 -0.003 0.1012 即 k = 0.1012, a = 0.007, b = -0.003, 第28頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三進(jìn)一步求得：2. 非線性規(guī)劃法 利用CB在時(shí)刻tj的測量數(shù)

18、據(jù)Cj(j=1,2,.,N)來辨識(shí) K 和 。問題可轉(zhuǎn)化為求函數(shù)即求函數(shù)的最小值點(diǎn)（K，a，b）。第29頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三3. 導(dǎo)函數(shù)擬合法前面得到的微分方程為：令上式變?yōu)椋哼@可以看作隨CB的變化規(guī)律(j=1, 2,.,N)若知道一組數(shù)據(jù)則可用最小二乘擬合的方法來求出函數(shù)中的未知參數(shù)K和h。第30頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三即為求參數(shù)K, a使下列誤差函數(shù)達(dá)到最小： 該問題等價(jià)于用函 數(shù) f(K,a,CBB)來擬合數(shù)據(jù)(j=1, 2,.,N)第31頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三用MATLAB軟件進(jìn)

19、行計(jì)算.求數(shù)據(jù)點(diǎn)(j=1, 2,.,N)tdata = linspace(100,1000,10);cdata = 1e-05.*454 499 535 565 590. 610 626 639 650 659;d,ifail=e01bef(tdata,cdata);cj,dcj=e01bgf(tdata,cdata,d,tdata);1）編寫函數(shù)M-文件 baomof.mfunction f=baomof(x,cdata)f=x(1)*(0.01*x(2)-0.02*cdata)其中 x(1) = K; x(2) =h2) 編寫命令M文件(baomo21.m)第32頁，共78頁，2022年，

20、5月20日，15點(diǎn)20分，星期三3) 輸出結(jié)果: x = 0.1009 0.014 即 k = 0.1009, h = 0.014作函數(shù)擬合x0=0.2,0.1;x=curvefit(baomof,x0,cdata,dcj)第33頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三4.線性化迭代法前面帶初始條件的一階線性微分方程的解為其中： 如果得到了參數(shù)K的一個(gè)較好的近似值K*，則將CB(t)關(guān)于K在K*處展開，略去K的二次及以上的項(xiàng)得CB(t)的一個(gè)近似式通過極小化第34頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三確定a, b, d, 再由K=d/0.02b得到K*的修

21、正值K。K*K*- K, 得到K的一個(gè)新的近似值，用同樣的方法再求新的修正值K。這個(gè)過程可以不斷重復(fù)，直到修正值足夠小為止。1）當(dāng)K的初值取為k=0.3時(shí)，出現(xiàn)奇異情況，迭代不收斂；2）當(dāng)K的初值取為k=0.2時(shí)，經(jīng)四次迭代，已經(jīng)收斂到一個(gè)很好的解。迭代結(jié)果如下表。第35頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三五、結(jié)果及誤差分析 幾種方法得出的結(jié)果及相應(yīng)的誤差總結(jié)于下表，誤差為計(jì)算數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)之差的平方和。注：導(dǎo)函數(shù)擬合法得出的參數(shù)值精度有限，線性化迭代法要求參數(shù)的初值比較接近精確值。因此可將導(dǎo)函數(shù)擬合法和線性化迭代法結(jié)合起來使用，把前者得到的參數(shù)K的值作為迭代法中K的初值

22、，這樣可使迭代法收斂或收斂更快。3）取K的初值為k=0.1009,只一次迭代就得到2）中的最后結(jié)果。第36頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三第37頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三函數(shù)擬合法的擬合效果第38頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三求解參數(shù)辨識(shí)模型的方法： 函數(shù)擬合； 非線性規(guī)劃； 導(dǎo)函數(shù)擬合； 線性化迭代； 其它方法 。 第39頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三用Logistic模擬水稻葉伸長生長時(shí)間11.82.63.44.14.85.46.16.87.48.1重量0.30.50.91.4

23、2.53.24.37.610.114.418.5時(shí)間8.89.410.110.811.712.413.114.415.115.7重量23.025.230.433.738.841.743.744.845.545.3生長觀測記錄數(shù)據(jù)第40頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三模型表達(dá)式：第41頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三程序！第42頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三關(guān)于polyfit命令命令：p=polyfit(x,y,n)（1）x與y為模擬數(shù)據(jù)（2）n為擬合多項(xiàng)式的次數(shù)（3）當(dāng)n=1時(shí)為用最小二乘法進(jìn)行 直線擬合（4）得到

24、的向量p為長度n+1向量， 對(duì)應(yīng)p的分量依次是次數(shù)從高 到底各多項(xiàng)式系數(shù)第43頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三用Richard模擬水稻葉伸長生長第44頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三關(guān)于inline函數(shù)例如：y=inline(sin(x)-cos(x),x)輸入y(0)，可得：-1作圖：x=0:0.1:2*pi;plot(x,y(x)第45頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三1、插值問題： 不知道某一函數(shù)f(x)在待定范圍a,b上 的具體表達(dá)式，而只能通過實(shí)驗(yàn)測量得到該 函數(shù)在一系列點(diǎn)ax1, x2 , ., xn b上

25、的值 y0, y1, y2, ., yn，需要找一個(gè)簡單的函數(shù)P(x) 來近似地代替f(x)，要求滿足： P(xi)=yi (i=1,2,.,n)2、插值函數(shù)：P(x) ，3、插值法：求插值函數(shù)P(x)的方法 插 值第46頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三二、常用插值函數(shù)1、多項(xiàng)式函數(shù)2、樣條函數(shù)第47頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三1、多項(xiàng)式插值方法（1）n次代數(shù)插值（2）拉格朗日插值幾點(diǎn)說明：（1）拉格朗日插值基函數(shù)僅與節(jié)點(diǎn)有關(guān)，而 與被插值函數(shù)f(x)無關(guān)。（2）拉格朗日插值多項(xiàng)式僅由數(shù)對(duì)(xi,yi)(i 1,2,n)確定，而與數(shù)對(duì)排列

26、次序無關(guān)。（3）多項(xiàng)式插值除了上述插值法外還有其它 插值法，如newton插值法、hermite插 值法等。第48頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三2、樣條插值方法（1）樣條函數(shù)m次半截冪函數(shù)（2）k次B樣條或k次基本樣條函數(shù)的定義第49頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三（一）廣泛使用的樣條函數(shù)（1）廣泛采用：二次樣條、三次樣條及B樣條。（2）力學(xué)意義： A:二次樣條在力學(xué)上解釋為集中力偶作用 下的彈性細(xì)梁撓度曲線。 B:彈性細(xì)梁受集中載荷作用形成的撓度曲 線，在小撓度的情況下，恰好表示為三 次樣條函數(shù)，集中載荷的作用點(diǎn)，恰好 就是三次樣條函數(shù)的

27、節(jié)點(diǎn)。第50頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三（1）二次樣條的定義 設(shè)a,b 的一個(gè)劃分：a=x0 x1, x2 , ., xn= b，函數(shù)f ( x )各節(jié)點(diǎn)的值分別為： f ( xi )=yi (i=1,2,.,n) 如果二次樣條函數(shù)：滿足： S ( xi )=yi (i=1,2,.,n) 第51頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三（2）三次樣條函數(shù)的定義 設(shè)a,b 的一個(gè)劃分：a=x0 x1, x2 , ., xn= b，函數(shù)f ( x )各節(jié)點(diǎn)的值分別為： f ( xi )=yi (i=1,2,.,n) 如果三次樣條函數(shù)：3滿足： S (

28、xi )=yi (i=1,2,.,n) 第52頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三例：某實(shí)驗(yàn)對(duì)一根長10米的鋼軌進(jìn)行熱源的溫度傳播測試。用x表示測量點(diǎn)0:2.5:10(米)，用h表示測量時(shí)間0:30:60(秒)，用T表示測試所得各點(diǎn)的溫度()。試用線性插值求出在一分鐘內(nèi)每隔20秒、鋼軌每隔1米處的溫度TI。命令如下：x=0:2.5:10;h=0:30:60;T=95,14,0,0,0;88,48,32,12,6;67,64,54,48,41;xi=0:10;hi=0:20:60;TI=interp2(x,h,T,xi,hi)第53頁，共78頁，2022年，5月20日，15

29、點(diǎn)20分，星期三 例：設(shè)z=x2+y2，對(duì)z函數(shù)在0,10,2區(qū)域內(nèi)進(jìn)行插值。 為了說明這個(gè)維數(shù)的插值，再考慮一個(gè)問題。設(shè)人們對(duì)平板上的溫度分布估計(jì)感興趣，給定的溫度值取自平板表面均勻分布的格柵。 采集了下列的數(shù)據(jù)： width=1:5; % index for width of plate (i.e.,the x-dimension) depth=1:3; % index for depth of plate (i,e,the y-dimension) temps=82 81 80 82 84; 79 63 61 65 81; 84 84 82 85 86 % temperature dat

30、a temps = 82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 第54頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 如同在標(biāo)引點(diǎn)上測量一樣，矩陣temps表示整個(gè)平板的溫度分布。temps的列與下標(biāo)depth或y-維相聯(lián)系，行與下標(biāo)width或x-維相聯(lián)系(見圖6)。為了估計(jì)在中間點(diǎn)的溫度，我們必須對(duì)它們進(jìn)行辨識(shí)。 wi=1:0.2:5; % estimate across width of plate d=2; % at a depth of 2 zlinear=interp2(width, depth, temps, wi, d)

31、 ; % linear interpolation zcubic=interp2(width, depth, temps, wi,d, cubic ) ; % cubic interpolation plot(wi, zlinear, - , wi, zcubic) % plot results xlabel( Width of Plate ), ylabel( Degrees Celsius ) title( Temperature at Depth = num2str(d) )第55頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三圖6 在深度d=2處的平板溫度第56頁，共78頁，

32、2022年，5月20日，15點(diǎn)20分，星期三 另一種方法，我們可以在兩個(gè)方向插值。先在三維坐標(biāo)畫出原始數(shù)據(jù)，看一下該數(shù)據(jù)的粗糙程度(見圖7)。 mesh(width, depth, temps) % use mesh plot xlabel( Width of Plate ), ylabel( Depth of Plate ) zlabel( Degrees Celsius ), axis( ij ), grid xi,yi=meshgrid(width, depth);zi=interp2(width, depth, temps, xi, yi, cubic) ;mesh(xi, yi, z

33、i)第57頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三圖7 平板溫度第58頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三然后在兩個(gè)方向上插值，以平滑數(shù)據(jù)(見圖8)。 di=1:0.2:3; % choose higher resolution for depth wi=1:0.2:5; % choose higher resolution for width zcubic=interp2(width, depth, temps, wi, di, cubic ) ; % cubic mesh(wi, di, zcubic) xlabel( Width of Plate

34、), ylabel( Depth of Plate ) zlabel( Degrees Celsius ), axis( ij ), grid第59頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三圖8 二維插值后的平板溫度第60頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 上面的例子清楚地證明了，二維插值更為復(fù)雜，只是因?yàn)橛懈嗟牧恳３指櫋nterp2的基本形式是interp2(x, y, z, xi, yi, method)。這里x和y是兩個(gè)獨(dú)立變量，z是一個(gè)應(yīng)變量矩陣。x和y對(duì)z的關(guān)系是 z(i, :) = f(x, y(i) 和 z(:, j) = f(x(

35、j), y). 也就是，當(dāng)x變化時(shí)，z的第i行與y的第i個(gè)元素y(i)相關(guān)，當(dāng)y變化時(shí)，z的第j列與x的第j個(gè)元素x(j)相關(guān)，。xi是沿x-軸插值的一個(gè)數(shù)值數(shù)組；yi是沿y-軸插值的一個(gè)數(shù)值數(shù)組。第61頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 可選的參數(shù)method可以是 linear，cubic或nearest。在這種情況下，cubic不意味著3次樣條，而是使用3次多項(xiàng)式的另一種算法。linear方法是線性插值，僅用作連接圖上數(shù)據(jù)點(diǎn)。nearest方法只選擇最接近各估計(jì)點(diǎn)的粗略數(shù)據(jù)點(diǎn)。 插值的優(yōu)點(diǎn)： 利用已知點(diǎn)確定未知點(diǎn) 粗糙 精確 集合大的 簡化的第62頁，共78頁，2

36、022年，5月20日，15點(diǎn)20分，星期三例：某觀測站測得某日6:00時(shí)至18:00時(shí)之間每隔2小時(shí)的室內(nèi)外溫度()，用3次樣條插值分別求得該日室內(nèi)外6:30至17:30時(shí)之間每隔2小時(shí)各點(diǎn)的近似溫度()。 設(shè)時(shí)間變量h為一行向量，溫度變量t為一個(gè)兩列矩陣，其中第一列存放室內(nèi)溫度，第二列儲(chǔ)存室外溫度。命令如下：h =6:2:18;t=18,20,22,25,30,28,24;15,19,24,28,34,32,30;XI =6.5:2:17.5YI=interp1(h,t,XI,spline) %用3次樣條插值計(jì)算第63頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 在討論二維插值

37、之前，強(qiáng)調(diào)interp1所強(qiáng)制的二個(gè)強(qiáng)約束是很重要的。 首先，人們不能要求有獨(dú)立變量范圍以外的結(jié)果，例如，interp1(hours, temps, 13.5)導(dǎo)致一個(gè)錯(cuò)誤，因?yàn)閔ours在1到12之間變化。 其次，獨(dú)立變量必須是單調(diào)的。即獨(dú)立變量在值上必須總是增加的或總是減小的。 第64頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三 在我們的例子里，hours是單調(diào)的。然而，如果我們已經(jīng)定義獨(dú)立變量為一天的實(shí)際時(shí)間， time_of_day=7:12 1:6 % start at 7AM,end at 6PM time_of_day = 7 8 9 10 11 12 1 2 3

38、4 5 6 則獨(dú)立變量將不是單調(diào)的，因?yàn)閠ime_of_day增加到12，然后跌到1，再然后增加。如果用time_of_day代替interp1中的hours，將會(huì)返回一個(gè)錯(cuò)誤。同樣的理由，人們不能對(duì)temps插值來找出產(chǎn)生某溫度的時(shí)間(小時(shí))，因?yàn)閠emps不是單調(diào)的。第65頁，共78頁，2022年，5月20日，15點(diǎn)20分，星期三案例3 估計(jì)水箱的水流量模型長度單位：E(3024cm)容積單位：G(=3.785L(升) 某些鎮(zhèn)的用水管理機(jī)構(gòu)需估計(jì)公眾的用水速度(單位是G/h)和每天總用水量的數(shù)據(jù)許多地方?jīng)]有測量流入或流出水箱流量的設(shè)備，而只能測量水箱中的水位(誤差不超過5%)當(dāng)水箱水位低于某最低水位L時(shí)，水泵抽水，灌入水箱，直至水位達(dá)到最高水位H為止。但這也無法測量水泵的流量，因此在水泵啟動(dòng)時(shí)不易建立水箱中水位和水泵工作時(shí)用水量之間關(guān)系。水泵一天灌水12次，每次約2h。試估計(jì)在任意時(shí)刻(包括水泵灌水期間)t流出水箱