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