哈工大-數(shù)值分析上機實驗報告

..實驗報告一題目：非線性方程求解摘要：非線性方程的解析解通常很難給出,因此線性方程的數(shù)值解法就尤為重要。本實驗采用兩種常見的求解方法二分法和Newton法及改進的Newton法。前言：〔目的和意義掌握二分法與Newton法的基本原理和應(yīng)用。數(shù)學(xué)原理：對于一個非線性方程的數(shù)值解法很多。在此介紹兩種最常見的方法：二分法和Newton法。對于二分法,其數(shù)學(xué)實質(zhì)就是說對于給定的待求解的方程f<x>,其在[a,b]上連續(xù),f<a>f<b><0,且f<x>在[a,b]內(nèi)僅有一個實根x*,取區(qū)間中點c,若,則c恰為其根,否則根據(jù)f<a>f<c><0是否成立判斷根在區(qū)間[a,c]和[c,b]中的哪一個,從而得出新區(qū)間,仍稱為[a,b]。重復(fù)運行計算,直至滿足精度為止。這就是二分法的計算思想。Newton法通常預(yù)先要給出一個猜測初值x0,然后根據(jù)其迭代公式產(chǎn)生逼近解x*的迭代數(shù)列{xk},這就是Newton法的思想。當(dāng)x0接近x*時收斂很快,但是當(dāng)x0選擇不好時,可能會發(fā)散,因此初值的選取很重要。另外,若將該迭代公式改進為其中r為要求的方程的根的重數(shù),這就是改進的Newton法,當(dāng)求解已知重數(shù)的方程的根時,在同種條件下其收斂速度要比Newton法快的多。程序設(shè)計：本實驗采用Matlab的M文件編寫。其中待求解的方程寫成function的方式,如下functiony=f<x>;y=-x*x-sin<x>;寫成如上形式即可,下面給出主程序。二分法源程序：clear%%%給定求解區(qū)間b=1.5;a=0;%%%誤差R=1;k=0;%迭代次數(shù)初值while<R>5e-6>;c=<a+b>/2;iff12<a>*f12<c>>0;a=c;elseb=c;endR=b-a;%求出誤差k=k+1;endx=c%給出解Newton法及改進的Newton法源程序：clear%%%%輸入函數(shù)f=input<'請輸入需要求解函數(shù)>>','s'>%%%求解f<x>的導(dǎo)數(shù)df=diff<f>;%%%改進常數(shù)或重根數(shù)miu=2;%%%初始值x0x0=input<'inputinitialvaluex0>>'>;k=0;%迭代次數(shù)max=100;%最大迭代次數(shù)R=eval<subs<f,'x0','x'>>;%求解f<x0>,以確定初值x0時否就是解while<abs<R>>1e-8>x1=x0-miu*eval<subs<f,'x0','x'>>/eval<subs<df,'x0','x'>>;R=x1-x0;x0=x1;k=k+1;if<eval<subs<f,'x0','x'>><1e-10>;breakendifk>max;%如果迭代次數(shù)大于給定值,認為迭代不收斂,重新輸入初值ss=input<'mayberesultiserror,chooseanewx0,y/n?>>','s'>;ifstrcmp<ss,'y'>x0=input<'inputinitialvaluex0>>'>;k=0;elsebreakendendendk;%給出迭代次數(shù)x=x0；%給出解結(jié)果分析和討論：用二分法計算方程在[1,2]內(nèi)的根。<,下同>計算結(jié)果為x=1.40441513061523；f<x>=-3.797205105904311e-007；k=18；由f<x>知結(jié)果滿足要求,但迭代次數(shù)比較多,方法收斂速度比較慢。用二分法計算方程在[1,1.5]內(nèi)的根。計算結(jié)果為x=1.32471847534180；f<x>=2.209494846194815e-006；k=17；由f<x>知結(jié)果滿足要求,但迭代次數(shù)還是比較多。用Newton法求解下列方程x0=0.5；計算結(jié)果為x=；f<x>=2.220446049250313e-016；k=4；由f<x>知結(jié)果滿足要求,而且又迭代次數(shù)只有4次看出收斂速度很快。x0=1；x0=0.45,x0=0.65；當(dāng)x0=0.45時,計算結(jié)果為x=0.49999999999983；f<x>=-8.362754932994584e-014；k=4；由f<x>知結(jié)果滿足要求,而且又迭代次數(shù)只有4次看出收斂速度很快,實際上該方程確實有真解x=0.5。當(dāng)x0=0.65時,計算結(jié)果為x=0.50000000000000；f<x>=0；k=9；由f<x>知結(jié)果滿足要求,實際上該方程確實有真解x=0.5,但迭代次數(shù)增多,實際上當(dāng)取x0〉0.68時,x≈1,就變成了方程的另一個解,這說明Newton法收斂與初值很有關(guān)系,有的時候甚至可能不收斂。用改進的Newton法求解,有2重根,取x0=0.55；并與3.中的c>比較結(jié)果。當(dāng)x0=0.55時,程序死循環(huán),無法計算,也就是說不收斂。改時,結(jié)果收斂為x=0.50000087704286；k=16；顯然這個結(jié)果不是很好,而且也不是收斂至方程的2重根上。當(dāng)x0=0.85時,結(jié)果收斂為x=1.00000000000489；f<x>=2.394337647718737e-023；k=4；這次達到了預(yù)期的結(jié)果,這說明初值的選取很重要,直接關(guān)系到方法的收斂性,實際上直接用Newton法,在給定同樣的條件和精度要求下,可得其迭代次數(shù)k=15,這說明改進后的Newton法法速度確實比較快。結(jié)論：對于二分法,只要能夠保證在給定的區(qū)間內(nèi)有根,使能夠收斂的,當(dāng)時收斂的速度和給定的區(qū)間有關(guān),二且總體上來說速度比較慢。Newton法,收斂速度要比二分法快,但是最終其收斂的結(jié)果與初值的選取有關(guān),初值不同,收斂的結(jié)果也可能不一樣,也就是結(jié)果可能不時預(yù)期需要得結(jié)果。改進的Newton法求解重根問題時,如果初值不當(dāng),可能會不收斂,這一點非常重要,當(dāng)然初值合適,相同情況下其速度要比Newton法快得多。實驗報告二題目：Gauss列主元消去法摘要：求解線性方程組的方法很多,主要分為直接法和間接法。本實驗運用直接法的Guass消去法,并采用選主元的方法對方程組進行求解。前言：〔目的和意義學(xué)習(xí)Gauss消去法的原理。了解列主元的意義。確定什么時候系數(shù)陣要選主元數(shù)學(xué)原理：由于一般線性方程在使用Gauss消去法求解時,從求解的過程中可以看到,若=0,則必須進行行交換,才能使消去過程進行下去。有的時候即使0,但是其絕對值非常小,由于機器舍入誤差的影響,消去過程也會出現(xiàn)不穩(wěn)定得現(xiàn)象,導(dǎo)致結(jié)果不正確。因此有必要進行列主元技術(shù),以最大可能的消除這種現(xiàn)象。這一技術(shù)要尋找行r,使得并將第r行和第k行的元素進行交換,以使得當(dāng)前的的數(shù)值比0要大的多。這種列主元的消去法的主要步驟如下：消元過程對k=1,2,…,n-1,進行如下步驟。選主元,記若很小,這說明方程的系數(shù)矩陣嚴重病態(tài),給出警告,提示結(jié)果可能不對。交換增廣陣A的r,k兩行的元素。<j=k,…,n+1>計算消元<i=k+1,…,n;j=k+1,……,n+1>回代過程對k=n,n-1,…,1,進行如下計算至此,完成了整個方程組的求解。程序設(shè)計：本實驗采用Matlab的M文件編寫。Gauss消去法源程序：cleara=input<'輸入系數(shù)陣：>>\n'>b=input<'輸入列陣b：>>\n'>n=length<b>;A=[ab]x=zeros<n,1>;%%%函數(shù)主體fork=1:n-1;%%%是否進行主元選取ifabs<A<k,k>><yipusilong;%事先給定的認為有必要選主元的小數(shù)yzhuyuan=1;elseyzhuyuan=0;endifyzhuyuan;%%%%選主元t=A<k,k>;forr=k+1:n;ifabs<A<r,k>>>abs<t>p=r;elsep=k;endend%%%交換元素ifp~=k;forq=k:n+1;s=A<k,q>;A<k,q>=A<p,q>;A<p,q>=s;endendend%%%判斷系數(shù)矩陣是否奇異或病態(tài)非常嚴重ifabs<A<k,k>><yipusilongdisp<‘矩陣奇異,解可能不正確’>end%%%%計算消元,得三角陣forr=k+1:n;m=A<r,k>/A<k,k>;forq=k:n+1;A<r,q>=A<r,q>-A<k,q>*m;endendend%%%%求解xx<n>=A<n,n+1>/A<n,n>;fork=n-1:-1:1;s=0;forr=k+1:n;s=s+A<k,r>*x<r>;endt=<A<k,n+1>-s>x<k>=<A<k,n+1>-s>/A<k,k>end結(jié)果分析和討論：例：求解方程。其中為一小數(shù),當(dāng)時,分別采用列主元和不列主元的Gauss消去法求解,并比較結(jié)果。記Emax為求出的解代入方程后的最大誤差,按要求,計算結(jié)果如下：當(dāng)時,不選主元和選主元的計算結(jié)果如下,其中前一列為不選主元結(jié)果,后一列為選主元結(jié)果,下同。0.999999347683910.999999347826512.000002174219722.000002173911632.999997608594512.99999760869721Emax=,0此時,由于不是很小,機器誤差就不是很大,由Emax可以看出不選主元的計算結(jié)果精度還可以,因此此時可以考慮不選主元以減少計算量。當(dāng)時,不選主元和選主元的計算結(jié)果如下1.000017846308770.999999999993481.999980097208072.000000000021743.000006634247312.99999999997609Emax=2.036758973744668e-005,0此時由Emax可以看出不選主元的計算精度就不好了,誤差開始增大。當(dāng)時,不選主元和選主元的計算結(jié)果如下1.666666666666662.00000000000000300000000000000Emax=0.70770085900503,0此時由Emax可以看出,不選主元的結(jié)果應(yīng)該可以說是不正確了,這是由機器誤差引起的。當(dāng)時,不選主元和選主元的計算結(jié)果如下NaN1NaN2NaN3Emax=NaN,0不選主元時,程序報錯：Warning:Dividebyzero.。這是因為機器計算的最小精度為10-15,所以此時的就認為是0,故出現(xiàn)了錯誤現(xiàn)象。而選主元時則沒有這種現(xiàn)象,而且由Emax可以看出選主元時的結(jié)果應(yīng)該是精確解。結(jié)論：采用Gauss消去法時,如果在消元時對角線上的元素始終較大〔假如大于10-5,那么本方法不需要進行列主元計算,計算結(jié)果一般就可以達到要求,否則必須進行列主元這一步,以減少機器誤差帶來的影響,使方法得出的結(jié)果正確。實驗報告三題目：Rung現(xiàn)象產(chǎn)生和克服摘要：由于高次多項式插值不收斂,會產(chǎn)生Runge現(xiàn)象,本實驗在給出具體的實例后,采用分段線性插值和三次樣條插值的方法有效的克服了這一現(xiàn)象,而且還取的很好的插值效果。前言：〔目的和意義深刻認識多項式插值的缺點。明確插值的不收斂性怎樣克服。明確精度與節(jié)點和插值方法的關(guān)系。數(shù)學(xué)原理：在給定n+1個節(jié)點和相應(yīng)的函數(shù)值以后構(gòu)造n次的Lagrange插值多項式,實驗結(jié)果表明〔見后面的圖這種多項式并不是隨著次數(shù)的升高對函數(shù)的逼近越來越好,這種現(xiàn)象就是Rung現(xiàn)象。解決Rung現(xiàn)象的方法通常有分段線性插值、三次樣條插值等方法。分段線性插值：設(shè)在區(qū)間[a,b]上,給定n+1個插值節(jié)點a=x0<x1<…<xn=b和相應(yīng)的函數(shù)值y0,y1,…,yn,,求作一個插值函數(shù),具有如下性質(zhì)：,j=0,1,…,n。在每個區(qū)間[xi,xj]上是線性連續(xù)函數(shù)。則插值函數(shù)稱為區(qū)間[a,b]上對應(yīng)n個數(shù)據(jù)點的分段線性插值函數(shù)。三次樣條插值：給定區(qū)間[a,b]一個分劃⊿：a=x0<x1<…<xN=b若函數(shù)S<x>滿足下列條件：S<x>在每個區(qū)間[xi,xj]上是不高于3次的多項式。S<x>及其2階導(dǎo)數(shù)在[a,b]上連續(xù)。則稱S<x>使關(guān)于分劃⊿的三次樣條函數(shù)。程序設(shè)計：本實驗采用Matlab的M文件編寫。其中待插值的方程寫成function的方式,如下functiony=f<x>;y=1/〔1+25*x*x;寫成如上形式即可,下面給出主程序Lagrange插值源程序：n=input<'將區(qū)間分為的等份數(shù)輸入：\n'>;s=[-1+2/n*[0:n]];%%%給定的定點,Rf為給定的函數(shù)x=-1:0.01:1;f=0;forq=1:n+1;l=1;%求插值基函數(shù)fork=1:n+1;ifk~=q;l=l.*<x-s<k>>./<s<q>-s<k>>;elsel=l;endendf=f+Rf<s<q>>*l;%求插值函數(shù)endplot<x,f,'r'>%作出插值函數(shù)曲線gridonholdon分段線性插值源程序clearn=input<'將區(qū)間分為的等份數(shù)輸入：\n'>;s=[-1+2/n*[0:n]];%%%給定的定點,Rf為給定的函數(shù)m=0;hh=0.001;forx=-1:hh:1;ff=0;fork=1:n+1;%%%求插值基函數(shù)switchkcase1ifx<=s<2>;l=<x-s<2>>./<s<1>-s<2>>;elsel=0;endcasen+1ifx>s<n>;l=<x-s<n>>./<s<n+1>-s<n>>;elsel=0;endotherwiseifx>=s<k-1>&x<=s<k>;l=<x-s<k-1>>./<s<k>-s<k-1>>;elseifx>=s<k>&x<=s<k+1>;l=<x-s<k+1>>./<s<k>-s<k+1>>;elsel=0;endendendff=ff+Rf<s<k>>*l;%%求插值函數(shù)值endm=m+1;f<m>=ff;end%%%作出曲線x=-1:hh:1;plot<x,f,'r'>;gridonholdon三次樣條插值源程序：〔采用第一邊界條件clearn=input<'將區(qū)間分為的等份數(shù)輸入：\n'>;%%%插值區(qū)間a=-1;b=1;hh=0.001;%畫圖的步長s=[a+<b-a>/n*[0:n]];%%%給定的定點,Rf為給定的函數(shù)%%%%第一邊界條件Rf"<-1>,Rf"<1>v=5000*1/<1+25*a*a>^3-50/<1+25*a*a>^4;fork=1:n;%取出節(jié)點間距h<k>=s<k+1>-s<k>;endfork=1:n-1;%求出系數(shù)向量lamuda,miula<k>=h<k+1>/<h<k+1>+h<k>>;miu<k>=1-la<k>;end%%%%賦值系數(shù)矩陣Afork=1:n-1;forp=1:n-1;switchpcasekA<k,p>=2;casek-1A<k,p>=miu<p+1>;casek+1A<k,p>=la<p-1>;otherwiseA<k,p>=0;endendend%%%%求出d陣fork=1:n-1;switchkcase1d<k>=6*f2c<[s<k>s<k+1>s<k+2>]>-miu<k>*v;casen-1d<k>=6*f2c<[s<k>s<k+1>s<k+2>]>-la<k>*v;otherwised<k>=6*f2c<[s<k>s<k+1>s<k+2>]>;endend%%%%求解M陣M=A\d';M=[v;M;v];%%%%m=0;f=0;forx=a:hh:b;ifx==a;p=1;elsep=ceil<<x-s<1>>/<<b-a>/n>>;endff1=0;ff2=0;ff3=0;ff4=0;m=m+1;ff1=1/h<p>*<s<p+1>-x>^3*M<p>/6;ff2=1/h<p>*<x-s<p>>^3*M<p+1>/6;ff3=<<Rf<s<p+1>>-Rf<s<p>>>/h<p>-h<p>*<M<p+1>-M<p>>/6>*<x-s<p>>;ff4=Rf<s<p>>-M<p>*h<p>*h<p>/6;f<m>=ff1+ff2+ff3+ff4;end%%%作出插值圖形x=a:hh:b;plot<x,f,'k'>holdongridon結(jié)果分析和討論：本實驗采用函數(shù)進行數(shù)值插值,插值區(qū)間為[-1,1],給定節(jié)點為xj=-1+jh,h=0.1,j=0,…,n。下面分別給出Lagrange插值,三次樣條插值,線性插值的函數(shù)曲線和數(shù)據(jù)表。圖中只標(biāo)出Lagrange插值的十次多項式的曲線,其它曲線沒有標(biāo)出,從數(shù)據(jù)表中可以看出具體的誤差。表中,L10<x>為Lagrange插值的10次多項式,S10<x>,S40<x>分別代表n=10,40的三次樣條插值函數(shù),X10<x>,X40<x>分別代表n=10,40的線性分段插值函數(shù)。xf<x>L10<x>S10<x>S40<x>X10<x>X40<x>-1.000000000000000.038461538461540.038461538461540.038461538461540.038461538461540.038461538461540.03846153846154-0.950000000000000.042440318302391.923631149719200.042408331510400.042440318302390.043552036199100.04244031830239-0.900000000000000.047058823529411.578720990349260.047096975854580.047058823529410.048642533936650.04705882352941-0.800000000000000.058823529411760.058823529411760.058823529411760.058823529411760.058823529411760.058823529411768660.075471698113210.079411764705880.07547169811321-0.650000000000000.08648648648649-0.072604203224180.085897763608490.086486486486490.089705882352940.08648648648649-0.400000000000000.200000000000000.200000000000000.200000000000000.200000000000000.200000000000000.20000000000000-0.300000000000000.307692307692310.235346591310800.297356916958600.307692307692310.350000000000000.30769230769231-0.250000000000000.390243902439020.342641234397890.380487381403270.390243902439020.425000000000000.39024390243902-0.200000000000000.500000000000000.500000000000000.500000000000000.500000000000000.500000000000000.500000000000007469693684310.640000000000000.625000000000000.64000000000000764705882401.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.0000000000000000000000000.678989577293400.657469693684310.640000000000000.625000000000000.640000000000000.200000000000000.500000000000000.500000000000000.500000000000000.500000000000000.500000000000000.500000000000000.250000000000000.390243902439020.342641234397890.380487381403270.390243902439020.425000000000000.390243902439020.300000000000000.307692307692310.235346591310800.297356916958600.307692307692310.350000000000000.3076923076923100.400000000000000.200000000000000.200000000000000.200000000000000.200000000000000.200000000000000.200000000000000.650000000000000.08648648648649-0.072604203224180.085897763608490.086486486486490.089705882352940.08648648648649113210.079411764705880.075471698113210.750000000000000.0660.800000000000000.058823529411760.058823529411760.058823529411760.058823529411760.058823529411760.058823529411760.85000000000.900000000000000.047058823529411.578720990349260.047096975854580.047058823529410.048642533936650.047058823529410.950000000000000.042440318302391.923631149719200.042408331510400.042440318302390.043552036199100.042440318302391.000000000000000.038461538461540.038461538461540.038461538461540.038461538461540.038461538461540.03846153846154從以上結(jié)果可以看到,用三次樣條插值和線性分段插值,不會出現(xiàn)多項式插值是出現(xiàn)的Runge現(xiàn)象,插值效果明顯提高。進一步說,為了提高插值精度,用三次樣條插值和線性分段插值是可以增加插值節(jié)點的辦法來滿足要求,而用多項式插值函數(shù)時,節(jié)點數(shù)的增加必然會使多項式的次數(shù)增加,這樣會引起數(shù)值不穩(wěn)定,所以說這兩種插值要比多項式插值好的多。而且在給定節(jié)點數(shù)的條件下,三次樣條插值的精度要優(yōu)于線性分段插值,曲線的光滑性也要好一些。實驗報告四題目：多項式最小二乘法摘要：對于具體實驗時,通常不是先給出函數(shù)的解析式,再進行實驗,而是通過實驗的觀察和測量給出離散的一些點,再來求出具體的函數(shù)解析式。又因為測量誤差的存在,實際真實的解析式曲線并不一定通過測量給出的所有點。最小二乘法是求解這一問題的很好的方法,本實驗運用這一方法實現(xiàn)對給定數(shù)據(jù)的擬合。前言：〔目的和意義學(xué)習(xí)使用最小二成法的原理了解法方程的特性數(shù)學(xué)原理：對于給定的測量數(shù)據(jù)<xi,fi><i=1,2,…,n>,設(shè)函數(shù)分布為特別的,取為多項式<j=0,1,…,m>則根據(jù)最小二乘法原理,可以構(gòu)造泛函令<k=0,1,…,m>則可以得到法方程求該解方程組,則可以得到解,因此可得到數(shù)據(jù)的最小二乘解程序設(shè)計：本實驗采用Matlab的M文件編寫。其中多項式函數(shù)寫成function的方式,如下functiony=fai<x,j>y=1;fori=1:jy=x.*y;end寫成如上形式即可,下面給出主程序。多項式最小二乘法源程序clear%%%給定測量數(shù)據(jù)點<s,f>s=[3456789];f=[2.012.983.505.025.476.027.05];%%%計算給定的數(shù)據(jù)點的數(shù)目n=length<f>;%%%給定需要擬合的數(shù)據(jù)的最高次多項式的次數(shù)m=10;%%%程序主體fork=0:m;g=zeros<1,m+1>;forj=0:m;t=0;fori=1:n;%計算內(nèi)積<fai<si>,fai<si>>t=t+fai<s<i>,j>*fai<s<i>,k>;endg<j+1>=t;endA<k+1,:>=g;%法方程的系數(shù)矩陣t=0;fori=1:n;%計算內(nèi)積<f<si>,fai<si>>t=t+f<i>*fai<s<i>,k>;endb<k+1,1>=t;enda=A\b%求出多項式系數(shù)x=[s<1>:0.01:s<n>]';y=0;fori=0:m;y=y+a<i+1>*fai<x,i>;endplot<x,y>%作出擬合成的多項式的曲線gridonholdonplot<s,f,'rx'>%在上圖中標(biāo)記給定的點結(jié)果分析和討論：例用最小二乘法處理下面的實驗數(shù)據(jù).xi3456789fi2.012.983.505.025.476.027.05并作出的近似分布圖。分別采用一次,二次和五次多項式來擬合數(shù)據(jù)得到相應(yīng)的擬合多項式為：y1=-0.38643+0.82750x；2；y5=-50.75309+51.53527x-19.65947x2+3.66585x3-0.32886x4+0.01137x5；分別作出它們的曲線圖,圖中點劃線為y1曲線,實線為y2曲線,虛線為y5曲線。’x’為給定的數(shù)據(jù)點。從圖中可以看出并不是多項式次數(shù)越高越好,次數(shù)高了,曲線越能給定點處和實際吻合,但別的地方就很差了。因此,本例選用一次和兩次的多項式擬合應(yīng)該就可以了。實驗報告五題目：Romberg積分法摘要：對于實際的工程積分問題,很難應(yīng)用Newton-Leibnitz公式去求解。因此應(yīng)用數(shù)值方法進行求解積分問題已經(jīng)有著很廣泛的應(yīng)用,本文基于Romberg積分法來解決一類積分問題。前言：〔目的和意義理解和掌握Romberg積分法的原理；學(xué)會使用Romberg積分法；明確Romberg積分法的收斂速度及應(yīng)用時容易出現(xiàn)的問題。數(shù)學(xué)原理：考慮積分,欲求其近似值,通常有復(fù)化的梯形公式、Simpsion公式和Cotes公式。但是給定一個精度,這些公式達到要求的速度很緩慢。如何提高收斂速度,自然是人們極為關(guān)心的課題。為此,記T1,k為將區(qū)間[a,b]進行2k等分的復(fù)化的梯形公式計算結(jié)果,記T2,k為將區(qū)間[a,b]進行2k等分的復(fù)化的Simpsion公式計算結(jié)果,記T3,k為將區(qū)間[a,b]進行2k等分的復(fù)化的Cotes公式計算結(jié)果。根據(jù)Richardson外推加速方法,可以得到收斂速度較快的Romberg積分法。其具體的計算公式為：準備初值,計算按梯形公式的遞推關(guān)系,計算按Romberg積分公式計算加速值m=2,…,k精度控制。對給定的精度,若則終止計算,并取為所求結(jié)果；否則返回2重復(fù)計算,直至滿足要求的精度為止。程序設(shè)計：本實驗采用Matlab的M文件編寫。其中待積分的函數(shù)寫成function的方式,例如如下functionyy=f<x,y>;yy=x.^3;寫成如上形式即可,下面給出主程序Romberg積分法源程序%%%Romberg積分法clear%%%積分區(qū)間b=3;a=1;%%%精度要求R=1e-5;%%%應(yīng)用梯形公式準備初值T<1,1>=<b-a>*<f<b>+f<a>>/2;T<1,2>=T<1,1>/2+<b-a>/2*f<<b+a>/2>;T<2,1>=<4*T<1,2>-T<1,1>>/<4-1>;j=2;m=2;%%%主程序體%%%while<abs<T<m,1>-T<m-1,1>>>R>;%%%精度控制j=j+1;s=0;forp=1:2^<j-2>;s=s+f<a+<2*p-1>*h/<2^<j-1>>>;endT<1,j>=T<1,j-1>/2+h*s/<2^<j-1>>;%%%梯形公式應(yīng)用form=2:j;k=<j-m+1>;T<m,k>=<<4^<m-1>>*T<m-1,k+1>-T<m-1,k>>/<4^<m-1>-1>;endend%%%給出Romberg積分法的函數(shù)表I=T<m,1>結(jié)果分析和討論：進行具體的積分時,精度取R=1e-8。求積分。精確解I=24999676。運行程序得Romberg積分法的函數(shù)表為1.0e+007*2.499967600000002.4999676000000002.4999676000000000由函數(shù)表知Romberg積分給出的結(jié)果為2.4999676*10^7,與精確沒有誤差,精度很高。求積分。精確解I=ln3=1.09861228866811。運行程序得Romberg積分法的函數(shù)表為1.111271.0986123199913001.099259259259261.098640371973711.098613022277491.098612302616251.09861228889937001.098630548366001.098612588155331.098612291193061.098612288681640001.098612517723131.098612290028501.0986122886717900001.098612289805931.09861228867046000001.09861228867019000000從積分表中可以看出程序運行的結(jié)果為1.09861228867019,取8位有效數(shù)字,滿足要求。求積分。直接按前面方法進行積分,會發(fā)現(xiàn)系統(tǒng)報錯,出現(xiàn)了0為除數(shù)的現(xiàn)象。出現(xiàn)這種情況的原因就是當(dāng)x=0時,被積函數(shù)分母出現(xiàn)了0,如果用一個適當(dāng)?shù)男?shù)〔最好不要小于程序給定的最小誤差值,但是不能小于機器的最大精度來代替,可以避免這個問題。本實驗取,可得函數(shù)表為：0.920735483196590.939793275001900.944513511714170.945690853594890.945985019937430.946082994063680.946083059350920.94608306035138000.946083060387220.946083060367260000.946083060367180000故該函數(shù)的積分為0.94608306036718,取8位有效數(shù)字。求積分本題的解析解很難給出,但運用Romberg積分可以很容易給出近似解,函數(shù)表為：0.420735492400000故該函數(shù)的積分為0,取8位有效數(shù)字。結(jié)論：Romberg積分通常要求被積函數(shù)在積分區(qū)間上沒有奇點。如有奇點,且奇點為第一間斷點,那么采用例3的方法,還是能夠求出來的,否則,必須采用其它的積分方法。當(dāng)然,Romberg積分的收斂速度還是比較快的。實驗報告六題目：常微分方程初值問題的數(shù)值解法摘要：本實驗主要采用經(jīng)典四階的R-K方法和四階Adams預(yù)測-校正方法來求解常微分方程的數(shù)值解。前言：〔目的和意義通過編寫程序,進行上機計算,使得對常微分方程初值問題的數(shù)值解法有更深刻的理解,掌握單步法和線性多步法是如何進行實際計算的及兩類方法的適用范圍和優(yōu)缺點,特別是對這兩類方法中最有代表性的方法:R-K方法和Adams方法及預(yù)測-校正方法有更好的理解。通過這兩種方法的配合使用,掌握不同的方法如何配合在一起,解決實際問題。數(shù)學(xué)原理：對于一階常微分方程初值問題〔1的數(shù)值解法是近似計算中很中的一部分。常微分方程的數(shù)值解法通常就是給出定義域上的n個等距節(jié)點,求出所對應(yīng)的函數(shù)值yn。通常其數(shù)值解法可分為兩大類：單步法：這類方法在計算yn+1的值時,只需要知道xn+1、xn和yn即可,就可算出。這類方法典型有歐拉法和R-K法。多步法：這類方法在計算yn+1的值時,除了需要知道xn+1、xn和yn值外,還需要知道前k步的值。典型的方法如Adams法。經(jīng)典的R-K法是一個四階的方法。它最大的優(yōu)點就是它是單步法,精度高,計算過程便于改變步長。其缺點也很明顯,計算量大,每前進一步就要計算四次函數(shù)值f。它的具體的計算公式如式〔2所示。四階Adams預(yù)測-校正方法是一個線性多步法,它是由Adams顯式公式〔2和隱式公式組成,其計算公式如式〔3所示。預(yù)測〔3a求導(dǎo)〔3b校正〔3c求導(dǎo)〔3d將局部截斷誤差用預(yù)測值和校正值來表示,在預(yù)測和校正的公式中分別以它們各自的階段誤差來進行彌補,可期望的到精度更高的修正的預(yù)測-校正公式為：預(yù)測修正求導(dǎo)校正修正求導(dǎo)由于開始時無預(yù)測值和校正值可以利用,故令p0=c0=0,以后按上面進行計算。此方法的優(yōu)點是可以減少計算量；缺點是它不是自開始的,需要先知道前面的四個點的值,因此不能單獨使用。另外,它也不便于改變步長。程序設(shè)計：本實驗采用Matlab的M文件編寫。其中待求的微分方程寫成function的方式,如下functionyy=g<x,y>;yy=-x*x-y*y;寫成如上形式即可,下面給出主程序。經(jīng)典四階的R-K方法源程序clear%%%步長選取h=0.1;%%%初始條件,即x=0時,y=1。y<1>=1;%%%求解區(qū)間a=0;b=2;%%%迭代公式forx=a:h:b-h;k1=g<x,y<<x-a>/h+1>>;%%,下同。k2=g<x+h/2,y<<x-a>/h+1>+h/2*k1>;k3=g<x+h/2,y<<x-a>/h+1>+h/2*k2>;k4=g<x+h,y<<x-a>/h+1>+h*k3>;y<<x-a>/h+2>=y<<x-a>/h+1>+h*<k1+2*k2+2*k3+k4>/6;end四階Adams預(yù)測-校正方法源程序%%%步長選取h=0.1;%%%初始條件y<1>=1;%%%求解區(qū)間a=0;b=2;%%%應(yīng)用RK迭代公式計算初始值y0,y1,y2,y3forx=a:h:a+2*h;k1=g<x,y<<x-a>/h+1>>;k2=g<x+h/2,y<<x-a>/h+1>+h/2*k1>;k3=g<x+h/2,y<<x-a>/h+1>+h/2*k2>;k4=g<x+h,y<<x-a>/h+1>+h*k3>;y<<x-a>/h+2>=y<<x-a>/h+1>+h*<k1+2*k2+2*k3+k4>/6;end%%%應(yīng)用預(yù)測校正法求解c<4>=0;%%%校正初值p<4>=0;%%%預(yù)測初值f<1>=g<a+0*h,y<1>>;,且將該值存在數(shù)組f中。f<2>=g<a