計算機應(yīng)用基礎(chǔ)積分方程及應(yīng)用

計算機應(yīng)用基礎(chǔ)積分方程及應(yīng)用第1頁，共40頁，2023年，2月20日，星期四一、插值基礎(chǔ)

設(shè)函數(shù)y=f(x)在區(qū)間[a,b]上有定義，且已知在點a≤x0<x1<…<xn≤b上的值y0,y1,…,yn,若存在一簡單函數(shù)P(x),使

P（xi）=yi(i=0,1,2,…,n)（4-1）

成立，就稱P(x)為f(x)的插值函數(shù)，點x0，x1，…，xn稱為插值節(jié)點，包含插值節(jié)點的區(qū)間[a,b]稱為插值區(qū)間，求插值函數(shù)P(x)的方法稱為插值法。

若P(x)為插值多項式，即

P(x)=a0+a1x+…+anxn

（4-2）

其中ai為實數(shù)，就稱P(x)為插值多項式，相應(yīng)的插值法稱為多項式插值。若P(x)為分段多項式，就稱分段插值。

4.1插值函數(shù)第2頁，共40頁，2023年，2月20日，星期四

從幾何圖形上看，插值法就是求曲線y=P(x),使其通過給定的n+1個點（xi,yi),i=0,1,…,n,并用它近似已知曲線y=f(x)。

插值多項式存在的唯一性

設(shè)P(x)是形如（4-2）的插值多項式，則滿足條件（4-1）的插值多項式（4-2）存在唯一。4.1插值函數(shù)第3頁，共40頁，2023年，2月20日，星期四一線性插值與拋物插值

通過平面的兩點可以確定一條直線經(jīng)過兩點，這就是拉格朗日線性插值；對于不在一條直線的3個點得到的插值多項式稱為拋物線插值。二、拉格朗日插值多項式插值多項式將上式展開，就得到插值的基函數(shù)，插值多項式為4.1插值函數(shù)第4頁，共40頁，2023年，2月20日，星期四拉格朗日插值函數(shù)通用程序functionf=Lang(x,y,x0)symst;if(length(x)==length(y))n=length(x);elsedisp('x和y的維數(shù)不相等！');return;end%檢錯f=0;for(i=1:n)l=y(i);for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end;%x為已知數(shù)據(jù)點的x坐標向量；%y為已知數(shù)據(jù)點的y坐標向量；%x0為插值點的x坐標；%f為求得的拉格朗日插值多項式；%f0為在x0處的插值4.1插值函數(shù)第5頁，共40頁，2023年，2月20日，星期四for(j=i+1:n)l=l*(t-x(j))/(x(i)-x(j));%計算拉格朗日基函數(shù)

end;f=f+l;%計算拉格朗日插值函數(shù)

simplify(f);%化簡

if(i==n)if(nargin==3)f=subs(f,'t',x0);%計算插值點的函數(shù)值

elsef=collect(f);%將插值多項式展開

f=vpa(f,6);%將插值多項式的系數(shù)化成6位精度的小數(shù)

endendend4.1插值函數(shù)第6頁，共40頁，2023年，2月20日，星期四三Matlab自帶函數(shù)插值

3.1一維插值

Matlab提供了interp1函數(shù)用于一維插值，其調(diào)用格式為yi=interp1(x,Y,xi,method)對節(jié)點（x,Y)進行插值，計算插值點xi的函數(shù)值；Method插值算法，默認的為線性插值：

nearest:線性最近插值

linear:線性插值（默認）

spline:三次樣條插值

pchip:分段三次埃爾米特插值

cubic:雙三次插值4.1插值函數(shù)第7頁，共40頁，2023年，2月20日，星期四例【4-1】已知函數(shù)f(x)=(x2-12x+9)e-4xcosx生成數(shù)據(jù)點，試根據(jù)生成的數(shù)據(jù)進行插值處理。比較不同算法，計算程序“ME_7_1.m”4.1插值函數(shù)例【4-2】：某觀測站測得某日6:00時至18:00時之間每隔2小時的室內(nèi)外溫度(℃)，用3次樣條插值分別求得該日室內(nèi)外6:30至17:30時之間每隔2小時各點的近似溫度(℃)。第8頁，共40頁，2023年，2月20日，星期四設(shè)時間變量h為一行向量，溫度變量t為一個兩列矩陣，其中第一列存放室內(nèi)溫度，第二列儲存室外溫度。命令如下：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次樣條插值計算4.1插值函數(shù)第9頁，共40頁，2023年，2月20日，星期四3.2二維插值

Matlab提供interp2用于二維插值，適用于二維網(wǎng)格數(shù)據(jù)插值，其調(diào)用格式如下：ZI=interp2(X,Y,Z,XI,YI,method)X,Y均為矩陣在2-D區(qū)域的點，這些點的數(shù)值Z已知，因此構(gòu)造插值函數(shù)Z=F(X,Y)4.1插值函數(shù)第10頁，共40頁，2023年，2月20日，星期四3.3二維隨機點的插值對于任意多組數(shù)據(jù)（xi,yi,zi)采用griddata函數(shù)，其調(diào)用格式：[XI,YI,ZI]=griddata(x,y,z,XI,YI,method)x,y,z為已知樣本點，[XI,YI]插值點的位置【例4-3】在區(qū)間【-2,2】上創(chuàng)建100個隨機點，用griddata函數(shù)對這些數(shù)據(jù)點進行插值計算程序“ME_7_2.m”4.1插值函數(shù)第11頁，共40頁，2023年，2月20日，星期四4.1插值函數(shù)3.4三次樣條插值對于插值基點很多時，采用多項式插值效果并不是很好，采用三次樣條插值效果更好一些，其調(diào)用格式：yy=spline(x,Y,xx),計算三次樣條插值的分段多項式，采用ppval(pp,x)計算多項式在x點的值。pp=spline(x,Y)例4-4利用spline函數(shù)對y=tan(pi*x/24)進行樣條插值。ME_7_6.m第12頁，共40頁，2023年，2月20日，星期四functionME_7_6clear;x=0:10;y=tan(pi*x/24);xi=linspace(0,10);yi=spline(x,y,xi);plot(x,y,'rp',xi,yi);title('插值效果圖');4.1插值函數(shù)第13頁，共40頁，2023年，2月20日，星期四3.5B樣條插值是一類常用的樣條函數(shù)，其調(diào)用格式為：Spline=spapi(knots,x,y)spapi(k,x,y),k為B樣條階數(shù)4.1插值函數(shù)第14頁，共40頁，2023年，2月20日，星期四一最小二乘曲線擬合4.2擬合函數(shù)有一組離散數(shù)據(jù)xi,yi,i=1,2…N,滿足某一函數(shù)的原型yˇ(x)=f(a,x),其中a為待定系數(shù)向量，則最小二乘曲線擬合的目標：求出這一組待定系數(shù)的值，使得目標函數(shù)為最小。Matlab調(diào)用格式如下：第15頁，共40頁，2023年，2月20日，星期四4.2擬合函數(shù)其中Fun為原型函數(shù)的Matlab表示，可以使M函數(shù)inline()函數(shù)；a0

為最優(yōu)化的初值；

x,y為原始輸入和輸出數(shù)據(jù)向量；

a為返回的待定系數(shù)向量；

Jm為在此待定系數(shù)下目標函數(shù)的值。第16頁，共40頁，2023年，2月20日，星期四4.2擬合函數(shù)

二最小二乘多項式擬合對于離散型函數(shù)，若數(shù)據(jù)點較多，若將每個數(shù)據(jù)點都當做插值節(jié)點，運算顯得非常復(fù)雜。在工程試驗中，常測得一組離散數(shù)據(jù)點（xi,yi),(i=1,2…N),要求y=(x），這種應(yīng)變量只有一個自變量的數(shù)據(jù)擬合方法稱之為直線擬合。(仍然采用最小二乘方法)p=polyfit(x,y,n)其中x和y為原始的樣本點構(gòu)成的向量n為選定的多項式階數(shù)p為多項式系數(shù)按降冪排列得出的行向量第17頁，共40頁，2023年，2月20日，星期四4.2擬合函數(shù)函數(shù)polyval()

根據(jù)由polyfit（）計算得到的多項式系數(shù)向量，計算在x處多項式函數(shù)的值y，它通常與polyfit（）聯(lián)合使用，調(diào)用格式：y=polyval(p,x)p:多項式函數(shù)；x:需要計算點的x值。第18頁，共40頁，2023年，2月20日，星期四4.2擬合函數(shù)三兩種擬合方法的比較最小二乘曲線擬合最小二乘多項式擬合p=polyfit(x,y,n)擬合的函數(shù)是固定的,且是降階的多項式,因此擬合調(diào)用不需要注明擬合函數(shù),不需要初值。相同之處:都是采用最小二乘法優(yōu)化獲得擬合曲線系數(shù)擬合的函數(shù)形式可任意,因此擬合調(diào)用需要注明擬合函數(shù),即需要建立一個Fun的函數(shù)，需要初值，而且結(jié)果與初值密切相關(guān)。第19頁，共40頁，2023年，2月20日，星期四【例4-5】進行多項式擬合x0=0:0.1:1;y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);p3=polyfit(x0,y0,3)vpa(poly2sym(p3),10)x=0:0.01:1;%繪圖ya=(x.^2-3*x+5).*exp(-5*x).*sin(x);y1=polyval(p3,x);plot(x,y1,x,ya,x0,y0,‘o’)擬合程序“ME_7_3”結(jié)果如下：4.2擬合函數(shù)

已知的數(shù)據(jù)點來自f(x)=(x2-3x+5)e-5xsinx第20頁，共40頁，2023年，2月20日，星期四p3=2.8400-4.78981.94320.0598f(x)=2.84x3-4.79x2+1.94x+0.05984.2擬合函數(shù)第21頁，共40頁，2023年，2月20日，星期四例【4-6】4.2擬合函數(shù)已知數(shù)據(jù)可能滿足y(x)=ax+bx2e-cx+dXi0.10.20.30.40.5Yi2.322.652.973.293.6Xi0.60.70.80.91Yi3.914.214.524.825.13求滿足最小二乘解a,b,c,d的值令a1=a,a2=b,a3=c,a4=d,則原型函數(shù)可寫成：

y(x)=a1x+a2x2e-a3x+a4第22頁，共40頁，2023年，2月20日，星期四functionME_7_4clear;x=0.1:0.1:1;y=[2.32,2.65,2.97,3.29,3.6,...3.91,4.21,4.52,4.82,5.13];a=lsqcurvefit(@c8f3,[1;2;2;3],x,y);a'y1=c8f3(a,x);plot(x,y,x,y1)計算程序4.2擬合函數(shù)functiony=c8f3(a,x)y=a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4);編寫原型函數(shù)程序第23頁，共40頁，2023年，2月20日，星期四4.2擬合函數(shù)四最小二乘擬合生成樣條曲線函數(shù)csaps()的用法:

平滑生成三次樣條曲線,調(diào)用格式:樣條擬合曲線不要求曲線通過全部實驗數(shù)據(jù)點，更適用實驗觀察得到的離散數(shù)據(jù)點sp=csaps(x,y,p)ys=casps(x,y,p,xx,w)函數(shù)span2()的用法:

用最小二乘法生成B樣條曲線,調(diào)用格式:sp=span2(knots,k,x,y)sp=span2(knots,k,x,y,w)第24頁，共40頁，2023年，2月20日，星期四4.2擬合函數(shù)函數(shù)spaps()的用法:

平滑生成B樣條曲線,調(diào)用格式:sp=spaps(x,y,tol)[sp,ys]=spaps(x,y,tol,m,w)函數(shù)fnval()的用法:計算函數(shù)f在給定點x處的函數(shù)值，調(diào)用格式:values=finval(f,x)第25頁，共40頁，2023年，2月20日，星期四t/h青霉素濃度t/h青霉素濃度001209430201061401095040160016010280603000180962080581020094001008600例4-6青霉素發(fā)酵實驗數(shù)據(jù)如表所示，請分別用插值方法和最小二乘法估算t=10,30,50,70,90,110,150,190h時青霉素的濃度。并分別畫出插值函數(shù)和擬合函數(shù)曲線，并標注插值點。（ME_7_5.m)4.2擬合函數(shù)第26頁，共40頁，2023年，2月20日，星期四4.3數(shù)值積分及應(yīng)用一復(fù)合辛普森法積分二梯形求積三單變量數(shù)值積分問題求解四雙重積分問題的數(shù)值解第27頁，共40頁，2023年，2月20日，星期四4.3數(shù)值積分問題對于一個數(shù)學積分關(guān)系式：

常用的解法有：矩形法、梯形法和辛普森法。

一復(fù)合辛普森法

基本思想：將函數(shù)f(x)的積分區(qū)間n等分（n為偶數(shù)），以兩個小區(qū)間作為一個計算單元，用一個簡單的二次函數(shù)近似被積函數(shù)，則函數(shù)f(x)在區(qū)間[a,b]積分可表述為：

S=h/3{f(a)-f(b)+2∑[2f(a+(2i-1)h)+f(a+2ih)]}

第28頁，共40頁，2023年，2月20日，星期四離散誤差：4.3數(shù)值積分問題調(diào)用格式：q=quad(@fun,a,b,tol)fun:被積函數(shù)，必須是向量表達式，表達式必

須用點運算符。a,b：積分限tol：絕對誤差限函數(shù)quad1()精度更高，使用方法同前。第29頁，共40頁，2023年，2月20日，星期四例4-4用復(fù)合simposn法計算以下積分，要求精度在1e-5內(nèi)，并計算誤差。第30頁，共40頁，2023年，2月20日，星期四要使則M>=5,取m=5,h=0.2程序ME_3_2.m,S=3.8875,err=2.72864.3數(shù)值積分問題第31頁，共40頁，2023年，2月20日，星期四二由給定數(shù)據(jù)進行梯形求積4.3數(shù)值積分問題第32頁，共40頁，2023年，2月20日，星期四4.3數(shù)值積分問題利用trapz()函數(shù)進行積分求解,格式如下:S=trapz(x,y)例4-7用定步長方法求解積分畫圖：>>x=[0:0.01:3*pi/2,3*pi/2];y=cos(15*x);plot(x,y)計算程序：ME_4_9第33頁，共40頁，2023年，2月20日，星期四求理論值：不同步距：>>symsx,A=int(cos(15*x),0,3*pi/2)h0=[0.1,0.01,0.001,0.0001,0.00001,0.000001];v=[];Forh=h0x=[0:h:3*pi/2,3*pi/2];y=cos(15*x);I=trapz(x,y);v=[v;h,I,1/15-I];end4.3數(shù)值積分問題

>>symsx,A=int(cos(15*x),0,3*pi/2)第34頁，共40頁，2023年，2月20日，星期四二雙重積分問題的數(shù)值解4.3數(shù)值積分問題求解格式:(1)y=dblquad(Fun,xm,xM,ym,yM)%矩形區(qū)域的雙重積分(2)y=dblquad(Fun,xm,xM,ym,yM,)%限定精度的雙重積分

第3