數學建模與實驗嚴可頌主講數學建模與實驗嚴可頌主講1教學項目四數值分析法模型教學項目四數值分析法模型24.1插值法建模

拉格朗日插值分段線性插值三次樣條插值一、插值的定義二、插值的方法三、用Matlab解插值問題4.1插值法建模拉格朗日插值分段線性插值三次樣條插值一、3已知n+1個節點其中互不相同，不妨設求任一插值點處的插值節點可視為由產生，表達式復雜，或無封閉形式，或未知。4.1.1插值法的定義

已知n+1個節點其中互不相同，不妨設求任一插值點處的插值4構造一個(相對簡單的)函數通過全部節點,即再用計算插值，即構造一個(相對簡單的)函數通過全部節點,即再用計算插值，54.1.2插值的方法

已知函數f(x)在n+1個點x0,x1,…,xn處的函數值為y0,y1,…,yn

。求一n次多項式函數Pn(x)，使其滿足：拉格朗日(Lagrange)插值4.1.2插值的方法已知函數f(x)在n+16稱為拉格朗日插值基函數。解決此問題的拉格朗日插值多項式公式如下其中Li(x)為n次多項式：稱為拉格朗日插值基函數。解決此問題的拉格朗日插值多項式公7特殊情況兩點一次(線性)插值多項式:三點二次(拋物)插值多項式:特殊情況兩點一次(線性)插值多項式:三點二次(拋物)插值多項8采用拉格朗日多項式插值：選取不同插值節點個數n+1，其中n為插值多項式的次數，當n取10時，繪出插值結果圖形.例采用拉格朗日多項式插值：選取不同插值節點個數9分段線性插值計算量與n無關;n越大，誤差越小.xjxj-1xj+1x0xnxoy分段線性插值計算量與n無關;xjxj-1xj+110比分段線性插值更光滑。xyxi-1xiab在數學上，光滑程度的定量描述是：函數(曲線)的k階導數存在且連續，則稱該曲線具有k階光滑性。光滑性的階次越高，則越光滑。是否存在較低次的分段多項式達到較高階光滑性的方法？三次樣條插值就是一個很好的例子。三次樣條插值比分段線性插值更光滑。xyxi-111

三次樣條插值g(x)為被插值函數。三次樣條插值g(x)為被插值函數。12用MATLAB作插值計算一維插值函數：yi=interp1(x，y，xi，'method')插值方法被插值點插值節點xi處的插值結果‘nearest’

：最鄰近插值‘linear’

：線性插值；‘spline’

：三次樣條插值；‘cubic’

：立方插值。缺省時：分段線性插值。注意：所有的插值方法都要求x是單調的，并且xi不能夠超過x的范圍。用MATLAB作插值計算一維插值函數：yi=interp1(13例：在1-12的11小時內，每隔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');(直接輸出數據將是很多的)plot(hours,temps,'+',h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius’)例：在1-12的11小時內，每隔1小時測量一次溫度，14xy實驗習題

已知飛機下輪廓線上數據如下，求x每改變0.1時的y值。機翼下輪廓線xy實驗習題機翼下輪廓線154.2擬合法簡介

2.擬合的基本原理1.擬合問題引例4.2擬合法簡介2.擬合的基本原理1.擬合問題引例16擬合問題引例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數據：求600C時的電阻R。

設

R=at+ba,b為待定系數擬合問題引例1溫度t(0C)20.517擬合問題引例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速靜脈注射下的血藥濃度數據(t=0注射300mg)求血藥濃度隨時間的變化規律c(t).作半對數坐標系(semilogy)下的圖形擬合問題引例2t(h)18曲線擬合問題的提法已知一組（二維）數據，即平面上n個點（xi,yi)i=1,…n,

尋求一個函數（曲線）y=f(x),

使f(x)

在某種準則下與所有數據點最為接近，即曲線擬合得最好。

+++++++++xyy=f(x)(xi,yi)ii為點（xi,yi)與曲線y=f(x)的距離曲線擬合問題的提法已知一組（二維）數據，即平19擬合與插值的關系

函數插值與曲線擬合都是要根據一組數據構造一個函數作為近似，由于近似的要求不同，二者的數學方法上是完全不同的。實例：下面數據是某次實驗所得，希望得到X和f之間的關系？問題：給定一批數據點，需確定滿足特定要求的曲線或曲面解決方案：若不要求曲線（面）通過所有數據點，而是要求它反映對象整體的變化趨勢，這就是數據擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）通過所給所有數據點，就是插值問題；擬合與插值的關系函數插值與曲線擬合都是要根據一組20最臨近插值、線性插值、樣條插值與曲線擬合結果：最臨近插值、線性插值、樣條插值與曲線擬合結果：21曲線擬合問題最常用的解法——線性最小二乘法的基本思路第一步:先選定一組函數

r1(x),r2(x),…rm(x),m<n,

令

f(x)=a1r1(x)+a2r2(x)+…+amrm(x)（1）其中

a1,a2,…am

為待定系數。

第二步:確定a1,a2,…am

的準則（最小二乘準則）：使n個點（xi,yi)與曲線y=f(x)的距離i的平方和最小

。記

問題歸結為，求

a1,a2,…am

使

J(a1,a2,…am)

最小。曲線擬合問題最常用的解法——線性最小二乘法的基本思路第一步:22線性最小二乘法的求解：預備知識超定方程組：方程個數大于未知量個數的方程組即Ra=y其中超定方程一般是不存在解的矛盾方程組。如果有向量a使得達到最小，則稱a為上述超定方程的最小二乘解。線性最小二乘法的求解：預備知識超定方程組：方程個數大于未知量23線性最小二乘法的求解

定理：當RTR可逆時，超定方程組（3）存在最小二乘解，且即為方程組

RTRa=RTy的解：a=(RTR)-1RTy

所以，曲線擬合的最小二乘法要解決的問題，實際上就是求以下超定方程組的最小二乘解的問題。其中Ra=y（3）線性最小二乘法的求解定理：當RTR可逆時，超定方程組24線性最小二乘擬合f(x)=a1r1(x)+…+amrm(x)中函數{r1(x),…rm(x)}的選取

1.通過機理分析建立數學模型來確定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx

2.將數據(xi,yi)i=1,…n作圖，通過直觀判斷確定f(x)：線性最小二乘擬合f(x)=a1r1(x)+…+amrm(25用MATLAB解擬合問題1、線性最小二乘擬合2、非線性最小二乘擬合用MATLAB解擬合問題1、線性最小二乘擬合2、非線性最小二26用MATLAB作線性最小二乘擬合1.作多項式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)2.對超定方程組可得最小二乘意義下的解。，用3.多項式在x處的值y可用以下命令計算：

y=polyval（a，x）輸入同長度的數組X，Y擬合多項式次數用MATLAB作線性最小二乘擬合1.作多項式f(x)=a127即要求出二次多項式:中的使得:例對下面一組數據作二次多項式擬合即要求出二次多項式:中的使得:例對下面一組數據作二次281）輸入以下命令：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'解法1．用解超定方程的方法2）計算結果：Ａ=-9.810820.1293-0.03171）輸入以下命令：解法1．用解超定方程的方法2）計算結果：291）輸入以下命令：

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')%作出數據點和擬合曲線的圖形2）計算結果：Ａ=-9.810820.1293-0.0317解法2．用多項式擬合的命令1）輸入以下命令：2）計算結果：Ａ=-9.810301.lsqcurvefit已知數據點：xdata=（xdata1，xdata2，…，xdatan），

ydata=（ydata1，ydata2，…，ydatan）

用MATLAB作非線性最小二乘擬合

Matlab的提供了兩個求非線性最小二乘擬合的函數：lsqcurvefit和lsqnonlin。兩個命令都要先建立M-文件fun.m，在其中定義函數f(x)，但兩者定義f(x)的方式是不同的,可參考例題.

lsqcurvefit用以求含參量x（向量）的向量值函數F(x,xdata)=（F（x，xdata1），…，F（x，xdatan））T中的參變量x(向量),使得

1.lsqcurvefit用MATLAB作非線性最小二乘擬31

輸入格式為:(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,…);fun是一個事先建立的定義函數F(x,xdata)

的M-文件,自變量為x和xdata說明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數據點選項見無約束優化輸入格式為:fun是一個事先建立的定義函數F32

lsqnonlin用以求含參量x（向量）的向量值函數

f(x)=(f1(x),f2(x),…,fn(x))T

中的參量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知數據點：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）lsqnonlin用以求含參量x（向量）的向33輸入格式為：

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）；fun是一個事先建立的定義函數f(x)的M-文件，自變量為x迭代初值選項見無約束優化輸入格式為：說明：x=lsqnonlin（‘fun’，x34

例2用下面一組數據擬合中的參數a，b，k該問題即解最優化問題：例2用下面一組數據擬合35

1）編寫M-文件curvefun1.m

functionf=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.用命令lsqcurvefit1）編寫M-文件curvefun1.m2）輸入命令363）運算結果為：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）結論：a=0.0063,b=-0.0034,k=0.25423）運算結果為：4）結論：a=0.0063,b=-0.0037

解法2

用命令lsqnonlin

f(x)=F(x,tdata,ctada)=x=（a，b，k）1）編寫M-文件curvefun2.m

functionf=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)-cdata2）輸入命令:

x0=[0.2,0.05,0.05];x=lsqnonlin('curvefun2',x0)f=curvefun2(x)解法2用383）運算結果為

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可以看出,兩個命令的計算結果是相同的.4）結論：即擬合得a=0.0063b=-0.0034k=0.25423）運算結果為可以看出,兩個命令的計算結果是相同的.4）結論39MATLAB解應用問題實例1、電阻問題2、給藥方案問題*3、水塔流量估計問題MATLAB解應用問題實例1、電阻問題2、給藥方案問題*3、40電阻問題溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032例.由數據擬合R=a1t+a2方法1.用命令polyfit(x,y,m)得到a1=3.3940,a2=702.4918方法2.直接用結果相同。電阻問題溫度t(0C)20.532.741一室模型：將整個機體看作一個房室，稱中心室，室內血藥濃度是均勻的。快速靜脈注射后，濃度立即上升；然后迅速下降。當濃度太低時，達不到預期的治療效果；當濃度太高，又可能導致藥物中毒或副作用太強。臨床上，每種藥物有一個最小有效濃度c1和一個最大有效濃度c2。設計給藥方案時，要使血藥濃度保持在c1~c2之間。本題設c1=10，c2=25(ug/ml).擬合問題實例2給藥方案——一種新藥用于臨床之前，必須設計給藥方案.藥物進入機體后血液輸送到全身，在這個過程中不斷地被吸收、分布、代謝，最終排出體外，藥物在血液中的濃度，即單位體積血液中的藥物含量，稱為血藥濃度。一室模型：將整個機體看作一個房室，稱中心室，室內血藥濃度是均42

在實驗方面,對某人用快速靜脈注射方式一次注入該藥物300mg后,在一定時刻t(小時)采集血藥,測得血藥濃度c(ug/ml)如下表:

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01要設計給藥方案,必須知道給藥后血藥濃度隨時間變化的規律。從實驗和理論兩方面著手：在實驗方面,對某人用快速靜脈注射方式一次注入該43給藥方案1.在快速靜脈注射的給藥方式下，研究血藥濃度（單位體積血液中的藥物含量）的變化規律。tc2cc10問題2.給定藥物的最小有效濃度和最大治療濃度，設計給藥方案：每次注射劑量多大；間隔時間多長。分析理論：用一室模型研究血藥濃度變化規律實驗：對血藥濃度數據作擬合，符合負指數變化規律給藥方案1.在快速靜脈注射的給藥方式下，研究血藥濃度（單443.血液容積v,t=0注射劑量d,血藥濃度立即為d/v.2.藥物排除速率與血藥濃度成正比，比例系數k(>0)模型假設1.機體看作一個房室，室內血藥濃度均勻——一室模型模型建立

在此，d=300mg，t及c（t）在某些點處的值見前表，需經擬合求出參數k、v3.血液容積v,t=0注射劑量d,血藥濃度立即為d/v.45用線性最小二乘擬合c(t)計算結果：d=300;t=[0.250.511.523468];c=[19.2118.1515.3614.1012.899.327.455.243.01];y=log(c);a=polyfit(t,y,1)k=-a(1)v=d/exp(a(2))程序：用線性最小二乘擬合c(t)計算結果：d=300;程序：46給藥方案設計cc2c10t設每次注射劑量D,間隔時間血藥濃度c(t)

應c1c(t)c2初次劑量D0應加大給藥方案記為：2、1、計算結果：給藥方案：c1=10,c2=25k=0.2347v=15.02給藥方案設計cc2c10t設每次注射劑量D,間隔時間47故可制定給藥方案：即:

首次注射375mg，其余每次注射225mg，注射的間隔時間為4小時。故可制定給藥方案：即:48估計水塔的流量2、解題思路3、算法設計與編程1、問題估計水塔的流量2、解題思路3、算法設計與編程1、問題49某居民區有一供居民用水的園柱形水塔，一般可以通過測量其水位來估計水的流量，但面臨的困難是，當水塔水位下降到設定的最低水位時，水泵自動啟動向水塔供水，到設定的最高水位時停止供水，這段時間無法測量水塔的水位和水泵的供水量．通常水泵每天供水一兩次，每次約兩小時.水塔是一個高12.2米，直徑17.4米的正園柱．按照設計，水塔水位降至約8.2米時，水泵自動啟動，水位升到約10.8米時水泵停止工作．表1是某一天的水位測量記錄，試估計任何時刻（包括水泵正供水時）從水塔流出的水流量，及一天的總用水量．某居民區有一供居民用水的園柱形水塔，一般可以通過測50教學項目4數值分析法模型課件51流量估計的解題思路擬合水位~時間函數確定流量~時間函數估計一天總用水量流量估計的解題思路擬合水位~時間函數確定流量~時間函數估計一52擬合水位~時間函數

測量記錄看，一天有兩個供水時段（以下稱第1供水時段和第2供水時段），和3個水泵不工作時段（以下稱第1時段t=0到t=8.97，第2次時段t=10.95到t=20.84和第3時段t=23以后）．對第1、2時段的測量數據直接分別作多項式擬合，得到水位函數．為使擬合曲線比較光滑，多項式次數不要太高，一般在3~6．由于第3時段只有3個測量記錄，無法對這一時段的水位作出較好的擬合．教學項目4數值分析法模型課件53

2、確定流量~時間函數

對于第1、2時段只需將水位函數求導數即可，對于兩個供水時段的流量，則用供水時段前后（水泵不工作時段）的流量擬合得到，并且將擬合得到的第2供水時段流量外推，將第3時段流量包含在第2供水時段內．2、確定流543、一天總用水量的估計

總用水量等于兩個水泵不工作時段和兩個供水時段用水量之和，它們都可以由流量對時間的積分得到。3、一天總用水量的估計55算法設計與編程1、擬合第1、2時段的水位，并導出流量2、擬合供水時段的流量3、估計一天總用水量4、流量及總用水量的檢驗算法設計與編程1、擬合第1、2時段的水位，并導出流量2、擬合56

1、擬合第1時段的水位，并導出流量設t，h為已輸入的時刻和水位測量記錄（水泵啟動的4個時刻不輸入），第1時段各時刻的流量可如下得：1）c1=polyfit（t（1：10），h（1：10），3）；

%用3次多項式擬合第1時段水位，c1輸出3次多項式的系數2）a1=polyder（c1）；

%a1輸出多項式（系數為c1）導數的系數

3）tp1=0：0.1：9；

x1=-polyval（a1，tp1）；

%x1輸出多項式（系數為a1）在tp1點的函數值（取負后邊為正值），即tp1時刻的流量

4）流量函數為：1、擬合第1時段的水位，并導57

2、擬合第2時段的水位，并導出流量設t，h為已輸入的時刻和水位測量記錄（水泵啟動的4個時刻不輸入），第2時段各時刻的流量可如下得：1）c2=polyfit(t(10.9:21),h(10.9:21),3)；

%用3次多項式擬合第2時段水位，c2輸出3次多項式的系數2）a2=polyder(c2)；

%a2輸出多項式（系數為c2）導數的系數

3）tp2=10.9:0.1:21;x2=-polyval(a2,tp2);%x2輸出多項式（系數為a2）在tp2點的函數值（取負后邊為正值），即tp2時刻的流量4）流量函數為：2、擬合第2時段的水位，并導58

3、擬合供水時段的流量在第1供水時段（t=9~11）之前（即第1時段）和之后（即第2時段）各取幾點，其流量已經得到，用它們擬合第1供水時段的流量．為使流量函數在t=9和t=11連續，我們簡單地只取4個點，擬合3次多項式（即曲線必過這4個點），實現如下：

xx1=-polyval（a1，[89]）；%取第1時段在t=8,9的流量

xx2=-polyval（a2，[1112]）；%取第2時段在t=11,12的流量

xx12=[xx1xx2]；

c12=polyfit（[891112]，xx12，3）；%擬合3次多項式

tp12=9：0.1：11；

x12=polyval（c12，tp12）；%x12輸出第1供水時段各時刻的流量擬合的流量函數為：3、擬合供水時段的流量擬合的59在第2供水時段之前取t=20，20.8兩點的流水量，在該時刻之后（第3時段）僅有3個水位記錄，我們用差分得到流量，然后用這4個數值擬合第2供水時段的流量如下：

dt3=diff（t(22：24)）；%最后3個時刻的兩兩之差

dh3=diff（h(22：24)）；%最后3個水位的兩兩之差

dht3=-dh3./dt3；%t(22)和t(23)的流量

t3=[2020.8t(22)t(23)]；

xx3=[-polyval(a2，t3(1：2)，dht3)]；%取t3各時刻的流量

c3=polyfit（t3，xx3，3）；%擬合3次多項式

t3=20.8：0.1：24；

x3=polyval（c3，tp3）；%x3輸出第2供水時段（外推至t=24）各時刻的流量擬合的流量函數為：在第2供水時段之前取t=20，20.8兩點的流60

3、一天總用水量的估計第1、2時段和第1、2供水時段流量的積分之和，就是一天總用水量．雖然諸時段的流量已表為多項式函數，積分可以解析地算出，這里仍用數值積分計算如下：

y1=0.1*trapz(x1)；%第1時段用水量（仍按高度計），0.1為積分步長

y2=0.1*trapz(x2)；%第2時段用水量

y12=0.1*trapz(x12)；%第1供水時段用水量

y3=0.1*trapz(x3)；%第2供水時段用水量

y=（y1+y2+y12+y3）*237.8*0.01；%一天總用水量（）計算結果：y1=146.2,y2=266.8,y12=47.4,y3=77.3，y=1250.43、一天總用水量的估計61

4、流量及總用水量的檢驗計算出的各時刻的流量可用水位記錄的數值微分來檢驗．用水量y1可用第1時段水位測量記錄中下降高度968-822=146來檢驗，類似地，y2用1082-822=260檢驗．供水時段流量的一種檢驗方法如下：供水時段的用水量加上水位上升值260是該時段泵入的水量，除以時段長度得到水泵的功率（單位時間泵入的水量），而兩個供水時段水泵的功率應大致相等．第1、2時段水泵的功率可計算如下：

p1=(y12+260)/2；%第1供水時段水泵的功率（水量仍以高度計）

tp4=20.8：0.1：23；

xp2=polyval（c3，tp4）；%xp2輸出第2供水時段各時刻的流量

p2=(0.1*trapz(xp2)+260)/2.2；%第2供水時段水泵的功率（水量仍以高度計）計算結果：p1=154.5，p2=140.14、流量及總用水量的檢驗62計算結果流量函數為：計算結果流量函數為：63流量曲線見圖n=(3,4)n=(5,6)流量曲線見圖n=(3,4)n=(5,6)64數學建模與實驗嚴可頌主講數學建模與實驗嚴可頌主講65教學項目四數值分析法模型教學項目四數值分析法模型664.1插值法建模

拉格朗日插值分段線性插值三次樣條插值一、插值的定義二、插值的方法三、用Matlab解插值問題4.1插值法建模拉格朗日插值分段線性插值三次樣條插值一、67已知n+1個節點其中互不相同，不妨設求任一插值點處的插值節點可視為由產生，表達式復雜，或無封閉形式，或未知。4.1.1插值法的定義

已知n+1個節點其中互不相同，不妨設求任一插值點處的插值68構造一個(相對簡單的)函數通過全部節點,即再用計算插值，即構造一個(相對簡單的)函數通過全部節點,即再用計算插值，694.1.2插值的方法

已知函數f(x)在n+1個點x0,x1,…,xn處的函數值為y0,y1,…,yn

。求一n次多項式函數Pn(x)，使其滿足：拉格朗日(Lagrange)插值4.1.2插值的方法已知函數f(x)在n+170稱為拉格朗日插值基函數。解決此問題的拉格朗日插值多項式公式如下其中Li(x)為n次多項式：稱為拉格朗日插值基函數。解決此問題的拉格朗日插值多項式公71特殊情況兩點一次(線性)插值多項式:三點二次(拋物)插值多項式:特殊情況兩點一次(線性)插值多項式:三點二次(拋物)插值多項72采用拉格朗日多項式插值：選取不同插值節點個數n+1，其中n為插值多項式的次數，當n取10時，繪出插值結果圖形.例采用拉格朗日多項式插值：選取不同插值節點個數73分段線性插值計算量與n無關;n越大，誤差越小.xjxj-1xj+1x0xnxoy分段線性插值計算量與n無關;xjxj-1xj+174比分段線性插值更光滑。xyxi-1xiab在數學上，光滑程度的定量描述是：函數(曲線)的k階導數存在且連續，則稱該曲線具有k階光滑性。光滑性的階次越高，則越光滑。是否存在較低次的分段多項式達到較高階光滑性的方法？三次樣條插值就是一個很好的例子。三次樣條插值比分段線性插值更光滑。xyxi-175

三次樣條插值g(x)為被插值函數。三次樣條插值g(x)為被插值函數。76用MATLAB作插值計算一維插值函數：yi=interp1(x，y，xi，'method')插值方法被插值點插值節點xi處的插值結果‘nearest’

：最鄰近插值‘linear’

：線性插值；‘spline’

：三次樣條插值；‘cubic’

：立方插值。缺省時：分段線性插值。注意：所有的插值方法都要求x是單調的，并且xi不能夠超過x的范圍。用MATLAB作插值計算一維插值函數：yi=interp1(77例：在1-12的11小時內，每隔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');(直接輸出數據將是很多的)plot(hours,temps,'+',h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius’)例：在1-12的11小時內，每隔1小時測量一次溫度，78xy實驗習題

已知飛機下輪廓線上數據如下，求x每改變0.1時的y值。機翼下輪廓線xy實驗習題機翼下輪廓線794.2擬合法簡介

2.擬合的基本原理1.擬合問題引例4.2擬合法簡介2.擬合的基本原理1.擬合問題引例80擬合問題引例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數據：求600C時的電阻R。

設

R=at+ba,b為待定系數擬合問題引例1溫度t(0C)20.581擬合問題引例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速靜脈注射下的血藥濃度數據(t=0注射300mg)求血藥濃度隨時間的變化規律c(t).作半對數坐標系(semilogy)下的圖形擬合問題引例2t(h)82曲線擬合問題的提法已知一組（二維）數據，即平面上n個點（xi,yi)i=1,…n,

尋求一個函數（曲線）y=f(x),

使f(x)

在某種準則下與所有數據點最為接近，即曲線擬合得最好。

+++++++++xyy=f(x)(xi,yi)ii為點（xi,yi)與曲線y=f(x)的距離曲線擬合問題的提法已知一組（二維）數據，即平83擬合與插值的關系

函數插值與曲線擬合都是要根據一組數據構造一個函數作為近似，由于近似的要求不同，二者的數學方法上是完全不同的。實例：下面數據是某次實驗所得，希望得到X和f之間的關系？問題：給定一批數據點，需確定滿足特定要求的曲線或曲面解決方案：若不要求曲線（面）通過所有數據點，而是要求它反映對象整體的變化趨勢，這就是數據擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）通過所給所有數據點，就是插值問題；擬合與插值的關系函數插值與曲線擬合都是要根據一組84最臨近插值、線性插值、樣條插值與曲線擬合結果：最臨近插值、線性插值、樣條插值與曲線擬合結果：85曲線擬合問題最常用的解法——線性最小二乘法的基本思路第一步:先選定一組函數

r1(x),r2(x),…rm(x),m<n,

令

f(x)=a1r1(x)+a2r2(x)+…+amrm(x)（1）其中

a1,a2,…am

為待定系數。

第二步:確定a1,a2,…am

的準則（最小二乘準則）：使n個點（xi,yi)與曲線y=f(x)的距離i的平方和最小

。記

問題歸結為，求

a1,a2,…am

使

J(a1,a2,…am)

最小。曲線擬合問題最常用的解法——線性最小二乘法的基本思路第一步:86線性最小二乘法的求解：預備知識超定方程組：方程個數大于未知量個數的方程組即Ra=y其中超定方程一般是不存在解的矛盾方程組。如果有向量a使得達到最小，則稱a為上述超定方程的最小二乘解。線性最小二乘法的求解：預備知識超定方程組：方程個數大于未知量87線性最小二乘法的求解

定理：當RTR可逆時，超定方程組（3）存在最小二乘解，且即為方程組

RTRa=RTy的解：a=(RTR)-1RTy

所以，曲線擬合的最小二乘法要解決的問題，實際上就是求以下超定方程組的最小二乘解的問題。其中Ra=y（3）線性最小二乘法的求解定理：當RTR可逆時，超定方程組88線性最小二乘擬合f(x)=a1r1(x)+…+amrm(x)中函數{r1(x),…rm(x)}的選取

1.通過機理分析建立數學模型來確定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx

2.將數據(xi,yi)i=1,…n作圖，通過直觀判斷確定f(x)：線性最小二乘擬合f(x)=a1r1(x)+…+amrm(89用MATLAB解擬合問題1、線性最小二乘擬合2、非線性最小二乘擬合用MATLAB解擬合問題1、線性最小二乘擬合2、非線性最小二90用MATLAB作線性最小二乘擬合1.作多項式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)2.對超定方程組可得最小二乘意義下的解。，用3.多項式在x處的值y可用以下命令計算：

y=polyval（a，x）輸入同長度的數組X，Y擬合多項式次數用MATLAB作線性最小二乘擬合1.作多項式f(x)=a191即要求出二次多項式:中的使得:例對下面一組數據作二次多項式擬合即要求出二次多項式:中的使得:例對下面一組數據作二次921）輸入以下命令：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'解法1．用解超定方程的方法2）計算結果：Ａ=-9.810820.1293-0.03171）輸入以下命令：解法1．用解超定方程的方法2）計算結果：931）輸入以下命令：

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')%作出數據點和擬合曲線的圖形2）計算結果：Ａ=-9.810820.1293-0.0317解法2．用多項式擬合的命令1）輸入以下命令：2）計算結果：Ａ=-9.810941.lsqcurvefit已知數據點：xdata=（xdata1，xdata2，…，xdatan），

ydata=（ydata1，ydata2，…，ydatan）

用MATLAB作非線性最小二乘擬合

Matlab的提供了兩個求非線性最小二乘擬合的函數：lsqcurvefit和lsqnonlin。兩個命令都要先建立M-文件fun.m，在其中定義函數f(x)，但兩者定義f(x)的方式是不同的,可參考例題.

lsqcurvefit用以求含參量x（向量）的向量值函數F(x,xdata)=（F（x，xdata1），…，F（x，xdatan））T中的參變量x(向量),使得

1.lsqcurvefit用MATLAB作非線性最小二乘擬95

輸入格式為:(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,…);fun是一個事先建立的定義函數F(x,xdata)

的M-文件,自變量為x和xdata說明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數據點選項見無約束優化輸入格式為:fun是一個事先建立的定義函數F96

lsqnonlin用以求含參量x（向量）的向量值函數

f(x)=(f1(x),f2(x),…,fn(x))T

中的參量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知數據點：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）lsqnonlin用以求含參量x（向量）的向97輸入格式為：

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）；fun是一個事先建立的定義函數f(x)的M-文件，自變量為x迭代初值選項見無約束優化輸入格式為：說明：x=lsqnonlin（‘fun’，x98

例2用下面一組數據擬合中的參數a，b，k該問題即解最優化問題：例2用下面一組數據擬合99

1）編寫M-文件curvefun1.m

functionf=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.用命令lsqcurvefit1）編寫M-文件curvefun1.m2）輸入命令1003）運算結果為：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）結論：a=0.0063,b=-0.0034,k=0.25423）運算結果為：4）結論：a=0.0063,b=-0.00101

解法2

用命令lsqnonlin

f(x)=F(x,tdata,ctada)=x=（a，b，k）1）編寫M-文件curvefun2.m

functionf=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)-cdata2）輸入命令:

x0=[0.2,0.05,0.05];x=lsqnonlin('curvefun2',x0)f=curvefun2(x)解法2用1023）運算結果為

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可以看出,兩個命令的計算結果是相同的.4）結論：即擬合得a=0.0063b=-0.0034k=0.25423）運算結果為可以看出,兩個命令的計算結果是相同的.4）結論103MATLAB解應用問題實例1、電阻問題2、給藥方案問題*3、水塔流量估計問題MATLAB解應用問題實例1、電阻問題2、給藥方案問題*3、104電阻問題溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032例.由數據擬合R=a1t+a2方法1.用命令polyfit(x,y,m)得到a1=3.3940,a2=702.4918方法2.直接用結果相同。電阻問題溫度t(0C)20.532.7105一室模型：將整個機體看作一個房室，稱中心室，室內血藥濃度是均勻的。快速靜脈注射后，濃度立即上升；然后迅速下降。當濃度太低時，達不到預期的治療效果；當濃度太高，又可能導致藥物中毒或副作用太強。臨床上，每種藥物有一個最小有效濃度c1和一個最大有效濃度c2。設計給藥方案時，要使血藥濃度保持在c1~c2之間。本題設c1=10，c2=25(ug/ml).擬合問題實例2給藥方案——一種新藥用于臨床之前，必須設計給藥方案.藥物進入機體后血液輸送到全身，在這個過程中不斷地被吸收、分布、代謝，最終排出體外，藥物在血液中的濃度，即單位體積血液中的藥物含量，稱為血藥濃度。一室模型：將整個機體看作一個房室，稱中心室，室內血藥濃度是均106

在實驗方面,對某人用快速靜脈注射方式一次注入該藥物300mg后,在一定時刻t(小時)采集血藥,測得血藥濃度c(ug/ml)如下表:

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01要設計給藥方案,必須知道給藥后血藥濃度隨時間變化的規律。從實驗和理論兩方面著手：在實驗方面,對某人用快速靜脈注射方式一次注入該107給藥方案1.在快速靜脈注射的給藥方式下，研究血藥濃度（單位體積血液中的藥物含量）的變化規律。tc2cc10問題2.給定藥物的最小有效濃度和最大治療濃度，設計給藥方案：每次注射劑量多大；間隔時間多長。分析理論：用一室模型研究血藥濃度變化規律實驗：對血藥濃度數據作擬合，符合負指數變化規律給藥方案1.在快速靜脈注射的給藥方式下，研究血藥濃度（單1083.血液容積v,t=0注射劑量d,血藥濃度立即為d/v.2.藥物排除速率與血藥濃度成正比，比例系數k(>0)模型假設1.機體看作一個房室，室內血藥濃度均勻——一室模型模型建立

在此，d=300mg，t及c（t）在某些點處的值見前表，需經擬合求出參數k、v3.血液容積v,t=0注射劑量d,血藥濃度立即為d/v.109用線性最小二乘擬合c(t)計算結果：d=300;t=[0.250.511.523468];c=[19.2118.1515.3614.1012.899.327.455.243.01];y=log(c);a=polyfit(t,y,1)k=-a(1)v=d/exp(a(2))程序：用線性最小二乘擬合c(t)計算結果：d=300;程序：110給藥方案設計cc2c10t設每次注射劑量D,間隔時間血藥濃度c(t)

應c1c(t)c2初次劑量D0應加大給藥方案記為：2、1、計算結果：給藥方案：c1=10,c2=25k=0.2347v=15.02給藥方案設計cc2c10t設每次注射劑量D,間隔時間111故可制定給藥方案：即:

首次注射375mg，其余每次注射225mg，注射的間隔時間為4小時。故可制定給藥方案：即:112估計水塔的流量2、解題思路3、算法設計與編程1、問題估計水塔的流量2、解題思路3、算法設計與編程1、問題113某居民區有一供居民用水的園柱形水塔，一般可以通過測量其水位來估計水的流量，但面臨的困難是，當水塔水位下降到設定的最低水位時，水泵自動啟動向水塔供水，到設定的最高水位時停止供水，這段時間無法測量水塔的水位和水泵的供水量．通常水泵每天供水一兩次，每次約兩小時.水塔是一個高12.2米，直徑17.4米的正園柱．按照設計，水塔水位降至約8.2米時，水泵自動啟動，水位升到約10.8米時水泵停止工作．表1是某一天的水位測量記錄，試估計任何時刻（包括水泵正供水時）從水塔流出的水流量，及一天的總用水量．某居民區有一供居民用水的園柱形水塔，一般可以通過測114教學項目4數值分析法模型課件115流量估計的解題思路擬合水位~時間函數確定流量~時間函數估計一天總用水量流量估計的解題思路擬合水位~時間函數確定流量~時間函數估計一116擬合水位~時間函數

測量記錄看，一天有兩個供水時段（以下稱第1供水時段和第2供水時段），和3個水泵不工作時段（以下稱第1時段t=0到t=8.97，第2次時段t=10.95到t=20.84和第3時段t=23以后）．對第1、2時段的測量數據直接分別作多項式擬合，得到水位函數．為使擬合曲線比較光滑，多項式次數不要太高，一般在3~6．由于第3時段只有3個測量記錄，無法對這一時段的水位作出較好的擬合．教學項目4數值分析法模型課件117

2、確定流量~時間函數

對于第1、2時段只需將水位函數求導數即可，對于兩個供水時段的流量，則用供水時段前后（水泵不工作時段）的流量擬合得到，并且將擬合得到的第2供水時段流量外推，將第3時段流量包含在第2供水時段內．2、確定流1183、一天總用水量的估計

總用水量等于兩個水泵不工作時段和兩個供水時段用水量之和，它們都可以由流量對時間的積分得到。3、一天總用水量的估計119算法設計與編程1、擬合第1、2時段的水位，并導出流量2、擬合供水時段的流量3、估計一天總用水量4、流量及總用水量的檢驗算法設計與編程1、擬合第1、2時段的水位，并導出流量2、擬合120

1、擬合第1時段的水位，并導出流量設t，h為已輸入的時刻和水位測量記錄（水泵啟動的4個時刻不輸入），第1時段各時刻的流量可如下得：1）c1=polyfit（t（1：10），h（1：10），3）；

%用3次多項式擬合第1時段水位，c1輸出3次多項式的系數2）a1=polyder（c1）；

%a1輸出多項式（系數為c1）導數的系數

3）tp1=0：0.1：9；

x1=-polyval（a1，tp1）；

%x1輸出多項式（系數為a1）在tp1點的函數值（取負后邊為正值），即tp1時刻的流量

4）流量函數為：1、擬合第1時段的水位，并導121

2、擬合第2時段的水位，并導出流量設t