數(shù)學建模課件——擬合與插值.pptVIP

數(shù)學建模教程 擬合與插值 在大量的應用領(lǐng)域中 人們經(jīng)常面臨這樣的問題 給定一批數(shù)據(jù)點 需確定滿足特定要求的曲線或曲面 對這個問題有兩種方法 一種是插值法 數(shù)據(jù)假定是正確的 要求以某種方法描述數(shù)據(jù)點之間所發(fā)生的情況 另一種方法是曲線擬合或回歸 人們設法找出某條光滑曲線 它最佳地擬合數(shù)據(jù) 但不必要經(jīng)過任何數(shù)據(jù)點 本專題的主要目的是 了解插值和擬合的基本內(nèi)容 掌握用Matlab求解插值與擬合問題的基本命令 函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個函數(shù)作為近似 由于近似的要求不同 二者的數(shù)學方法上是完全不同的 內(nèi)容提綱 1 擬合問題引例及基本理論2 Matlab求解擬合問題3 應用實例4 插值問題引例及基本理論5 Maltab求解插值問題6 應用實例 擬合問題引例1 溫度t 0C 20 532 751 073 095 7電阻R 7658268739421032 已知熱敏電阻數(shù)據(jù) 求600C時的電阻R 設R at ba b為待定系數(shù) 一 擬合問題 擬合問題引例2 求血藥濃度隨時間的變化規(guī)律c t 作半對數(shù)坐標系 semilogy 下的圖形 曲線擬合問題的提法 已知一組 二維 數(shù)據(jù) 即平面上n個點 xi yi i 1 n 尋求一個函數(shù) 曲線 y f x 使f x 在某種準則下與所有數(shù)據(jù)點最為接近 即曲線擬合得最好 y f x i為點 xi yi 與曲線y f x 的距離 線性最小二乘擬合f x a1r1 x amrm x 中函數(shù) r1 x rm x 的選取 1 通過機理分析建立數(shù)學模型來確定f x 2 將數(shù)據(jù) xi yi i 1 n作圖 通過直觀判斷確定f x 曲線擬合問題最常用的解法 線性最小二乘法的基本思路 第一步 先選定一組函數(shù)r1 x r2 x rm x m n 令f x a1r1 x a2r2 x amrm x 1 其中a1 a2 am為待定系數(shù) 第二步 確定a1 a2 am的準則 最小二乘準則 使n個點 xi yi 與曲線y f x 的距離 i的平方和最小 記 問題歸結(jié)為 求a1 a2 am使J a1 a2 am 最小 線性最小二乘法的求解 預備知識 超定方程組 方程個數(shù)大于未知量個數(shù)的方程組 超定方程一般是不存在解的矛盾方程組 如果有向量a使得達到最小 則稱a為上述超定方程的最小二乘解 線性最小二乘法的求解 定理 當RTR可逆時 超定方程組 3 存在最小二乘解 且即為方程組RTRa RTy 正則 正規(guī) 方程組的解 a RTR 1RTy 所以 曲線擬合的最小二乘法要解決的問題 實際上就是求以下超定方程組的最小二乘解的問題 用MATLAB解擬合問題 1 線性最小二乘擬合 2 非線性最小二乘擬合 用MATLAB作線性最小二乘擬合 1 作多項式f x a1xm amx am 1擬合 可利用已有命令 a polyfit x y m 3 對超定方程組 2 多項式在x處的值y的計算命令 y polyval a x 例對下面一組數(shù)據(jù)作二次多項式擬合 1 輸入命令 x 0 0 1 1 y 0 447 1 978 3 28 6 16 7 08 7 34 7 66 9 56 9 48 9 30 11 2 R x 2 x ones 11 1 A R y MATLAB zxec1 解法1 解超定方程的方法 2 計算結(jié)果 9 8108 20 1293 0 0317 2 計算結(jié)果 9 8108 20 1293 0 0317 解法2 用多項式擬合的命令 MATLAB zxec2 1 輸入命令 x 0 0 1 1 y 0 447 1 978 3 28 6 16 7 08 7 34 7 66 9 56 9 48 9 30 11 2 A polyfit x y 2 z polyval A x plot x y k x z r 作出數(shù)據(jù)點和擬合曲線的圖形 1 lsqcurvefit已知數(shù)據(jù)點 xdata xdata1 xdata2 xdatan ydata ydata1 ydata2 ydatan 用MATLAB作非線性最小二乘擬合 兩個求非線性最小二乘擬合的函數(shù) lsqcurvefit lsqnonlin 相同點和不同點 兩個命令都要先建立M 文件fun m 定義函數(shù)f x 但定義f x 的方式不同 lsqcurvefit用以求含參量x 向量 的向量值函數(shù)F x xdata F x xdata1 F x xdatan T使得 輸入格式 1 x lsqcurvefit fun x0 xdata ydata 2 x lsqcurvefit fun x0 xdata ydata lb ub 3 x lsqcurvefit fun x0 xdata ydata lb ub options 4 x resnorm lsqcurvefit fun x0 xdata ydata 5 x resnorm residual lsqcurvefit fun x0 xdata ydata 說明 x lsqcurvefit fun x0 xdata ydata options lsqnonlin用以求含參量x 向量 的向量值函數(shù)f x f1 x f2 x fn x T 使得最小 其中fi x f x xdatai ydatai F x xdatai ydatai 2 lsqnonlin 已知數(shù)據(jù)點 xdata xdata1 xdata2 xdatan ydata ydata1 ydata2 ydatan 輸入格式 1 x lsqnonlin fun x0 2 x lsqnonlin fun x0 lb ub 3 x lsqnonlin fun x0 lb ub options 4 x resnorm lsqnonlin fun x0 5 x resnorm residual lsqnonlin fun x0 說明 x lsqnonlin fun x0 options 例2用下面一組數(shù)據(jù)擬合中的參數(shù)a b k 該問題即解的最優(yōu)化問題 1 編寫M 文件curvefun1 mfunctionf curvefun1 x tdata f x 1 x 2 exp 0 02 x 3 tdata 其中x 1 a x 2 b x 3 k 2 輸入命令tdata 100 100 1000cdata 1e 03 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 lsqcurvefit curvefun1 x0 tdata cdata f curvefun1 x tdata F x tdata x a b k 解法1 用命令lsqcurvefit 3 運算結(jié)果 f 0 00430 00510 00560 00590 00610 00620 00620 00630 00630 0063x 0 0063 0 00340 2542 4 結(jié)論 a 0 0063 b 0 0034 k 0 2542 1 編寫M 文件curvefun2 mfunctionf curvefun2 x tdata 100 100 1000 cdata 1e 03 4 54 4 99 5 35 5 65 5 90 6 10 6 26 6 39 6 50 6 59 f x 1 x 2 exp 0 02 x 3 tdata cdata 2 輸入命令 x0 0 2 0 05 0 05 x lsqnonlin curvefun2 x0 f curvefun2 x 函數(shù)curvefun2的自變量是x cdata和tdata是已知參數(shù) 故應將cdatatdata的值寫在curvefun2 m中 解法2用命令lsqnonlin x a b k 3 運算結(jié)果為f 1 0e 003 0 2322 0 1243 0 2495 0 2413 0 1668 0 07240 02410 11590 20300 2792 x 0 0063 0 00340 2542 可以看出 兩個命令的計算結(jié)果是相同的 4 結(jié)論 即擬合得a 0 0063b 0 0034k 0 2542 插值問題 擬合與插值的關(guān)系 說明 函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個函數(shù)作為近似 由于近似的要求不同 二者的數(shù)學方法上完全不同 實例 下面數(shù)據(jù)是某次實驗所得 希望得到x和f之間的關(guān)系 MATLAB cn 問題 給定一批數(shù)據(jù)點 需確定滿足特定要求的曲線或曲面 解決方案 若不要求曲線 面 通過所有數(shù)據(jù)點 而是要求它反映對象整體的變化趨勢 就是數(shù)據(jù)擬合 又稱曲線擬合或曲面擬合 若要求所求曲線 面 通過所給所有數(shù)據(jù)點 就是插值問題 最臨近插值 線性插值 樣條插值與曲線擬合結(jié)果 拉格朗日插值 分段線性插值 一維插值 一 插值的定義 二 插值的方法 三 用Matlab解插值問題 返回 二維插值 一 二維插值定義 二 網(wǎng)格節(jié)點插值法 三 用Matlab解插值問題 最鄰近插值 分片線性插值 雙線性插值 網(wǎng)格節(jié)點數(shù)據(jù)的插值 散點數(shù)據(jù)的插值 一維插值的定義 返回 拉格朗日插值問題 的函數(shù)值 已知y f x 在n 1個點 構(gòu)造n次多項式pn x 使得 從而得到f x 的近似計算式 一 線性插值 n 1 一次插值 求解L1 x a1x a0 已知 使得L1 xi yi i 0 1 如果令 則稱l0 x l1 x 為x0 x1上的線性插值基函數(shù) 這時 f x L1 x y0l0 x y1l1 x 二 拋物線插值 n 2 二次插值 求解L2 x a2x2 a1x a0 使得L2 xi yi i 0 1 2 關(guān)于二次多項式的構(gòu)造采用如下方法 令 已知 并由插值條件 得到 L2 x A x x1 x x2 B x x0 x x2 C x x0 x x1 L2 x0 y0 L2 x1 y1 L2 x y2 于是得到 這時 f x L2 x y0l0 x y1l1 x y2l2 x 如果令 則有 并稱其為二次Lagrange插值多項式 基函數(shù)表示 x0 x1 x2上的二次插值基函數(shù) 稱為拉格朗日插值基函數(shù) 已知函數(shù)f x 在n 1個點x0 x1 xn處的函數(shù)值為y0 y1 yn 求一n次多項式函數(shù)Pn x 使其滿足 Pn xi yi i 0 1 n 解決此問題的拉格朗日插值多項式公式如下 其中Li x 為n次多項式 拉格朗日 Lagrange 插值 拉格朗日 Lagrange 插值 特別地 兩點一次 線性 插值多項式 三點二次 拋物 插值多項式 拉格朗日多項式插值的這種振蕩現(xiàn)象叫Runge現(xiàn)象 采用拉格朗日多項式插值 選取不同插值節(jié)點個數(shù)n 1 其中n為插值多項式的次數(shù) 當n分別取2 4 6 8 10時 繪出插值結(jié)果圖形 例 分段線性插值 計算量與n無關(guān) n越大 誤差越小 例 用分段線性插值法求插值 并觀察插值誤差 在 6 6 中平均選取41個點作插值 結(jié)果如圖示 比分段線性插值更光滑 在數(shù)學上 光滑程度的定量描述是 函數(shù) 曲線 的k階導數(shù)存在且連續(xù) 則稱該曲線具有k階光滑性 光滑性的階次越高 則越光滑 是否存在較低次的分段多項式達到較高階光滑性的方法 三次樣條插值就是一個很好的例子 三次樣條插值 g x 為被插值函數(shù) 三次樣條插值 例 用三次樣條插值選取11個基點計算插值 用MATLAB作插值計算 一維插值函數(shù) yi interp1 x y xi method nearest 最鄰近插值 linear 線性插值 spline 三次樣條插值 cubic 立方插值 缺省時 分段線性插值 注意 所有的插值方法都要求x是單調(diào)的 并且xi不能夠超過x的范圍 例 在1 12的11小時內(nèi) 每隔1小時測量一次溫度 測得的溫度依次為 5 8 9 15 25 29 31 30 22 25 27 24 試估計每隔1 10小時的溫度值 hours 1 12 temps 589152529313022252724 h 1 0 1 12 t interp1 hours temps h spline 直接輸出數(shù)據(jù)將是很多的 plot hours temps h t hours temps r 作圖xlabel Hour ylabel DegreesCelsius 例已知飛機下輪廓線上數(shù)據(jù)如下 求x每改變0 1時的y值 二維插值的定義 第一種 網(wǎng)格節(jié)點 已知m n個節(jié)點 第二種 散亂節(jié)點 注意 最鄰近插值一般不連續(xù) 具有連續(xù)性的最簡單的插值是分片線性插值 最鄰近插值 二維或高維情形的最鄰近插值 與被插值點最鄰近的節(jié)點的函數(shù)值即為所求 將四個插值點 矩形的四個頂點 處的函數(shù)值依次簡記為 分片線性插值 f xi yj f1 f xi 1 yj f2 f xi 1 yj 1 f3 f xi yj 1 f4 插值函數(shù)為 第二片 上三角形區(qū)域 x y 滿足 插值函數(shù)為 注意 x y 當然應該是在插值節(jié)點所形成的矩形區(qū)域內(nèi) 顯然 分片線性插值函數(shù)是連續(xù)的 分兩片的函數(shù)表達式如下 第一片 下三角形區(qū)域 x y 滿足 雙線性插值是一片一片的空間二次曲面構(gòu)成 雙線性插值函數(shù)的形式如下 其中有四個待定系數(shù) 利用該函數(shù)在矩形的四個頂點 插值節(jié)點 的函數(shù)值 得到四個代數(shù)方程 正好確定四個系數(shù) 雙線性插值 要求x0 y0單調(diào) x y可取為矩陣 或x取行向量 y取為列向量 x y的值分別不能超出x0 y0的范圍 z interp2 x0 y0 z0 x y method 用MATLAB作網(wǎng)格節(jié)點數(shù)據(jù)的插值 nearest 最鄰近插值 linear 雙線性插值 cubic 雙三次插值缺省時 雙線性插值 例 測得平板表面3 5網(wǎng)格點處的溫度分別為 828180828479636165818484828586試作出平板表面的溫度分布曲面z f x y 的圖形 輸入以下命令 x 1 5 y 1 3 temps 8281808284 7963616581 8484828586 mesh x y temps 1 先在三維坐標畫出原始數(shù)據(jù) 畫出粗糙的溫度分布曲圖 再輸入以下命令 xi 1 0 2 5 yi 1 0 2 3 zi interp2 x y temps xi yi cubic mesh xi yi zi 畫出插值后的溫度分布曲面圖 2 以平滑數(shù)據(jù) 在x y方向上每隔0 2個單位的地方進行插值 通過此例對最近鄰點插值 雙線性插值方法和雙三次插值方法的插值效果進行比較 原始數(shù)據(jù)圖 最鄰近插值 雙線性插值 雙三次插值 等高線圖 插值函數(shù)griddata格式為 cz griddata x y z cx cy method 用MATLAB作散點數(shù)據(jù)的插值計算 要求cx取行向量 cy取為列向量 nearest 最鄰近插值 linear 雙線性插值 cubic 雙三次插值 v4 Matlab提供的插值方法缺省時 雙線性插值 一室模型 將整個機體看作一個房室 稱中心室 室內(nèi)血藥濃度是均勻的 快速靜脈注射后 濃度立即上升 然后迅速下降 當濃度太低時 達不到預期的治療效果 當濃度太高 又可能導致藥物中毒或副作用太強 臨床上 每種藥物有一個最小有效濃度c1和一個最大有效濃度c2 設計給藥方案時 要使血藥濃度保持在c1 c2之間 本題設c1 10 c2 25 ug ml 一種新藥用于臨床之前 必須設計給藥方案 藥物進入機體后血液輸送到全身 在這個過程中不斷地被吸收 分布 代謝 最終排出體外 藥物在血液中的濃度 即單位體積血液中的藥物含量 稱為血藥濃度 要設計給藥方案 必須知道給藥后血藥濃度隨時間變化的規(guī)律 從實驗和理論兩方面著手 給藥方案 1 在快速靜脈注射的給藥方式下 研究血藥濃度 單位體積血液中的藥物含量 的變化規(guī)律 t 問題 2 給定藥物的最小有效濃度和最大治療濃度 設計給藥方案 每次注射劑量多大 間隔時間多長 分析 理