數(shù)值積分與求解

數(shù)值積分與求解第1頁，課件共32頁，創(chuàng)作于2023年2月

MATLAB求解連續(xù)函數(shù)積分第2頁，課件共32頁，創(chuàng)作于2023年2月引言我們知道,若函數(shù)f(x)在區(qū)間[a,b]上連續(xù)且其原函數(shù)為F(x),則可用Newton-Leibnitz公式求定積分的值,Newton-Leibnitz公式無論在理論上還是在解決實(shí)際問題上都起了很大作用，但它并不能完全解決定積分的計(jì)算問題，因?yàn)榉e分學(xué)涉及的實(shí)際問題極為廣泛，而且極其復(fù)雜，在實(shí)際計(jì)算中經(jīng)常遇到以下三種情況：1Matlab求解連續(xù)函數(shù)積分第3頁，課件共32頁，創(chuàng)作于2023年2月

(1)被積函數(shù)f(x)并不一定能夠找到用初等函數(shù)的有限形式表示的原函數(shù)F(x)，例如：Newton-Leibnitz公式就無能為力了(2)還有被積函數(shù)f(x)的原函數(shù)能用初等函數(shù)表示，但表達(dá)式太復(fù)雜，例如函數(shù)并不復(fù)雜，但積分后其表達(dá)式卻很復(fù)雜，積分后其原函數(shù)F(x)為：第4頁，課件共32頁，創(chuàng)作于2023年2月(3)被積函數(shù)f(x)沒有具體的解析表達(dá)式,其函數(shù)關(guān)系由表格或圖形表示。對(duì)于這些情況,要計(jì)算積分的準(zhǔn)確值都是十分困難的。由此可見,通過原函數(shù)來計(jì)算積分有它的局限性,因而研究一種新的積分方法來解決Newton-Leibniz公式所不能或很難解決的積分問題,這時(shí)需要用數(shù)值解法來建立積分的近似計(jì)算方法。

將積分區(qū)間細(xì)分,在每一個(gè)小區(qū)間內(nèi)用簡(jiǎn)單函數(shù)代替復(fù)雜函數(shù)進(jìn)行積分，這就是數(shù)值積分的思想，用代數(shù)插值多項(xiàng)式去代替被積函數(shù)發(fā)f(x)進(jìn)行積分是本章討論數(shù)值積分的主要內(nèi)容。

第5頁，課件共32頁，創(chuàng)作于2023年2月1.1定積分的Matlab符號(hào)計(jì)算例1由y=sinx，y=cosx,x=-1/2,x=3/2所圍成的平面區(qū)域D.求平面區(qū)域D的面積S.解輸入作函數(shù)圖形的程序>>x=-1:0.001:2;F1=sin(x);F2=cos(x);plot(x,F1,'b-',x,F2,'g-'),axis([-1,pi/4+1,-1.3,1.3]),xlabel('x'),ylabel('y'),title('y=sinx,y=cosx和x=-0.5及x=1.5所圍成的平面區(qū)域的圖形')運(yùn)行后屏幕顯示圖形.求平面區(qū)域D的面積S.輸入計(jì)算面積S的程序>>symsxf1=cos(x)-sin(x);f2=-f1;S1=int(f1,x,-0.5,pi/4);S2=int(f2,x,pi/4,1.5);S=S1+S2,Sj=double(S)運(yùn)行后屏幕顯示計(jì)算面積的值S及其近似值Sj如下S=2*2^(1/2)+sin(1/2)-cos(1/2)-sin(3/2)-cos(3/2)Sj=1.36203791318826坐標(biāo)調(diào)整第6頁，課件共32頁，創(chuàng)作于2023年2月1.2變限積分的Matlab符號(hào)計(jì)算例2已知,求F′(x)解：輸入程序：>>symsxtF1=int(exp(t)*sin(2+sqrt(t^3)),x,0);F2=int(exp(t)*sin(2+sqrt(t^3)),0,x^2);Fi=F1+F2;dF=diff(Fi)運(yùn)行后屏幕顯示計(jì)算變限積分F(x)的導(dǎo)數(shù).第7頁，課件共32頁，創(chuàng)作于2023年2月

建立數(shù)值積分公式的途徑比較多,其中最常用的有兩種：（1）由積分中值定理可知，對(duì)于連續(xù)函數(shù)f(x)，在積分區(qū)間[a,b]內(nèi)存在一點(diǎn)ξ，使得即所求的曲邊梯形的面積恰好等于底為(b-a)，高為的矩形面積。但是點(diǎn)ξ的具體位置一般是未知的,因而的值也是未知的,稱為f(x)在區(qū)間[a,b]上的平均高度。那么只要對(duì)平均高度提供一種算法，相應(yīng)地就獲得一種數(shù)值求積方法1.3數(shù)值求積方法第8頁，課件共32頁，創(chuàng)作于2023年2月①梯形公式②矩形公式按照這種思想，可構(gòu)造出一些求積分值的近似公式。例如分別取和則分別得到中矩形公式和梯形公式。第9頁，課件共32頁，創(chuàng)作于2023年2月③

Simpson公式

矩形公式把[a,b]的中點(diǎn)處函數(shù)值作為平均高度f()的近似值而獲得的一種數(shù)值積分方法。

在這三個(gè)公式中,梯形公式把f(a),f(b)的加權(quán)平均值

作為平均高度f()的近似值而獲得的一種數(shù)值積分方法。

第10頁，課件共32頁，創(chuàng)作于2023年2月Simpson公式是以函數(shù)f(x)在a,b,(a+b)/2這三點(diǎn)的函數(shù)值f(a),f(b),的加權(quán)平均值似值而獲得的一種數(shù)值積分方法。作為平均高度f()的近（2）先用某個(gè)簡(jiǎn)單函數(shù)近似逼近f(x),用代替原被積函數(shù)f(x)，即以此構(gòu)造數(shù)值算法。從數(shù)值計(jì)算的角度考慮,函數(shù)應(yīng)對(duì)f(x)有充分的逼近程度,并且容易計(jì)算其積分。由于多項(xiàng)式能很好地逼近連續(xù)函數(shù),且又容易計(jì)算積分,因此將選取為插值多項(xiàng)式,這樣f(x)的積分就可以用其插值多項(xiàng)式的積分來近似代替第11頁，課件共32頁，創(chuàng)作于2023年2月(一)用函數(shù)trapz計(jì)算定積分調(diào)用格式一：Z=trapz(Y)調(diào)用格式二：Z=trapz(X,Y)調(diào)用格式三：Z=trapz(X,Y,DIM)或trapz(Y,DIM)(二)用函數(shù)cumtrapz計(jì)算定積分調(diào)用格式一：Z=cumtrapz(Y)調(diào)用格式二：Z=cumtrapz(X,Y)調(diào)用格式三：Z=cumtrapz(X,Y,DIM)或cumtrapz(Y,DIM)1.3.1梯形公式的Matlab程序第12頁，課件共32頁，創(chuàng)作于2023年2月梯形數(shù)值積分的MATLAB主程序functionT=rctrap(fun,a,b,m)n=1;h=b-a;T=zeros(1,m+1);x=a;T(1)=h*(feval(fun,a)+feval(fun,b))/2;fori=1:mh=h/2;n=2*n;s=0;fork=1:n/2x=a+h*(2*k-1);s=s+feval(fun,x);endT(i+1)=T(i)/2+h*s;endT=T(1:m);第13頁，課件共32頁，創(chuàng)作于2023年2月例3用MATLAB的函數(shù)trapz和cumtrapz分別計(jì)算

精確到10-4

，并與矩形公式比較.解：將[0,

π/2]分成20等份，步長(zhǎng)為π/40，輸入程序如下（注意trapz(y)是單位步長(zhǎng),trapz(y)*h=trapz(x,y)）：>>h=pi/40;x=0:h:pi/2;y=exp(-x).*sin(x);z1=sum(y(1:20))*h,z2=sum(y(2:21))*h,z=(z1+z2)/2z3=trapz(y)*h,z3h=trapz(x,y),z3c=cumtrapz(y)*h,運(yùn)行后屏幕顯示用矩形公式（9.3），（9.4）計(jì)算結(jié)果z1、z2和二者的平均數(shù)z、函數(shù)trapz和cumtrapz分別計(jì)算結(jié)果z3、z3c.第14頁，課件共32頁，創(chuàng)作于2023年2月（一）函數(shù)sum的調(diào)用格式調(diào)用格式一：sum(X)調(diào)用格式二：sum(X,DIM)（二）函數(shù)cumsum的調(diào)用格式調(diào)用格式一：cumsum(X)調(diào)用格式二：cumsum(X,DIM)1.3.2矩形公式的Matlab程序第15頁，課件共32頁，創(chuàng)作于2023年2月例4用MATLAB的函數(shù)sum和cumsum及矩形公式計(jì)算

，并與精確值比較.解：將[0,/2]分成20等份，步長(zhǎng)為/40，輸入程序如下（注意sum和cumsum的用法）>>h=pi/40;x=0:h:pi/2;y=exp(-x).*sin(x);z1=sum(y(1:20))*h,z2=sum(y(2:21))*h,z=cumsum(y);z11=z(20)*h,z12=(z(21)-z(1))*h,運(yùn)行后屏幕顯示計(jì)算結(jié)果分別如下z1=z2=z11=z12=0.38730.40360.38730.4036求定積分的精確值，輸入程序>>symsxF=int(exp(-x)*sin(x),x,0,pi/2)Fs=double(F),wz1=abs(Fs-z1),wz2=abs(Fs-z2)運(yùn)行后屏幕顯示定積分的精確值Fs和用矩形公式計(jì)算結(jié)果的絕對(duì)誤差wz1、wz2分別如下F=Fs=1/2*(-1+exp(pi)^(1/2))/exp(pi)^(1/2)0.3961wz1=wz2=0.00880.0075第16頁，課件共32頁，創(chuàng)作于2023年2月調(diào)用格式一：quad(‘fun’,a,b)調(diào)用格式二：quad(‘fun’,a,b,tol)調(diào)用格式三：[Q,FCNT]=quad(...)調(diào)用格式四：quad(‘fun’,a,b,tol,TRACE)調(diào)用格式五：quad(‘fun’,a,b,tol,TRACE,P1,P2,…)復(fù)合辛普森（Simpson）數(shù)值積分的MATLAB主程序functiony=comsimpson(fun,a,b,n)z1=feval(fun,a)+feval(fun,b);m=n/2;h=(b-a)/(2*m);x=a;z2=0;z3=0;x2=0;x3=0;fork=2:2:2*mx2=x+k*h;z2=z2+2*feval(fun,x2);endfork=3:2:2*mx3=x+k*h;z3=z3+4*feval(fun,x3);endy=(z1+z2+z3)*h/3;1.3.3辛普森公式的Matlab程序第17頁，課件共32頁，創(chuàng)作于2023年2月例5用comsimpson.m和quad.m分別計(jì)算定積分，取精度為10-4

，并與精確值比較.解：輸入程序>>[Q1,FCNT14]=quad(@fun,0,1,1.e-4,3),Q2=comsimpson(@fun,0,1,10000)symsxfi=int(exp((-x.^2)./2)./(sqrt(2*pi)),x,0,1);Fs=double(fi)wQ1=double(abs(fi-Q1)),wQ2=double(abs(fi-Q2))運(yùn)行后屏幕顯示I的精確值Fs，用comsimpson.m和quad.m分別計(jì)算I的近似值Q2、Q1和迭代次數(shù)FCNT14，取精度分別為104，Q2、Q1分別與精確值Fs的絕對(duì)誤差wQ2,wQ1第18頁，課件共32頁，創(chuàng)作于2023年2月

MATLAB求解離散函數(shù)積分第19頁，課件共32頁，創(chuàng)作于2023年2月

1、S=trapz(X,Y)等距梯形法求數(shù)值積分調(diào)用格式一：Z=trapz(Y)調(diào)用格式二：Z=trapz(X,Y)調(diào)用格式三：Z=trapz(X,Y,DIM)或trapz(Y,DIM)

2、S=cumsum(Y)歐拉法求數(shù)值積分調(diào)用格式一：cumsum(X)調(diào)用格式二：cumsum(X,DIM)

3、對(duì)于離散點(diǎn)的積分，先對(duì)其擬合獲得函數(shù)表達(dá)式，再作為連續(xù)被積函數(shù)求積分。2Matlab求解離散函數(shù)積分第20頁，課件共32頁，創(chuàng)作于2023年2月例6求解：（1）輸入程序>>symsxF=limit(1/(sqrt(1-x^2)),x,1,'left')運(yùn)行后屏幕顯示結(jié)果為F=inf即當(dāng)x→1-時(shí)，被積函數(shù)→∞（2）輸入程序>>symsxF1=int(1/(sqrt(1-x^2)),x,0,1),LimF1=double(F1)運(yùn)行后屏幕顯示結(jié)果為F1=LimF1=1/2*pi1.5708第21頁，課件共32頁，創(chuàng)作于2023年2月例7求解：（1）因?yàn)楸环e函數(shù)在[-1,1]上除x=0外連續(xù)，輸入程序>>symsxF=limit((3*x^5-8)/(x^2),x,0)運(yùn)行后屏幕顯示結(jié)果為F=-inf（2）輸入程序>>symsxF1=int((3*x^5-8)/(x^2),x,-1,0),F2=int((3*x^5-8)/(x^2),x,0,1),F=F1+F2運(yùn)行后屏幕顯示結(jié)果為F1=-infF2=35/4F=-inf第22頁，課件共32頁，創(chuàng)作于2023年2月例8討論反常積分

的斂散性.解：（1）因?yàn)楸环e函數(shù)1/xq在(0,1]上連續(xù)，在x=0無定義，輸入>>symsxLF05=limit(1/(x^0.5),x,0,'right'),LF1=limit(1/x,x,0,'right')LF2=limit(1/(x^2),x,0,'right')運(yùn)行后屏幕顯示結(jié)果為L(zhǎng)F05=InfLF1=InfLF2=Inf（2）輸入程序>>symsxF05=int(1/(x^0.5),x,0,1),F1=int(1/x,x,0,1),F2=int(1/(x^2),x,0,1)運(yùn)行后屏幕顯示結(jié)果為F05=2F1=InfF2=Inf第23頁，課件共32頁，創(chuàng)作于2023年2月MATLAB求解多重積分第24頁，課件共32頁，創(chuàng)作于2023年2月3.1二重積分的數(shù)值求解

使用MATLAB提供的dblquad函數(shù)就可以直接求出一般域上二重定積分的數(shù)值解。該函數(shù)的調(diào)用格式為：

I=dblquad(f,a,b,c,d,tol,trace)該函數(shù)求f(x,y)在[a,b]×[c,d]區(qū)域上的二重定積分。參數(shù)tol，trace的用法與函數(shù)quad完全相同。3Matlab求解多重積分第25頁，課件共32頁，創(chuàng)作于2023年2月例9用MATLAB函數(shù)dblquad求直徑為8的半球的體積V球，誤差為10-4.解:在MATLAB工作窗口輸入下列MATLAB程序>>a=-4;b=4;c=-4;d=4;V1=dblquad(inline('sqrt(max(4^2-(x.^2+y.^2),0))'),a,b,c,d,1.e-4,@quadl)V=dblquad(inline('sqrt(4^2-(x.^2+y.^2)).*(x.^2+y.^2<=4^2)'),a,b,c,d,1.e-4)symstwbjh=sqrt(4^2-(w.^2+t.^2));y1=-sqrt(4^2-t.^2);y2=sqrt(4^2-t.^2);jfx=int(bjh,w,y1,y2);jfy=int(jfx,t,a,b);I2=double(jfy),JuewuL1=abs(I2-V1)Juewu1=abs(I2-V),ezplot('x^2+y^2-16',[-5,5]);axisequal第26頁，課件共32頁，創(chuàng)作于2023年2月運(yùn)行后屏幕顯示如下V1=V=1.340434607608882e+0021.340477821376317e+002Warning:Explicitintegralcouldnotbefound.I2=JuewuL1=1.340412865531645e+0020.00217420772367Juewu1=0.00649558446713第27頁，課件共32頁，創(chuàng)作于2023年2月3.2三重積分的數(shù)值求解(一)用MATLAB符號(hào)計(jì)算三重積分（1）畫出積分區(qū)域V的草圖.輸入程序（2）確定積分限.輸入程序（3）輸入計(jì)算程序例10計(jì)算，其中積分區(qū)域V是由旋轉(zhuǎn)拋物面Z=8-x2-y2，圓柱x2+y2=4和z=0所圍成的空間閉區(qū)域.第28頁，課件共32頁，創(chuàng)作于2023年2月解：（1）畫出積分區(qū)域V的草圖.輸入程序>>[x,y]=meshgrid(-2:0.01:2);z1=8-(x.^2+y.^2);figure(1)meshc(x,y,z1)holdonx=-2:0.01:2;r=2;[x,y,z]=cylinder(r,30)mesh(x,y,z)holdofftitle('由旋轉(zhuǎn)拋物面z=8-(x^2+y^2)，圓柱面x^2+y^2=4和z=0所圍成的積分區(qū)域V')figure(2)contour(x,y,z,10)title('由z=8-(x^2+y^2)，圓柱面x^2+y^2=4和z=0所圍成區(qū)域V在x0y面上的投影區(qū)域Dxy')運(yùn)行后屏幕顯示圖形.第29頁，課件共32頁，創(chuàng)作于2023年2月（2）確定積分限.輸入程序>>symsxyzf1=('z=8-(x^2+y^2)');f2=('x^2+y^2=4');[x,y,z]=solve(f1,f2,x,y,z)運(yùn)行后屏幕顯示旋轉(zhuǎn)拋物面和圓柱面的交線如下x=y=z=[(4-y^2)^(1/2)][y][4][-(4-y^2)^