MATLAB-第14講_擬合.ppt

1 數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn) 經(jīng)濟(jì)數(shù)學(xué)系數(shù)學(xué)建模研究室 擬合 2 實(shí)驗(yàn)?zāi)康?實(shí)驗(yàn)內(nèi)容 2 掌握用數(shù)學(xué)軟件求解擬合問題 1 直觀了解擬合基本內(nèi)容 1 擬合問題引例及基本理論 4 實(shí)驗(yàn)作業(yè) 2 用數(shù)學(xué)軟件求解擬合問題 3 應(yīng)用實(shí)例 3 擬合 2 擬合的基本原理 1 擬合問題引例 4 擬合問題引例1 求600C時(shí)的電阻R 設(shè)R at ba b為待定系數(shù) 5 擬合問題引例2 求血藥濃度隨時(shí)間的變化規(guī)律c t 作半對數(shù)坐標(biāo)系 semilogy 下的圖形 MATLAB aa1 6 曲線擬合問題的提法 已知一組 二維 數(shù)據(jù) 即平面上n個(gè)點(diǎn) xi yi i 1 n 尋求一個(gè)函數(shù) 曲線 y f x 使f x 在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近 即曲線擬合得最好 y f x i為點(diǎn) xi yi 與曲線y f x 的距離 7 擬合與插值的關(guān)系 函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近似 由于近似的要求不同 二者的數(shù)學(xué)方法上是完全不同的 實(shí)例 下面數(shù)據(jù)是某次實(shí)驗(yàn)所得 希望得到X和f之間的關(guān)系 MATLAB cn 問題 給定一批數(shù)據(jù)點(diǎn) 需確定滿足特定要求的曲線或曲面 解決方案 若不要求曲線 面 通過所有數(shù)據(jù)點(diǎn) 而是要求它反映對象整體的變化趨勢 這就是數(shù)據(jù)擬合 又稱曲線擬合或曲面擬合 若要求所求曲線 面 通過所給所有數(shù)據(jù)點(diǎn) 就是插值問題 8 最臨近插值 線性插值 樣條插值與曲線擬合結(jié)果 9 曲線擬合問題最常用的解法 線性最小二乘法的基本思路 第一步 先選定一組函數(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的準(zhǔn)則 最小二乘準(zhǔn)則 使n個(gè)點(diǎn) xi yi 與曲線y f x 的距離 i的平方和最小 記 問題歸結(jié)為 求a1 a2 am使J a1 a2 am 最小 10 線性最小二乘法的求解 預(yù)備知識(shí) 超定方程組 方程個(gè)數(shù)大于未知量個(gè)數(shù)的方程組 超定方程一般是不存在解的矛盾方程組 如果有向量a使得達(dá)到最小 則稱a為上述超定方程的最小二乘解 11 線性最小二乘法的求解 定理 當(dāng)RTR可逆時(shí) 超定方程組 3 存在最小二乘解 且即為方程組RTRa RTy的解 a RTR 1RTy 所以 曲線擬合的最小二乘法要解決的問題 實(shí)際上就是求以下超定方程組的最小二乘解的問題 12 線性最小二乘擬合f x a1r1 x amrm x 中函數(shù) r1 x rm x 的選取 1 通過機(jī)理分析建立數(shù)學(xué)模型來確定f x 2 將數(shù)據(jù) xi yi i 1 n作圖 通過直觀判斷確定f x 13 用MATLAB解擬合問題 1 線性最小二乘擬合 2 非線性最小二乘擬合 14 用MATLAB作線性最小二乘擬合 1 作多項(xiàng)式f x a1xm amx am 1擬合 可利用已有程序 a polyfit x y m 2 對超定方程組 3 多項(xiàng)式在x處的值y可用以下命令計(jì)算 y polyval a x 15 例對下面一組數(shù)據(jù)作二次多項(xiàng)式擬合 16 1 輸入以下命令 x 0 0 1 1 y 0 4471 9783 286 167 087 347 669 569 489 3011 2 R x 2 x ones 11 1 A R y MATLAB zxec1 解法1 用解超定方程的方法 2 計(jì)算結(jié)果 9 810820 1293 0 0317 17 1 輸入以下命令 x 0 0 1 1 y 0 4471 9783 286 167 087 347 669 569 489 3011 2 A polyfit x y 2 z polyval A x plot x y k x z r 作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形 2 計(jì)算結(jié)果 9 810820 1293 0 0317 解法2 用多項(xiàng)式擬合的命令 MATLAB zxec2 18 1 lsqcurvefit已知數(shù)據(jù)點(diǎn) xdata xdata1 xdata2 xdatan ydata ydata1 ydata2 ydatan 用MATLAB作非線性最小二乘擬合 Matlab的提供了兩個(gè)求非線性最小二乘擬合的函數(shù) lsqcurvefit和lsqnonlin 兩個(gè)命令都要先建立M 文件fun m 在其中定義函數(shù)f x 但兩者定義f x 的方式是不同的 可參考例題 lsqcurvefit用以求含參量x 向量 的向量值函數(shù)F x xdata F x xdata1 F x xdatan T中的參變量x 向量 使得 19 輸入格式為 1 x lsqcurvefit fun x0 xdata ydata 2 x lsqcurvefit fun x0 xdata ydata options 3 x lsqcurvefit fun x0 xdata ydata options grad 4 x options lsqcurvefit fun x0 xdata ydata 5 x options funval lsqcurvefit fun x0 xdata ydata 6 x options funval Jacob lsqcurvefit fun x0 xdata ydata 說明 x lsqcurvefit fun x0 xdata ydata options 20 lsqnonlin用以求含參量x 向量 的向量值函數(shù)f x f1 x f2 x fn x T中的參量x 使得最小 其中fi x f x xdatai ydatai F x xdatai ydatai 2 lsqnonlin 已知數(shù)據(jù)點(diǎn) xdata xdata1 xdata2 xdatan ydata ydata1 ydata2 ydatan 21 輸入格式為 1 x lsqnonlin fun x0 2 x lsqnonlin fun x0 options 3 x lsqnonlin fun x0 options grad 4 x options lsqnonlin fun x0 5 x options funval lsqnonlin fun x0 說明 x lsqnonlin fun x0 options 22 例2用下面一組數(shù)據(jù)擬合中的參數(shù)a b k 該問題即解最優(yōu)化問題 23 MATLAB fzxec1 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 24 3 運(yùn)算結(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 25 MATLAB fzxec2 解法2用命令lsqnonlinf x F x tdata ctada x a b k 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ù) 故應(yīng)將cdatatdata的值寫在curvefun2 m中 26 3 運(yùn)算結(jié)果為f 1 0e 003 0 2322 0 1243 0 2495 0 2413 0 1668 0 07240 02410 11590 20300 2792x 0 0063 0 00340 2542 可以看出 兩個(gè)命令的計(jì)算結(jié)果是相同的 4 結(jié)論 即擬合得a 0 0063b 0 0034k 0 2542 27 MATLAB解應(yīng)用問題實(shí)例 1 電阻問題 2 給藥方案問題 3 水塔流量估計(jì)問題 28 MATLAB dianzu1 電阻問題 得到a1 3 3940 a2 702 4918 方法2 直接用 結(jié)果相同 MATLAB dianzu2 29 一室模型 將整個(gè)機(jī)體看作一個(gè)房室 稱中心室 室內(nèi)血藥濃度是均勻的 快速靜脈注射后 濃度立即上升 然后迅速下降 當(dāng)濃度太低時(shí) 達(dá)不到預(yù)期的治療效果 當(dāng)濃度太高 又可能導(dǎo)致藥物中毒或副作用太強(qiáng) 臨床上 每種藥物有一個(gè)最小有效濃度c1和一個(gè)最大有效濃度c2 設(shè)計(jì)給藥方案時(shí) 要使血藥濃度保持在c1 c2之間 本題設(shè)c1 10 c2 25 ug ml 一種新藥用于臨床之前 必須設(shè)計(jì)給藥方案 藥物進(jìn)入機(jī)體后血液輸送到全身 在這個(gè)過程中不斷地被吸收 分布 代謝 最終排出體外 藥物在血液中的濃度 即單位體積血液中的藥物含量 稱為血藥濃度 30 要設(shè)計(jì)給藥方案 必須知道給藥后血藥濃度隨時(shí)間變化的規(guī)律 從實(shí)驗(yàn)和理論兩方面著手 給藥方案 1 在快速靜脈注射的給藥方式下 研究血藥濃度 單位體積血液中的藥物含量 的變化規(guī)律 t 問題 2 給定藥物的最小有效濃度和最大治療濃度 設(shè)計(jì)給藥方案 每次注射劑量多大 間隔時(shí)間多長 分析 理論 用一室模型研究血藥濃度變化規(guī)律 實(shí)驗(yàn) 對血藥濃度數(shù)據(jù)作擬合 符合負(fù)指數(shù)變化規(guī)律 3 血液容積v t 0注射劑量d 血藥濃度立即為d v 2 藥物排除速率與血藥濃度成正比 比例系數(shù)k 0 模型假設(shè) 1 機(jī)體看作一個(gè)房室 室內(nèi)血藥濃度均勻 一室模型 模型建立 在此 d 300mg t及c t 在某些點(diǎn)處的值見前表 需經(jīng)擬合求出參數(shù)k v 用線性最小二乘擬合c t MATLAB lihe1 計(jì)算結(jié)果 用非線性最小二乘擬合c t 給藥方案設(shè)計(jì) 設(shè)每次注射劑量D 間隔時(shí)間 血藥濃度c t 應(yīng)c1 c t c2 初次劑量D0應(yīng)加大 給藥方案記為 2 1 計(jì)算結(jié)果 給藥方案 c1 10 c2 25k 0 2347v 15 02 35 故可制定給藥方案 即 首次注射375mg 其余每次注射225mg 注射的間隔時(shí)間為4小時(shí) 36 估計(jì)水塔的流量 2 解題思路 3 算法設(shè)計(jì)與編程 1 問題 37 某居民區(qū)有一供居民用水的園柱形水塔 一般可以通過測量其水位來估計(jì)水的流量 但面臨的困難是 當(dāng)水塔水位下降到設(shè)定的最低水位時(shí) 水泵自動(dòng)啟動(dòng)向水塔供水 到設(shè)定的最高水位時(shí)停止供水 這段時(shí)間無法測量水塔的水位和水泵的供水量 通常水泵每天供水一兩次 每次約兩小時(shí) 水塔是一個(gè)高12 2米 直徑17 4米的正園柱 按照設(shè)計(jì) 水塔水位降至約8 2米時(shí) 水泵自動(dòng)啟動(dòng) 水位升到約10 8米時(shí)水泵停止工作 表1是某一天的水位測量記錄 試估計(jì)任何時(shí)刻 包括水泵正供水時(shí) 從水塔流出的水流量 及一天的總用水量 38 39 流量估計(jì)的解題思路 擬合水位 時(shí)間函數(shù) 確定流量 時(shí)間函數(shù) 估計(jì)一天總用水量 40 擬合水位 時(shí)間函數(shù)測量記錄看 一天有兩個(gè)供水時(shí)段 以下稱第1供水時(shí)段和第2供水時(shí)段 和3個(gè)水泵不工作時(shí)段 以下稱第1時(shí)段t 0到t 8 97 第2次時(shí)段t 10 95到t 20 84和第3時(shí)段t 23以后 對第1 2時(shí)段的測量數(shù)據(jù)直接分別作多項(xiàng)式擬合 得到水位函數(shù) 為使擬合曲線比較光滑 多項(xiàng)式次數(shù)不要太高 一般在3 6 由于第3時(shí)段只有3個(gè)測量記錄 無法對這一時(shí)段的水位作出較好的擬合 41 2 確定流量 時(shí)間函數(shù)對于第1 2時(shí)段只需將水位函數(shù)求導(dǎo)數(shù)即可 對于兩個(gè)供水時(shí)段的流量 則用供水時(shí)段前后 水泵不工作時(shí)段 的流量擬合得到 并且將擬合得到的第2供水時(shí)段流量外推 將第3時(shí)段流量包含在第2供水時(shí)段內(nèi) 42 3 一天總用水量的估計(jì)總用水量等于兩個(gè)水泵不工作時(shí)段和兩個(gè)供水時(shí)段用水量之和 它們都可以由流量對時(shí)間的積分得到 43 算法設(shè)計(jì)與編程 1 擬合第1 2時(shí)段的水位 并導(dǎo)出流量 2 擬合供水時(shí)段的流量 3 估計(jì)一天總用水量 4 流量及總用水量的檢驗(yàn) 44 1 擬合第1時(shí)段的水位 并導(dǎo)出流量設(shè)t h為已輸入的時(shí)刻和水位測量記錄 水泵啟動(dòng)的4個(gè)時(shí)刻不輸入 第1時(shí)段各時(shí)刻的流量可如下得 1 c1 polyfit t 1 10 h 1 10 3 用3次多項(xiàng)式擬合第1時(shí)段水位 c1輸出3次多項(xiàng)式的系數(shù)2 a1 polyder c1 a1輸出多項(xiàng)式 系數(shù)為c1 導(dǎo)數(shù)的系數(shù)3 tp1 0 0 1 9 x1 polyval a1 tp1 x1輸出多項(xiàng)式 系數(shù)為a1 在tp1點(diǎn)的函數(shù)值 取負(fù)后邊為正值 即tp1時(shí)刻的流量 MATLAB llgj1 4 流量函數(shù)為 45 2 擬合第2時(shí)段的水位 并導(dǎo)出流量設(shè)t h為已輸入的時(shí)刻和水位測量記錄 水泵啟動(dòng)的4個(gè)時(shí)刻不輸入 第2時(shí)段各時(shí)刻的流量可如下得 1 c2 polyfit t 10 9 21 h 10 9 21 3 用3次多項(xiàng)式擬合第2時(shí)段水位 c2輸出3次多項(xiàng)式的系數(shù)2 a2 polyder c2 a2輸出多項(xiàng)式 系數(shù)為c2 導(dǎo)數(shù)的系數(shù)3 tp2 10 9 0 1 21 x2 polyval a2 tp2 x2輸出多項(xiàng)式 系數(shù)為a2 在tp2點(diǎn)的函數(shù)值 取負(fù)后邊為正值 即tp2時(shí)刻的流量 MATLAB llgj2 4 流量函數(shù)為 46 3 擬合供水時(shí)段的流量在第1供水時(shí)段 t 9 11 之前 即第1時(shí)段 和之后 即第2時(shí)段 各取幾點(diǎn) 其流量已經(jīng)得到 用它們擬合第1供水時(shí)段的流量 為使流量函數(shù)在t 9和t 11連續(xù) 我們簡單地只取4個(gè)點(diǎn) 擬合3次多項(xiàng)式 即曲線必過這4個(gè)點(diǎn) 實(shí)現(xiàn)如下 xx1 polyval a1 89 取第1時(shí)段在t 8 9的流量xx2 polyval a2 1112 取第2時(shí)段在t 11 12的流量xx12 xx1xx2 c12 polyfit 891112 xx12 3 擬合3次多項(xiàng)式tp12 9 0 1 11 x12 polyval c12 tp12 x12輸出第1供水時(shí)段各時(shí)刻的流量 MATLAB llgj3 擬合的流量函數(shù)為 47 在第2供水時(shí)段之前取t 20 20 8兩點(diǎn)的流水量 在該時(shí)刻之后 第3時(shí)段 僅有3個(gè)水位記錄 我們用差分得到流量 然后用這4個(gè)數(shù)值擬合第2供水時(shí)段的流量如下 dt3 diff t 22 24 最后3個(gè)時(shí)刻的兩兩之差dh3 diff h 22 24 最后3個(gè)水位的兩兩之差dht3 dh3 dt3 t 22 和t 23 的流量t3 2020 8t 22 t 23 xx3 polyval a2 t3 1 2 dht3 取t3各時(shí)刻的流量c3 polyfit t3 xx3 3 擬合3次多項(xiàng)式t3 20 8 0 1 24 x3 polyval c3 tp3 x3輸出第2供水時(shí)段 外推至t 24 各時(shí)刻的流量 MATLAB llgj4 擬合的流量函數(shù)為 48 3 一天總用水量的估計(jì)第1 2時(shí)段和第1 2供水時(shí)段流量的積分之和 就是一天總用水量 雖然諸時(shí)段的流量已表為多項(xiàng)式函數(shù) 積分可以解析地算出 這里仍用數(shù)值積分計(jì)算如下 y1 0 1 trapz x1 第1時(shí)段用水量 仍按高度計(jì) 0 1為積分步長y2 0 1 trapz x2 第2時(shí)段用水量y12 0 1 trapz x12 第1供水時(shí)段用水量y3 0 1 trapz x3 第2供水時(shí)段用水量y y1 y2 y12 y3 237 8 0 01 一天總用水量 計(jì)算結(jié)果 y1 146 2 y2 266 8 y12 47 4 y3 77 3 y 1250 4 MATLAB llgjz 49 4 流量及總用水量的檢驗(yàn)計(jì)算出的各時(shí)刻的流量可用水位記錄的數(shù)值微分來檢驗(yàn) 用水量y1可用第1時(shí)段水位測量記錄中下降高度968 822 146來檢驗(yàn) 類似地 y2用1082 822 260檢驗(yàn) 供水時(shí)段流量的一種檢驗(yàn)方法如下 供水時(shí)段的用水量加上水位上升值260是該時(shí)段泵入的水量 除以時(shí)段長度得到水泵的功率 單位時(shí)間泵入的水量