matlab積分例子

1、matlab積分例子【篇一：matlab積分例子】符號(hào)積分由函數(shù)int來(lái)實(shí)現(xiàn).該函數(shù)的一般調(diào)用格式為：int(s):沒(méi)有指定積分變量和積分階數(shù)時(shí),系統(tǒng)按findsym函數(shù)指示的默認(rèn)變量對(duì)被積函數(shù)或符號(hào)表達(dá)式s求不定積分；int(s,v):以v為自變量,對(duì)被積函數(shù)或符號(hào)表達(dá)式s求不定積分；int(s,v,a,b):求定積分運(yùn)算.a,b分別表示定積分的下限和上限.該函數(shù)求被積函數(shù)在區(qū)間a,b上的定積分.a和b可以是兩個(gè)具體的數(shù),也可以是一個(gè)符號(hào)表達(dá)式,還可以是無(wú)窮(inf)o當(dāng)函數(shù)f關(guān)于變量x在閉區(qū)間a,b上可積時(shí),函數(shù)返回一個(gè)定積分結(jié)果.當(dāng)a,b中有一個(gè)是inf時(shí),函數(shù)返回一個(gè)廣義積分.當(dāng)a,b

2、中有一個(gè)符號(hào)表達(dá)式時(shí),函數(shù)返回一個(gè)符號(hào)函數(shù).例：求函數(shù)x八2+y八2+z八2的三重積分.內(nèi)積分上下限都是函數(shù),對(duì)z積分下限是sqrt(x*y),積分上限是x八2*y;對(duì)y積分下限是sqrt(x),積分上限是x八2;對(duì)x的積分下限1,上限是2,求解如下：symsxyz%定義符號(hào)變量f2=int(int(int(xA2+yA2+zA2,z,sqrt(x*y),xA2*y),y,sqrt(x),xA2),x,1,2)%注意定積分的書(shū)寫(xiě)格式f2=1610027357/6563700-6072064/348075*2八(1/2)+14912/4641*2八(1/4)+64/225*2八(3/4)%給出有

3、理數(shù)解vf2=vpa(f2)%給出默認(rèn)精度的數(shù)值解vf2=224.92153573331143159790710032805二、數(shù)值積分1 .數(shù)值積分根本原理求解定積分的數(shù)值方法多種多樣,如簡(jiǎn)單的梯形法、辛普生(simpson)?法、牛頓柯特斯(newton-cotes)法等都是經(jīng)常采用的方法.它們的根本思想都是將整個(gè)積分區(qū)間a,b分成n個(gè)子區(qū)間xi,xi+1,i=1,2,n,其中x1=a,xn+1=b.這樣求定積分問(wèn)題就分解為求和問(wèn)題.2 .數(shù)值積分的實(shí)現(xiàn)方法基于變步長(zhǎng)辛普生法,matlab給出了quad函數(shù)來(lái)求定積分.該函數(shù)的調(diào)用格式為：i,n=quad(fname,a,b,tol,tra

4、ce)基于變步長(zhǎng)、牛頓柯特斯(newton-cotes)法,matlab給出了quadl函數(shù)來(lái)求定積分.該函數(shù)的調(diào)用格式為：i,n=quadl(fname,a,b,tol,trace)其+fname是被積函數(shù)名.a和b分別是定積分的下限和上限.tol用來(lái)限制積分精度,缺省時(shí)取tol=0.001otrace限制是否展現(xiàn)積分過(guò)程,假設(shè)取非0那么展現(xiàn)積分過(guò)程,取0那么不展現(xiàn),缺省時(shí)取trace=0返回參數(shù)i即定積分值,n為被積函數(shù)的調(diào)用次數(shù).例：求函數(shù)exp(-x*x)的定積分,積分下限為0,積分上限為1.fun=inline(exp(-x.*x),x);%用內(nèi)聯(lián)函數(shù)定義被積函數(shù)fnameisim=

5、quad(fun,0,1)%辛普生法isim=0.746824180726425il=quadl(fun,0,1)%牛頓柯特斯法il=0.746824133988447三、梯形法求向量積分trapz(x,y)梯形法沿列方向求函數(shù)y關(guān)于自變量x的積分(向量形式,數(shù)值方京).d=0.001;x=0:d:1;s=d*trapz(exp(-x.八2)s=0.7468或：formatlonggx=0:0.001:1;%x向量,也可以是不等間距y=exp(-x.八2);%y向量,也可以不是由函數(shù)生成的向量s=trapz(x,y);%求向量積分s=0.746824071499185int的積分可以是定積分,

6、也可以是不定積分(即有沒(méi)有積分上下限都可以積)可以得到解析的解,比方你對(duì)x八2積分,得到的結(jié)果是1/3*xA3,這是通過(guò)解析的方法來(lái)解的.如果int(x八2,x,1,2)得到的結(jié)果是7/3quad是數(shù)值積分,它只能是定積分(就是有積分上下限的積分),它是通過(guò)simpson數(shù)值積分來(lái)求得的(并不是通過(guò)解析的方法得到解析解,再將上下限代入,而是用小梯形的面積求和得到的).如果f=inline(x.A2);quad(f,1,2)得到的結(jié)果是2.333333,這個(gè)數(shù)并不是7/3int是符號(hào)解,無(wú)任何誤差,唯一問(wèn)題是計(jì)算速度；quad是數(shù)值解,有計(jì)算精度限制,優(yōu)點(diǎn)是總是能有一定的速度,即總能在一定時(shí)間內(nèi)

7、給出一個(gè)一定精度的解.from:58.192.116.*對(duì)于y=exp(-(x.A2+x+1)/(1+x),被積函數(shù)之原函數(shù)無(wú)封閉解析表達(dá)式,符號(hào)計(jì)算無(wú)法解題,這是符號(hào)計(jì)算有限性,結(jié)果如下：symsxy=exp(-(x.A2+x+1)/(1+x)s=int(y,x,0,inf)y=exp(-xA2-x-1)/(1+x)warning:at58s=int(exp(-xA2-x-1)/(1+x),x=0.inf)只有通過(guò)數(shù)值計(jì)算解法dx=0.05;%采樣間隔x=0:dx:1000;%數(shù)值計(jì)算適合于有限區(qū)間上,取有限個(gè)

8、采樣點(diǎn),只要終值足夠大,精度不受影響y=exp(-(x.A2+x+1)./(1+x);s=dx*cumtrapz(y);%計(jì)算區(qū)間內(nèi)曲線下列圖形面積,為小矩形面積累加得s(end)ans=0.5641%所求定積分值或進(jìn)行編程,積分上限人工輸入,程序如下：%表達(dá)式保存為函數(shù)文件functiony=fxy(x)y=exp(-(x.A2+x+1)./(1+x);%savefxy.m%main主程序clear,clch=.001;p=0;a=0;r=input(請(qǐng)輸入積分上限,r=)whilearp=p+(fxy(a)+fxy(a+h)*h/2;a=a+h;endp=vpa(p,10)運(yùn)行主程序后得到

9、結(jié)果：請(qǐng)輸入積分上限,r=1000r=1000p=.5641346055其它結(jié)果如下：0-1:int=.30676016860-2:int=.45996331590-5:int=.55830682170-10:int=.56409289750-100:int=.56413460550-1000:int=.5641346055from:211.65.33.*在積分函數(shù)中,sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos(n3*pi*z/9);變量e1,e2,e3,n1,n2,n3通過(guò)函數(shù)參數(shù)輸入,如果直接用inline或字符串的形式,那么表

10、達(dá)式中的未知數(shù)有9個(gè),分別是e1,e2,e3,n1,n2,n3,x,y,z.而用匿名函數(shù)時(shí),變量e1,e2,e3,n1,n2,n3就會(huì)以常數(shù)看待,未知數(shù)就只有x,y,z了,可以求三重積分了.完整函數(shù)程序：functionfn(n1,n2,n3)ifn1=0e1=1;elseifn10e1=2;endendifn2=0e2=1;elseifn20e2=2;endendifn3=0e3=1;elseifn30e3=2;endendf=(x,y,z)sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos(n3*pi*z/9);s=triplequad

11、(f,-6,6,-5.5,5.5,-4.5,4.5)%求三重?cái)?shù)值積分將以上保存為fn.m程序文件,即m文件,然后運(yùn)行：fn(1,1,1)s=866.9655from:211.65.33.*三重積分請(qǐng)用三重積分函數(shù)triplequad,與三個(gè)積分上下限對(duì)應(yīng),即x=triplequad(f,-6,6,-5.5,5.5,-4.5,4.5)其中被積函數(shù)f用匿名函數(shù)來(lái)表達(dá),即f=(x,y,z)sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos(n3*pi*z/9);如果直接用inline或字符串的形式,那么表達(dá)式中的未知數(shù)有9個(gè),分別是e1,e2,e

12、3,n1,n2,n3,x,y,z.而用匿名函數(shù)時(shí),變量e1,e2,e3,n1,n2,n3就會(huì)以常數(shù)看待,未知數(shù)就只有x,y,z了.完整函數(shù)程序：functionfn(n1,n2,n3)ifn1=0e1=1;elseifn10e1=2;endendifn2=0e2=1;elseifn20e2=2;endendifn3=0e3=1;elseifn30e3=2;endendf=(x,y,z)sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos(n3*pi*z/9);x=triplequad(f,-6,6,-5.5,5.5,-4.5,4.5)fn(1

13、,1,1)x=866.9655from:58.192.116.*【篇二：matlab積分例子】前言本帖將通過(guò)一個(gè)典型的數(shù)值積分案例來(lái)介紹如何在matlab中有效地計(jì)算特殊函數(shù)(specialfunctions)的數(shù)值積分.特殊函數(shù)的數(shù)值積分(尤其是無(wú)窮積分)比擬常見(jiàn)的一個(gè)問(wèn)題是：積分在滿足收斂條件下積分結(jié)果為nan.如果我們使用的是integral或quadgk,這一問(wèn)題通常會(huì)以警告的形式給出:“warning:infiniteornot-a-numbervalueencountered,伴隨該警告的惡果是,數(shù)值積分的結(jié)果為nan.為了直觀地說(shuō)明問(wèn)題,本帖將針對(duì)一個(gè)在數(shù)學(xué)、物理學(xué)、以及工程應(yīng)用

14、里較為常見(jiàn)的特殊函數(shù)modifiedbesselfunctionofthefirstkind(數(shù)學(xué)符號(hào)記作：,matlab中對(duì)應(yīng)的函數(shù)是besseli),介紹如何克服積分結(jié)果為nan的問(wèn)題.本帖的方法雖然只是針對(duì)這一種特殊函數(shù),但解決問(wèn)題的思路適用于這一類數(shù)值計(jì)算問(wèn)題,希望能對(duì)有興趣的人有幫助！最后,在分享本案例的過(guò)程中,我還會(huì)順便指出matlabr2021a符號(hào)積分系統(tǒng)的一個(gè)計(jì)算錯(cuò)誤,應(yīng)該算是一個(gè)bug.特殊函數(shù)積分案例先給出我要計(jì)算的積分：首先,我們有必要先對(duì)的函數(shù)特性做一個(gè)根本了解(參考).可以看出：是一個(gè)單調(diào)增函數(shù)從(2)式可以得到,類比,我們可以粗略的認(rèn)為,近似服從指數(shù)漸進(jìn)性(嚴(yán)格的

15、說(shuō),當(dāng),).所以,是一個(gè)近似服從指數(shù)增長(zhǎng)的函數(shù).要保證這樣一個(gè)無(wú)窮積分收斂,被積函數(shù)里必須要有一個(gè)衰減速度極快的因子,該因子衰減速度遠(yuǎn)快于的增長(zhǎng)的速率,使得被積函數(shù)整體是一個(gè)快速衰減的函數(shù),積分才有可能收斂.所以,(1)式的被積函數(shù)表達(dá)式里有一個(gè)快速衰減項(xiàng),該項(xiàng)使得被積函數(shù)整體快速衰減(注：在面前,平方項(xiàng)的遞增影響可以忽略不記).下面,我們來(lái)看看直接用matlabintegral函數(shù)計(jì)算的結(jié)果(matlab版本低于2021a的話可以用quadgk代替).integral(x)x.A2,*exp(-x.A2),*besseli(0,x),0,inf)復(fù)制代碼warning:infiniteorn

16、ot-a-numbervalueencountered.infunfunprivateintegralcalciteratescalarvaluedat349infunfunprivateintegralcalcvadaptat132infunfunprivateintegralcalcat83inintegralat88ans=nan可以看到,直接計(jì)算的結(jié)果為nan.產(chǎn)生nan的原因產(chǎn)生nan的原因其實(shí)很簡(jiǎn)單,由于按matlab的運(yùn)算規(guī)那么,一旦計(jì)算過(guò)程中出現(xiàn)了nan,最終結(jié)果就為nan.所以,這里一定是由于計(jì)算被積函數(shù)的函數(shù)值時(shí)得到了nan.簡(jiǎn)單的說(shuō),當(dāng)x很大時(shí)(如1000),快速遞增單獨(dú)

17、計(jì)算的結(jié)果為inf,而快速衰減項(xiàng)的結(jié)果為0,二者相乘,matlab里會(huì)得到inf*0=nan.為了詳細(xì)說(shuō)明問(wèn)題,我們先來(lái)計(jì)算一下被積函數(shù)當(dāng)x=1000的取值.直接用matlab數(shù)值計(jì)算f=(x)xA2.*exp(-xA2).*besseli(0,x)f(1000)復(fù)制代碼ans=nan數(shù)值計(jì)算得到結(jié)果為nan.這正是由于單獨(dú)的兩項(xiàng)besseli(0,x)和exp(-xA2)分別計(jì)算時(shí)得到了inf和0,如下所示：exp(-1000A2)復(fù)制代碼ans=0besseli(0,1000)復(fù)制代碼ans=inf二者再相乘得到nan.既然在計(jì)算函數(shù)值的過(guò)程中就出現(xiàn)了nan,難怪積分名果也是nan.我們?cè)?/p>

18、看看符號(hào)計(jì)算的結(jié)果：vpa(sym(1000A2*exp(-1000A2)*besseli(0,1000),4)復(fù)制代碼ans=8.195e-433857根據(jù)雙精度數(shù)值計(jì)算的精度,這個(gè)結(jié)果在matlab數(shù)值計(jì)算里可以認(rèn)為是0.所以,數(shù)值計(jì)算中如果能將這個(gè)結(jié)果作為0處理,也就能避免積分結(jié)果為nan了.解決nan的方法這里介紹的方法適用于一類具備共同特征的特殊函數(shù)積分：被積函數(shù)里有一個(gè)無(wú)上界的快速遞增項(xiàng)和一個(gè)快速衰減項(xiàng),二者共同作用的結(jié)果使得被積函數(shù)在積分區(qū)間上或者積分變量x取值較大時(shí)快速衰減,從而積分收斂.最簡(jiǎn)單的解決方法：既然問(wèn)題出現(xiàn)在大x下,可以適當(dāng)減小積分上限,將上限inf換成一個(gè)有限的數(shù),該數(shù)必須滿足兩個(gè)條件：1)必須足夠大使得被積函數(shù)整體的取值很小(近似為0),2)又不能太大而導(dǎo)致數(shù)值計(jì)算里出現(xiàn)inf*0=nan.我們可以嘗試將積分上限設(shè)置為100:integral(x)x,A2.*exp(-x.A2).*besseli(0,x),0,100)復(fù)制代碼ans