數值積分的matlab實現

(10.2)(10.2)....3．會用數值積分方法解決一些實際問題。以仍然得不到積分的精確值，如不定積分j1sinxdx。這時我們一般考慮用數值方法計算其0x近似值，稱為數值積分。設函數y=f(x)在x*可導，則其導數為f,(x*)=limf(x*+h)f(x*)h0h .1)如果函數y=f(x)以列表形式給出(見表10-1)，則其精確值無法求得，但可由下式求得其近似值f,(x*)必f(x*+h)f(x*)h表10-1xyxxxn2xyxxxn2yyy2n2一般的，步長h越小，所得結果越精確。(10.2)式右端項的分子稱為函數y=f(x)在x*的差分，分母稱為自變量在x*的差分，所以右端項又稱為差商。數值微分即用差商近似用的差商公式為：2hf(x)必2h 0.4)......f(x)必n2h .5)其誤差均為O(h2)，稱為統(tǒng)稱三點公式。MATLAB提供了一個指令求解一階向前差分，其使用格式為：dx=diff(x)2132n1m例1用三點公式計算y=f(x)在x=1.0,1.2,1.4處的導數值，f(x)的值由下表給xf(x)functionf=diff3(x,y)n=length(x);h=x(2)-x(1);f(1)=(-3*y(1)+4*y(2)-y(3))/(2*h);forj=2:n-1f(j)=(y(j+1)-y(j-1))/(2*h);f(n)=(y(n-2)-4*y(n-1)+3*y(n))/(2*h);x=[1.0,1.1,1.2,1.3,1.4];y=[0.2500,0.2268,0.2066,0.1890,0.1736];diff3(x,y)step1：對給定數據點(x,y)，利用指令pp=spline(x,y)，獲得三次樣條函數數據pp，的階數、段數、節(jié)點的橫坐標值和各段多項式的系數。step2：對于上面所求的數據向量pp，利用指令[breaks,coefs,m,n]=unmkpp(pp)進行段多項式pp。functiondy=ppd(pp)[breaks,coefs,m]=unmkpp(pp);........fori=1:mcoefsm(i,:)=polyder(coefs(i,:));dy=mkpp(breaks,coefsm);xyxx處的導數dyy：pp=spline(x,y),dy=ppd(pp);dyy=ppval(dy,xx);例2基于正弦函數ysinx的數據點，利用三點公式和三次樣條插值分別求導，并與解析所求得的導數進行比較。h=0.1*pi;x=0:h:2*pi;y=sin(x);dy1=diff3(x,y);pp=spline(x,y);dy=ppd(pp);dy2=ppval(dy,x);z=cos(x);error1=norm(dy1-z),error2=norm(dy2-z)plot(x,dy1,'k:',x,dy2,'r--',x,z,'b')ror顯然利用三次樣條插值求導所得誤差比三點公式求導小很多，同時由圖2.15可知利用三次樣條插值求導所得曲線與解析求導曲線基本重合，而三點公式在極值點附近和兩個端點附近誤差較大，其它點吻合的較好。問題：湖水在夏天會出現分層現象，其特點是接近湖面的水的溫度較高，越往下水的溫度越低。這種現象會影響水的對流和混合過程，使得下層水域缺氧，導致水生魚類死亡。個湖的水溫進行觀測得數據見表10-2。深度(m)溫度(℃)0試找出湖水溫度變化最大的深度。湖水的溫度可視為關于深度的函數，于是湖水溫度的變化問題便轉化為溫度函數的導........的數據較少，由此計算的深度不夠精確，所以采用插值的方法計算加密深度數據的導數值，2．模型的建立及求解記湖水的深度為h(m)，相應的溫度為T(℃)，且有T=T(h)，并假定函數T(h)可對給定的數據進行三次樣條插值，并對其求導，得到T(h)的插值導函數；然后將給定h=[02.34.99.113.718.322.927.2];T=[22.822.822.820.613.911.711.111.1];hh=0:0.1:27.2;pp=spline(h,T);dT=ppd(pp);dTT=ppval(dT,hh);[dTTmax,i]=max(abs(dTT)),hh(i)plot(hh,dTT,'b',hh(i),dTT(i),'r.'),gridon時溫度變化最大，如圖10.2所示(黑點為溫度變化最大的點)?？紤]定積分jbfxdx(10.6)a如果被積函數f(x)是以列表形式給出，則其求解思想同數值微分類似，即用逼近多項式P(x)近似地代替被積函數f(x)，然后計算積分jbP(x)dx，得(10.6)式的近似值；anna函數代替被積函數f(x)，然后積分得(10.6)式的近似值。這兩種類型最終都可歸結為函數f(x)在節(jié)點x上的函數值f(x)的某種線性組合，即下面數值求積公式：kkI=jbf(x)dxnAf(x)或(10.7)akkak=0......aak=0 .8)其中R[f]為截斷誤差。此誤差可用代數精度衡量，代數精度越高，誤差越?。环粗`代數多項式，(10.7)式等號成立；而當f(x)是一個m+1次多項式時，(10.7)式不能精確成立，則稱(10.7)式的代數精度為m。選取不同的近似函數，可產生不同的數值求積公式，常見的有：梯形公式、辛普森公式MATLAB提供了下面幾個函數計算積分，其使用格式分別為： (1)trapz(x)采用梯形公式計算積分(h=1),x為f(k=0,1,,n)k(2)quad('fun',a,b,tol)采用自適應Simpson法計算積分(3)quadl('fun',a,b,tol)采用自適應Gauss-Lobatto法計算積分0解：先對積分作符號運算，然后將其計算結果轉換為數值型，再將其與這三種方法求symsxxz0=simple(int('sqrt(1+xx^2)',0,1))z=double(z0);z=vpa(z,8)x=0:0.01:1;y=sqrt(1+x.^2);z1=trapz(y)*0.01;z1=vpa(z1,8),err1=z-z1;err1=vpa(err1,8)z2=quad('sqrt(1+x.^2)',0,1);z2=vpa(z2,8),err2=z-z2;err2=vpa(err2,8)z3=quadl('sqrt(1+x.^2)',0,1);z3=vpa(z3,8),err3=z-z3;err3=vpa(err3,8)運行得精確值為1(2ln(21))=1.1477936，三種公式計算得數值積分值分別為2軸，則根據所給數據知a=6371+2384=8755,b=6371+439=6810。由對弧長的曲線積分知識知，橢圓的長度為......1L=4j2(a2sin2t+b2cos2t)2dt0上積分稱為橢圓積分，它無法用解析方法計算，可用計算其數值解，編寫M函數文件如下：functiony=y(t)a=8755;b=6810;y=4*sqrt(a^2*sin(t).^2+b^2*cos(t).^2);l=quad('y',0,pi/2)對于用列表形式給出的函數，上述方法不再適用，可利用指令spline構造三次樣條插值函數，再計算積分，具體步驟可參考例2。例3在橋梁的一端每隔一段時間記錄1min有幾輛車過橋，得到數據見表10-3：22025857924:00車輛數/輛9893試估計一天通過橋梁的車流量。解：記記錄時刻為t時，相應的車輛數為x(t)，則所求車流量即為計算積分j24x(t)dt，0x=[0,2,4,5,6,7,8,9,10.3,11.3,12.3,14,16,17,18,19,20,21,22,23,24];y=[2,2,0,2,5,8,25,12,5,10,12,7,9,28,22,10,9,11,8,9,3];pp=spline(x,y);s1=quadl('fun',0,24,[],[],pp)%利用三次樣條插值計算積分s2=trapz(x,y)%利用梯形公式計算積分functionvf=fun(x,pp)vf=ppval(pp,x);儲量計算問題4。試估算出該礦區(qū)(1x4,1y5)煤炭的儲量。123456789x坐標(km)1111122222y坐標(km)123451234526.19x坐標(km)3333344444........y坐標(km)1234512345煤炭厚度h(m)23.2826.4829.1412.0414.5819.9523.7315.3518.0116.29圖10.3，此圖經過插值得到)，而煤炭的儲量即為此立體的體積。要計算此立體的體積，可曲頂柱體，用數值方法計算這些細曲頂柱體的體積，再對其求和即得原曲頂柱體的體積。2．模型的建立及求解zh(10.9)(10.9)D由于函數Q(x,y)只給出了一些離散點上的函數值，無法直接計算上述二重積分，所以下面采用數值積分的方法計算其值。由數值積分的知識知，計算定積分有復合梯形公式為jbf(x)dx~h[f(a)+f(b)+2f(x)]a2k(10.10)kk由(10.9)式得(10.11)(10.11)......其中g(x)=jdQ(x,y)dy，則由(10.10)式可得c(10.12)W=jbg(x)dxhx[g(a)+g(b)+2g(x)](10.12)a2jj=1而c2kk=1c2kk=1kck2kkkkk=1有4kkjjjkj=1k=1(10.13)kkjjjkk=0j=0j=0k=0考慮到給定的數據較少，由此產生的誤差較大，所以利用插值后的數據計算(10.13)x=[11111222223333344444]*1000;y=[12345123451234512345]*1000;z=-[13.7225.808.4725.2722.3215.4721.3314.4924.8326.1923.2826.4829.1412.0414.5819.9523.7315.3518.0116.29];hx=10;hy=10;cx=1000:hx:4000;cy=1000:hy:5000;[X,Y]=meshgrid(cx,cy);n=length(cx);m=l