數(shù)值分析實(shí)習(xí)報(bào)告

計(jì)算機(jī)科學(xué)與技術(shù)CS0809班陶松橋數(shù)值分析實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)要求:實(shí)驗(yàn)一：Lagrange插值對[-5，5]作等距分劃xi=-5+ih，h=10/n，i=0，1，…，n，并對Runge給出的函數(shù)y=作Lagrange插值，取n=10，20計(jì)算插值多項(xiàng)式Pn(x)在x=4.8處的誤差，并作分析。實(shí)驗(yàn)二：Romberg算法用Romberg外推方法計(jì)算圓周率π=0141+實(shí)驗(yàn)三：經(jīng)典的四階Runge-Kutta法用經(jīng)典的四階Runge-Kutta法解初值問題y求y(2)的數(shù)值解，取h=0.1,0.05,0.01,…計(jì)算，并分析其收斂性.(精確解yx公式與算法描述:實(shí)驗(yàn)一：已知xx0x1x2…xnyy0y1y2…ynLagrange插值公式：P算法如下：輸入（xi,yi）i=0,1,2,…,n;x輸出yStep1.y=0；Step2.Fori=0tonDo{L=1.0；Forj=0tonDo{if(j≠i)thenL=L(x-xj)/(xi–xj)；}y=y+yiL；}Step3.Print(x,y)實(shí)驗(yàn)二：Romberg公式I=abfxI≈TnI-TnT2n=12Sn=43算法如下：輸入：a，b，ε，f(x)；輸出：IStep1.置k=0；精度要求ε；h=b-a；T1Step2.h=h/2;k=k+1;T2=TStep3.h=h/2;k=k+1;T4S2=Step4.h=h/2;k=k+1;T8S4=43Step5.h=h/2;k=k+1;T2S2k-1=4Step6.if|R2k-3elsegotoStep5實(shí)驗(yàn)三：經(jīng)典的四階Runge-Kutta公式：yn+1=ykkkk算法如下：輸入:f(x,y)，a，b，y0，h；輸出:yStep1:輸入a,b,y0,h；Step2:(b-a)/h=>N,a=>x,y0=>yStep3:ForI=1T0Nhf(x,y)=>k1hf(x+h/2,y+k1/2)=>k2hf(x+h/2,y+k2/2)=>k3hf(x+h,y＋k3)=>k4y+(k1+2k2+2k3+k4)/6=>yx+h=>xStep4:輸出(x,y)源程序?qū)嶒?yàn)一:lagrange.m文件如下：%lagrange.mfunctiony=lagrange(x,n)%x為所求插值結(jié)點(diǎn)y返回插值結(jié)果n為差值多項(xiàng)式最高項(xiàng)次數(shù)formatlongs=0;n=n+1;xx=linspace(-5,5,n);yy=1./(1+xx.^2);fori=1:nt=ones(1,length(x));forj=1:nifj~=it=t.*(x-xx(j))/(xx(i)-xx(j));endends=s+t*yy(i);endy=s;實(shí)驗(yàn)二：romberg.m文件如下：functionI=romberg(f,a,b,e)%a是積分下限b是積分上限e是要求誤差n=4;%n的初值為4formatlongwhile(1)T=zeros(n+1,4);%構(gòu)造（n+1）*4的矩陣TT(1,1)=(b-a)*(f(a)+f(b))/2;%f為被積函數(shù)fork=1:nh=(b-a)/(2^k);m=2^(k-1);s=0;fori=1:ms=s+f(a+h*(2*i-1));endT(k+1,1)=T(k,1)/2+h*s;%算出矩陣T的第一列endforp=1:3forj=p:nT(j+1,p+1)=((4^p)*T(j+1,p)-T(j,p))/(4^p-1);endendifabs(T(n+1,4)-T(n,4))<e%如果差值小于誤差，跳出，結(jié)束。否則進(jìn)入下次循環(huán)break;endn=n+1;endI=T(n+1,4);%輸出T（n+1,4）的值，即為所要求的值實(shí)驗(yàn)三：runge_kuta.m文件如下:functiony=runge_kutta(f,a,c,xx,h)%f表示y’,已知y(a)=c,求y(xx),h表示步長formatlongx0=a;y0=c;x=x0:h:xx;n=length(x);yt=y0;fori=0:n-2xt=x0+i*h;k1=h*f(xt,yt);k2=h*f(xt+h/2,yt+k1/2);k3=h*f(xt+h/2,yt+k2/2);k4=h*f(xt+h,yt+k3);ytt=yt+1/6*(k1+2*k2+2*k3+k4);yt=ytt;endy=yt;四、開發(fā)與運(yùn)算環(huán)境:本程序在matlabR2009b中編輯與調(diào)試,運(yùn)算環(huán)境為windows7。五、計(jì)算結(jié)果及分析:實(shí)驗(yàn)一:操作及結(jié)果如下：>>clear>>y=lagrange(4.8,10)R=-1.762788118357620y=1.804385456128002>>y=lagrange(4.8,20)R=50.906012520135377y=-50.864415182364993用以下命令做出插值圖形進(jìn)行比較:>>symsx>>y0=1/(1+x^2);>>y1=lagrange(x,10);>>y2=lagrange(x,20);>>ezplot(y2,[-5,5])>>holdon>>ezplot(y1,[-5,5])>>ezplot(y0,[-5,5])>>holdoff對圖形稍加修改結(jié)果如下:由圖可知,插值多項(xiàng)式Pn(10)與Pn(20)與原函數(shù)出現(xiàn)了較大偏差，尤其是在附近，龍格現(xiàn)象非常明顯。在x=4.8處誤差分別為-1.76和50.906。因此此函數(shù)不適合用高次的拉格朗日插值法進(jìn)行插值計(jì)算。實(shí)驗(yàn)二：操作及結(jié)果如下：>>clear>>f=inline('4/(1+x^2)')f=Inlinefunction:f(x)=4/(1+x^2)>>e=1/2*10^(-8);>>I=romberg(f,0,1,e)I=3.141592653589794計(jì)算結(jié)果表明：用Romberg外推方法應(yīng)用π=實(shí)驗(yàn)三：操作及結(jié)果如下:>>clear>>f=inline('x*y^(1/3)','x','y')f=Inlinefunction:f(x,y)=x*y^(1/3)>>y=runge_kutta(f,1,1,2,0.1)y=2.828426783860100>>y=runge_kutta(f,1,1,2,0.05)y=2.828427103015112>>y=runge_kutta(f,1,1,2,0.01)y=2.828427124710866標(biāo)準(zhǔn)四階龍格-庫塔公式的局部截?cái)嗾`差為O（h5）,其增量函數(shù)φx=16h[k1x,y,h+2取步h=0.1，0.05，0.01都達(dá)到了很好的計(jì)算結(jié)果。六、體會與問題:通過本實(shí)驗(yàn)，我很好的掌握了拉格朗插值、龍貝格求積和標(biāo)準(zhǔn)四階龍格-庫塔公式的算法?？吹皆搶?shí)驗(yàn)，我第一個(gè)想到的計(jì)算機(jī)高級語言就是Matlab語言，一方面是由于matlab語言的簡單性，另一方面還因?yàn)槠鋸?qiáng)大的計(jì)算能力和它的靈活性。在實(shí)現(xiàn)的過程中，可以方便的通過圖形的對比，數(shù)值與字符的轉(zhuǎn)化以及計(jì)算精度的轉(zhuǎn)化來很好的了解和分析計(jì)算結(jié)果。另外，程序還具有普遍性，就拿Lagrange插值程序來說，由于傳遞參數(shù)未設(shè)置類型，當(dāng)給的參數(shù)是一個(gè)數(shù)時(shí)，程序能計(jì)算出該值對應(yīng)的結(jié)果，當(dāng)接收到的參數(shù)值是一個(gè)向量即一組數(shù)時(shí)，程序能一次性上的計(jì)算出各值對應(yīng)的近似函數(shù)值，當(dāng)給的參數(shù)是字符形時(shí)，程序便能計(jì)算出插值公式及余項(xiàng)公式，這是其他計(jì)算機(jī)語言較