中北大學 數(shù)值分析1

中北大學理學院實驗報告實驗類別：數(shù)值分析專業(yè)：信息與計算科學班 級：13080241學 號：1308024120姓名：楊燕中北大學理學院實驗一函數(shù)插值方法實驗內(nèi)容】給定一元函數(shù)y二f(x)的n+1個節(jié)點值y二f(x)(j二0,1,n)，數(shù)據(jù)如下:jjxj0.40.550.650.800.951.05yj0.410750.578150.696750.901.001.25382求五次Lagrange多項式L5(x)或分段三次插值多項式或Newton插值多項式，并計f(0.596),f(0.99)的值。(提示：結(jié)果為f(0.596)~0.625732,f(0.99)~1.05423)實驗方法與步驟】利用Lagrange插值公式x利用Lagrange插值公式x-x

i

x-xI=0k,i豐kyk，用C語言編寫出插值多項式程序如下：#include<stdio.h>#defineN5floatx[]={0.4,0.55,0.65,0.80,0.95,1.05};floaty[]={0.41075,0.57815,0.69675,0.90,1.00,1.25382};floatp(floatxx){inti,k;floatpp=0,m1,m2;for(i=0;i<=N;i++){m1=1;m2=1;for(k=0;k<=N;k++)if(k!=i){m1*=xx-x[k];m2*=x[i]-x[k];}pp+=y[i]*m1/m2;}returnpp;}main(){printf("f(0.596)=%lf\n",p(0.596));printf("f(0.99)=%lf\n",p(0.99));}【實驗結(jié)果】【思考】1、給出的程序求f(1.06)行不行，精度高不高？2、五次Lagrange多項式與Newton插值多項式是同一個多項式嗎？五次Lagrange多項式與Newton插值多項式是同一個多項式。3、為什么高次插值不能令人滿意？一般來說，節(jié)點個數(shù)越多，插值函數(shù)和被插值函數(shù)就有越多的地方相等。但是隨著插值節(jié)點個數(shù)的增加，兩個插值節(jié)點之間插值函數(shù)并不一定能夠很好地逼近被插值函數(shù)。再次，從舍入誤差看，高次插值由于計算量大，可能會產(chǎn)生更嚴重的誤差積累，所以，穩(wěn)定性得不到保證。此時就會出現(xiàn)龍格現(xiàn)象。中北大學理學院實驗報告實驗類別：數(shù)值分析專業(yè)：信息與計算科學班 級：13080241學 號：1308024120姓名：楊燕中北大學理學院實驗二函數(shù)逼近與曲線擬合實驗內(nèi)容】1、編寫出Legendre、Chebyshev多項式的程序；2、從隨機的數(shù)據(jù)中找出其規(guī)律性，給出其近似表達式，在生產(chǎn)實踐和科學實驗中大量存在，通常利用數(shù)據(jù)的最小二乘法求得擬合曲線。例如在某冶煉過程中根據(jù)統(tǒng)計數(shù)據(jù)的含碳量與時間關(guān)系，試求含碳量y與時間t的擬合曲線。（t分）0510152025303540455055yC10一4)01.72.162.863.443.874.154.374.514.584.024.64i111111111111實驗方法與步驟】1、用C編寫語言出Legendre多項式的程序如下:#include<stdio.h>doublep(intn,doublex){if(n==0)return1;elseif(n==1)returnx;elsereturn((2*n-1)*x*p(n-1,x)-(n-1)*p(n-2,x))/n;}intmain(){intn;doublex;doubley;printf("inputn,x:\n");scanf("%d%lf",&n,&x);y=p(n,x);printf("%lf\n",y);return0;}2、給出表格數(shù)據(jù)的近似解析表達式為；申（t）=a+bt，選取基函數(shù)申o（x）二1,申1（x）=x，則得到(p,p)=12，(甲,*)=(*,*)=丈x,(甲,甲)=藝x2，(cp,f)=藝f，G,f)=￡fx‘001001i11i0i1iii=0 i=0 i=0 i=0于是得方程組険P0)a+%P1)b=%f),用C語言編寫曲線擬合的程序如I(p,p丿a+(p,p)b=(p,f)10111下：#include<stdio.h>#defineMax_N25main(){inti,n;doublex[Max_N],y[Max_N];doublea11,a12,a21,a22,d1,d2;doublea,b;printf("\nInputnvalue:");do{scanf("%d",&n);if(n>Max_N)printf("\npleasere_inputnvalue:");}while(n>Max_N||n<=0);printf("inputx[i],i=0,...%d:\n",n-1);for(i=0;i<n;i++)scanf("%lf",&x[i]);printf("inputy[i],i=0,...%d:\n",n-1);for(i=0;i<n;i++)scanf("%lf",&y[i]);for(i=0;i<n;i++){a21+=x[i];a22+=x[i]*x[i];d1+=y[i];d2+=x[i]*y[i];}a12=a21;a11=n;a=(d1*a22-d2*a12)/(a11*a22-a12*a21);b=(d1*a21-d2*a11)/(a21*a12-a22*a11);printf("slove:P(x)=%f+%fx\n",a,b);getchar();return0;}3、為作比較，用MATLAB進行曲線擬合，編寫的程序如下:x=[0510152025303540455055];y=[01.272.162.863.443.874.154.374.514.584.024.64];p=polyfit(x,y,1);disp([num2str(p(1)),'*x+',num2str(p(2))]);xx=linspace(0,60,60);yy=polyval(p,xx);plot(x,y,'rx',xx,yy)【實驗結(jié)果】1、C語言編寫的Legendre多項式的程序結(jié)果如下inputn,x：43321.000000Pressanykeytocontinue2、C語言擬合的曲線結(jié)果如下：Inputnualue-12inputxti1,1=0,...11：0ii1315翎2S303540455055inputyti.],1=0,...11：Q1.272.162.863.443.874.154.374.514.584.924.64Sloue：P<x>=l.321539+0.672762x3、MATLAB擬合的曲線結(jié)果和誤差分析結(jié)果如下:擬合的方程為：y=0.072762x+1.3215【思考】1、 最佳一致逼近與最佳平方逼近的區(qū)別是什么？最佳一致逼近中的范數(shù)取的是無窮大范數(shù)，公式為||f(x)-P*(x) =minmaxf(x)-P(x)。gP^Hna<x<b最佳平方逼近中的范數(shù)取的是二范數(shù)，公式為||f(x)—P*(x)|2=minJb[f(x)—P(x)1dx。2PwHan2、 曲線擬合與最佳平方逼近的區(qū)別是什么？曲線擬合中的f(x)是[a,b]上的一個列表函數(shù)，公式為：||f-P*||2=min￡[f(x)-P(x)1。2唄i=o i i最佳平方逼近中的f(x)是一個未知函數(shù)，公式為：f(x-P*x(2=) Jmfixn-Px(I2dx。()2PeHan3、能用什么方法確定最小二乘法的擬合函數(shù)？擬合的方法除了最小二乘法外，還有拉格朗日插值法、牛頓插值法、區(qū)間二分法、雅克比迭代法和牛頓科特斯數(shù)值積分法等，都可以確定擬合函數(shù)。中北大學理學院實驗報告實驗類別：數(shù)值分析專業(yè)：信息與計算科學班 級：13080241學 號：1308024120姓名：楊燕中北大學理學院實驗三數(shù)值積分與數(shù)值微分實驗內(nèi)容】選用復合梯形公式，復合Simpson公式，Romberg算法、高斯算法計算(1)I二宀小-sin2xdx(2)I=J1 ―dx(3)I=11旦匕衛(wèi)dXo o4+x2 o1+x2實驗方法與步驟】1、用C語言編寫Romberg算法數(shù)值積分的程序如下：#include<stdio.h>#include<math.h>#defineepsilono.ooo1/*求基函數(shù)*/floatf(floatx){return(sqrt(4-sin(x)*sin(x)));} /*梯形公式*/floatRomberg(floata,floatb){intk=1;floatS,x,T1,T2,S1,S2,C1,C2,R1,R2,h=b-a;T1=h/2*(f(a)+f(b));while(1){S=o;x=a+h/2;do{S+=f(x);x+=h;}while(x<b);T2=(T1+h*S)/2.o;if(fabs(T2-T1)<epsilon)return(T2);S2=T2+(T2-T1)/3.0;if(k==1){T1=T2;S1=S2;h/=2;k+=1;continue;}if(fabs(S2-S1)<epsilon)return(S2);C2=S2+(S2-S1)/15.0;if(k==2){C1=C2;T1=T2;S1=S2;h/=2;k+=1;continue;}if(fabs(C2-C1)<epsilon)return(C2);R2=C2+(C2-C1)/63.0;if(k==3){R1=R2;C1=C2;T1=T2;S1=S2;h/=2;k+=1;continue;}if(fabs(R2-R1)<epsilon)return(R2);R1=R2;C1=C2;T1=T2;S1=S2;h/=2;k+=1;}}main(){inti;floata,b,S;printf("\nInputbeginandend:");scanf("%f%f",&a,&b);S=Romberg(a,b);printf("Solveis:%f",S);scanf("%",S);}2、再用復合Simpson公式計算同一個積分，并比較其結(jié)果，且用C語言編寫的程序如下：#include<stdio.h>#include<math.h>floatf(floatx){return(sqrt(4-sin(x)*sin(x)));}floatSimpson(floata,floatb){intk=1,n=20;floats,x1,x2,h=(b-a)/n,c=a+h/2;s=h/6*(f(a)+f(b)+4*f(c));for(k=1;k<=n-1;k++){x1=a+k*h;x2=x1+h/2;s=s+(2*h/3)*f(x2)+(h/3)*f(x1);

return(s);main(){inti;floata,b,s;printf(“請輸入a,b:");scanf("%f%f",&a,&b);s=Simpson(a,b);printf("解為:%f",s);scanf("%",s);}3、分別取不同步長h二(b-a)/n，試比較計算結(jié)果(如n二10,20等);4、給定精度要求e，試用變步長算法，確定最佳步長。實驗結(jié)果】1、所給三個定積分的結(jié)果如下：Inputbeginandend：00.25Solueis=0.498711Inputbeginandend：01Bolueis：0.2721972、對第一個定積分I=I‘4J4-sin2xdx用Romberg算法和復合Simpson公式計0算的結(jié)果分別如下：=00-25解為：0?498711Inputbeginandend=00-25解為：0?4987113、步長為h二(b-a)/10時，所給三個定積分的結(jié)果如下:Inputbeginandend：00.25Solueis：0.498708Inputbeginandend：01Solveis：0.2721784、步長為h二(b-a)/20時，所給三個定積分的結(jié)果如下:Inputbeginandend：Inputbeginandend：00.25Solueis=0.498520Inputbeginandend：01olueis：0.390654Inputbeginandend：01Solueis=0.272082【思考】1、 復合求積公式的優(yōu)點是什么？復合求積公式通過把區(qū)間分成若干個子區(qū)間，再在每個子區(qū)間上用低階求積公式來提高求積的精度。2、 Romberg算法與與牛頓-柯特斯求積算法的聯(lián)系是什么？Romberg算法與與牛頓-柯特斯求積算法的節(jié)點都是等距選取均勻分布的，且n>8的牛頓-柯特斯公式因為誤差太大而不能使用。3、 高斯算法的求積節(jié)點如何確定？以這些節(jié)點為零的多項式①(x)二(x-x)(x-x)…(x-x)與任何次數(shù)不超過TOC\o"1-5"\h\zn+1 0 1 nn的多項式p(x)帶權(quán)p(x)正交，即：Ibp(x)e(x)p(x)dx=0。a n+14、 牛頓-柯特斯求積與高斯算法的節(jié)點分布有什么不同？牛頓-柯特斯求積的節(jié)點是等距選取均勻分布的，而高斯算法的節(jié)點是滿足以這些節(jié)點為零的多項式3(x)二(x-x)(x-x)???(x-x)與任何次數(shù)不超過n的多n+1 0 1 n項式p(x)帶權(quán)p(x)正交，即是不均勻分布的。中北大學理學院實驗報告實驗類別：數(shù)值分析專業(yè)：信息與計算科學班 級：13080241學 號：1308024120姓名：楊燕中北大學理學院

實驗四線性方程組的直接解法實驗內(nèi)容】2x一x=312用追趕法）一x+3x一2x=1用追趕法）1 2 3一2x+4x一2x=0234一2x+5x=一534實驗方法與步驟】1、對上述方程組用追趕法求解，用C語言編程如下：#include<stdio.h>#include<math.h>#include<string.h>#defineN5main(){floata[N]={0,0,-1,-2,-2};floatb[N]={0,2,3,4,5};floatc[N]={0,-1,-2,-2,0};floatd[N]={0,3,1,0,-5};floatx[N]={0,0,0,0};floatr[N]={0,0,0,0};floaty[N]={0,0,0,0};floatq;intk;r[1]=c[1]/b[1];y[1]=d[1]/b[1];for(k=2;k<=N-1;k++){q=b[k]-r[k-1]*a[k];r[k]=c[k]/q;y[k]=(d[k]-y[k-1]*a[k])/q;}y[N-1]=(d[N-1]-y[N-2]*a[N-1])/(b[N-1]-r[N-2]*a[N-1]);x[N-1]=y[N-1];for(k=N-2;k>=1;k--)x[k]=y[k]-r[k]*x[k+1];for(k=1;k<N;k++)printf("x[%d]=%lf\n",