計算機科學與技術CS0809班陶松橋數(shù)值分析實驗報告實驗要求:實驗一：Lagrange插值對[-5，5]作等距分劃xi=-5+ih，h=10/n，i=0，1，…，n，并對Runge給出的函數(shù)y=作Lagrange插值，取n=10，20計算插值多項式Pn(x)在x=4.8處的誤差，并作分析。實驗二：Romberg算法用Romberg外推方法計算圓周率π=0141+實驗三：經典的四階Runge-Kutta法用經典的四階Runge-Kutta法解初值問題y求y(2)的數(shù)值解，取h=0.1,0.05,0.01,…計算，并分析其收斂性.(精確解yx公式與算法描述:實驗一：已知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)實驗二：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實驗三：經典的四階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)源程序實驗一:lagrange.m文件如下：%lagrange.mfunctiony=lagrange(x,n)%x為所求插值結點y返回插值結果n為差值多項式最高項次數(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;實驗二：romberg.m文件如下：functionI=romberg(f,a,b,e)%a是積分下限b是積分上限e是要求誤差n=4;%n的初值為4formatlongwhile(1)T=zeros(n+1,4);%構造（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%如果差值小于誤差，跳出，結束。否則進入下次循環(huán)break;endn=n+1;endI=T(n+1,4);%輸出T（n+1,4）的值，即為所要求的值實驗三：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ā)與運算環(huán)境:本程序在matlabR2009b中編輯與調試,運算環(huán)境為windows7。五、計算結果及分析:實驗一:操作及結果如下：>>clear>>y=lagrange(4.8,10)R=-1.762788118357620y=1.804385456128002>>y=lagrange(4.8,20)R=50.906012520135377y=-50.864415182364993用以下命令做出插值圖形進行比較:>>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對圖形稍加修改結果如下:由圖可知,插值多項式Pn(10)與Pn(20)與原函數(shù)出現(xiàn)了較大偏差，尤其是在附近，龍格現(xiàn)象非常明顯。在x=4.8處誤差分別為-1.76和50.906。因此此函數(shù)不適合用高次的拉格朗日插值法進行插值計算。實驗二：操作及結果如下：>>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計算結果表明：用Romberg外推方法應用π=實驗三：操作及結果如下:>>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標準四階龍格-庫塔公式的局部截斷誤差為O（h5）,其增量函數(shù)φx=16h[k1x,y,h+2取步h=0.1，0.05，0.01都達到了很好的計算結果。六、體會與問題:通過本實驗，我很好的掌握了拉格朗插值、龍貝格求積和標準四階龍格-庫塔公式的算法?？吹皆搶嶒?，我第一個想到的計算機高級語言就是Matlab語言，一方面是由于matlab語言的簡單性，另一方面還因為其強大的計算能力和它的靈活性。在實現(xiàn)的過程中，可以方便的通過圖形的對比，數(shù)值與字符的轉化以及計算精度的轉化來很好的了解和分析計算結果。另外，程序還具有普遍性，就拿Lagrange插值程序來說，由于傳遞參數(shù)未設置類型，當給的參數(shù)是一個數(shù)時，程序能計算出該值對應的結果，當接收到的參數(shù)值是一個向量即一組數(shù)時，程序能一次性上的計算出各值對應的近似函數(shù)值，當給的參數(shù)是字符形時，程序便能計算出插值公式及余項公式，這是其他計算機語言較