2023年計(jì)算機(jī)數(shù)值方法實(shí)驗(yàn)報(bào)告

TAIYUANUNIVERSITYOFTECHNOLOGY本科實(shí)驗(yàn)報(bào)告課程名稱：計(jì)算機(jī)數(shù)值方法實(shí)9僉項(xiàng)目：方程求根線性方程組的直接解法線性方程組的迭代解法

代數(shù)插值和最小二乘法擬合多項(xiàng)式實(shí)驗(yàn)地點(diǎn)：逸夫302專業(yè)班級(jí)：軟件學(xué)號(hào):學(xué)生姓名：指導(dǎo)教師:田華2023年4月24日

實(shí)驗(yàn)結(jié)果和分析12314012824113b[01=1.000000b[11=2.00000011高斯消元法1高斯消元法b[2]=3.1高斯消元法請(qǐng)輸入矩陣階數(shù)n：蓍輸入矩陣各項(xiàng)：42158721048361261120請(qǐng)輸入方程組的常數(shù)項(xiàng):2LU分解2773x[l]=-1.00xL2J=1.00xt3]=-1.00x[4]=1.0032LU分解21-7121-5121-5121-512-5H0]=-4.166667Hl]=l.333333i[2]=-3.5000001C3]=0.66666?U4]=-2.833333,=4/~r=上科11H實(shí)臉?lè)治觯焊咚瓜?是先消元,再回帶的過(guò)程。由程序段可以發(fā)現(xiàn),始終消去對(duì)角線下方的元素。從消元過(guò)程可以看出，對(duì)于n階線性方程組，只要各步主元素不為零，通過(guò)n-1步消元，就可以得到一個(gè)等價(jià)的系數(shù)矩陣為上三角形陣的方程組，然后再運(yùn)用回代過(guò)程可求得原方程組的解。LU分解法，分解矩陣為單位下三角陣L與上三角陣U的乘積，然后解方程組Ly=b,回代，解方程組Ux=y。其中的L為n階單位下三角陣、U為上三角陣.對(duì)于追趕法，追趕法是合用于三角矩陣的線性方程組的求解的方法，并不合用于其他類型矩陣。心得體會(huì)（碰到的問(wèn)題和解決方法）本證實(shí)驗(yàn)難度比較式在編譯時(shí)經(jīng)常出現(xiàn)各種錯(cuò)誤「程序代碼也比較繁瑣「諫深感覺(jué)到自己的上機(jī)操作能力有限，應(yīng)加強(qiáng)自己的編程能力，以后?要繼續(xù)努力。

實(shí)驗(yàn)名稱實(shí)驗(yàn)三線性方程組的迭代求解實(shí)驗(yàn)?zāi)康暮鸵?guī)定使用雅可比迭代法或高斯-賽德?tīng)柕▽?duì)下列方程組進(jìn)行求解。jlO不一“一2巧=72(一巧+10巧-2巧=83°［一天一巧+54=42實(shí)驗(yàn)內(nèi)容設(shè)線性方程組Ax=b的系數(shù)矩陣A可逆,且主對(duì)角元素aIt,a22,???，a而均不為零,令D=diag(an,a22,?,,,ann)并將A分解成A=(A-D)+D從而線性方程組可寫成Dx=(D-A)x+b則有迭代公式x(k+l)=B1x(k)+fl其中，B1=I-D-1A,fl=D-lbo重要儀器設(shè)備臺(tái)式或筆記本計(jì)算機(jī)實(shí)驗(yàn)記錄(寫出實(shí)驗(yàn)內(nèi)容中的程序代碼和運(yùn)營(yíng)結(jié)果)(可分欄或加頁(yè))雅可比迭代法#include<stdio.h>#inc1ude<math.h>ttdefinen3voidf(float*b,f1oatx[1){f1oaty[n+1]={0,0,0,1}；inti,j,k;do{k=0;for(i=0；i<n+1;i++)x[i]=y[i]；for(i=0；i<n;i++){yLi]=0;for(j=0;j<n+1;j++)y[i]+=*(b+(n+D*i+j)*x[j];}for(i=0;i<n;i++)if(fabs(y[i]-x[i])<=0.5e-3)k++;if(k==3)break}while(1);for(i=0;i<n;i++)printf(〃y[%d]=%f\n〃,i,y[i]);}main(){floatb[n][n+l]={0,0.1,0.2,0.72,0.1,0,0.2,0.83,0.2,0.2,0,0.84)；floatx[n+l]={0,0,0,1};f(b[0],x);)高斯一賽德?tīng)柕?includeiostream”#include"iomanip”usingnamespacestd;intmain(){intiJ,k=(),m,n;doublee1,e2=0.0；coutvv”請(qǐng)輸入精度e：H；，cin?el;coutvv”請(qǐng)輸入系數(shù)矩陣行數(shù)：:dcin?m;coutvv”請(qǐng)輸入系數(shù)矩陣列數(shù)：”；in?n；cout<<endl；doubIe(**a)=newdouble*[m]；〃生成二維動(dòng)態(tài)數(shù)組or(i=0；i<=m;i++)(ga[i]=newdouble[n];)，double(*b)=newdouble[m]；ddouble(*x)=newdouble[nJ;cout?"請(qǐng)輸入系數(shù)矩陣:"vVendl;cout?”——--"————一—?…n?endl;for(intnum1=0;numl<m;numl++)。for(intnum2=0;num2<n；num2++)(g，cin>>a[num1][num2];gcoutvvend1;cout<<”輸入的系數(shù)矩陣為了VVendl；dfor(intnum3=0;num3<m；num3++){gfor(intnum4=0;num4<n;num4++)(??dcout?a[num3][num4]?H";?cout<<endl;cout-------------------M?endl;ocoutvv”請(qǐng)輸入矩陣b：n<<endl;?coutVv”——■-■?-■?■?■-...—....一——-一"VVendl;for(intnum5=0;num5<m;num5++){cin>>bEnum5]；)cout?”輸入的矩陣b為:"VVend1；dfor(intnum6=0；num6<m;num6++)(cout<<b[num6]?nM;)cout?endl;cout<<H————----———--11VVend);for(intnum7=0;num7<n;num7++)，(，x[num7]=0.0000;，do(，coUtvv"第"VvkVv"次迭代值:“；de2=0.0;for(i=0;i<m；i++)°{doublesum=0.0;6forQ=0；j<n;j++)(if(j!=i)sum+=a[i][j]*x[j]；d}，11=xfi];g“2=e2；^x[i]=(b[i]-sum)/a[i][i]；^e2=(x[i])-tl>=0?(x[i])-11：tl-(x[i]);?de2=(e2>=t2?e2：t2);cout?setprecision(8)?x[i]?H}cout?endl;，k++；}whi1e(e2>=e1&&k<30)；虱outVV”共迭代了“VvkvV"次'delete[]a;?delete[]b;delete[]x;?return0;}實(shí)驗(yàn)結(jié)果和分析1雅克比迭代k?[0]=l.099811&[U=1.199811&[23=1.299777信按任意鍵繼續(xù)一.2高斯一一賽德?tīng)柕?qǐng)輸入系數(shù)矩陣:10-1-2-110-2-1-15請(qǐng)輸入矩陣b,7.28.34.2值值值值值值值值fi:值值值次KK-弋代弋代弋M-?8啾欹■欲漱欹滋的第電r￥2ly.3l￥4ly5說(shuō)曾士值值值值fi:值值值次KK-弋代弋代弋值值值值fi:值值值次KK-弋代弋代弋M-?8啾欹■欲漱欹滋的第電r￥2ly.3l￥4ly5說(shuō)曾士0.721.043081.09312951.09912651.09989041.09998621.09999831.09999980.9021.1671881.19572371.19946681.19993281.19999151.19999891.1999999?16441.28205361.29777061.29971871.29996461.29999561.29999941.2999999心得體會(huì)(碰到的問(wèn)題和解決方法)本次實(shí)驗(yàn)，讓我對(duì)這兩種方法更加理解，在編程操作上也更加純熟，此后繼續(xù)努力，不斷豐富自己的知識(shí)，增強(qiáng)操作能力。實(shí)驗(yàn)名稱實(shí)驗(yàn)四代數(shù)插值和最小二乘法擬合實(shí)驗(yàn)?zāi)康暮鸵?guī)定1使用拉格朗日插值法或牛頓插值法求解:已知f(X)在6個(gè)點(diǎn)的函數(shù)值如下表所示，運(yùn)用插值方法，求f(().596)的近似值。X0.400.550.650.800.901.05f(X)0.410.5780.696().888111.02651.25386學(xué)生姓名實(shí)驗(yàn)成績(jī)實(shí)驗(yàn)名稱實(shí)驗(yàn)一方程求根學(xué)生姓名實(shí)驗(yàn)?zāi)康暮鸵?guī)定(必填)熟悉使用、迭代法、牛頓法、割線法等方法對(duì)給定的方程進(jìn)行根的求解。選擇上述方法中的兩種方法求方程:二分法f(x)=x3+4x2-10=0在[1,2]內(nèi)的一個(gè)實(shí)根，且規(guī)定滿足精度IX*-Xn|<0.5X1O-代碼1.二分法:#include<stdio.h>#indude<std1ib.h>#inc1ude<math.h>ntmain()(doubIea=1.0,b=2.0;doubiex,s;whi1e(1)(代碼1.二分法:#include<stdio.h>#indude<std1ib.h>#inc1ude<math.h>ntmain()(doubIea=1.0,b=2.0;doubiex,s;whi1e(1)(x=(a+b)/2；=pow(x,3)+4*x*x-10;if(-0.000005<s&&s<0.000005)(break;)eIseif(s<0){a=x；)函數(shù)f(x)在區(qū)間(x,y)上連續(xù)，先在區(qū)間(x,y)擬定a與b,若f(a),f(b)異號(hào),說(shuō)明在區(qū)間(a,b)內(nèi)存在零點(diǎn)，然后求f[(a+b)/2]。假設(shè)F(a)V0,F(b)>0,a<b,①假如f[(a+b)/2]=0,該點(diǎn)即為零點(diǎn)；②假如f[(a+b)/2]<0,則區(qū)間((a+b)/2,b)內(nèi)存在零點(diǎn)，(a+b)/2Na；③假如f[(a+b)/2]>0,則區(qū)間(a,(a+b)/2)內(nèi)存在零點(diǎn),(a+b)/2Wb；返回①重新循環(huán)，不斷接近零點(diǎn)。通過(guò)每次把f(x)的零點(diǎn)所在區(qū)間收縮一半的方法，使區(qū)間內(nèi)的兩個(gè)端點(diǎn)逐步逼近函數(shù)零點(diǎn)，最終求得零點(diǎn)近似值。重要儀器設(shè)備臺(tái)式或筆記本計(jì)算機(jī)實(shí)驗(yàn)記錄(寫出實(shí)驗(yàn)內(nèi)容中的程序代碼和運(yùn)營(yíng)結(jié)果)(可分欄或加頁(yè))075157522給定數(shù)據(jù)點(diǎn)(H,y),用最小二乘法擬合數(shù)據(jù)的多項(xiàng)式，并求平方誤差。Xi0().50.60.70.80.91.011.751.962.192.442.713.00實(shí)驗(yàn)內(nèi)容1設(shè)函數(shù)在區(qū)間［a,b］上n+1互異節(jié)點(diǎn)Xo,X1,???,Xn上的函數(shù)值分別為yo,yi,???，y2求n次插值多項(xiàng)式p“(x),滿足條件Pn(Xj)=yj,j=0,1,?-?,n令Ln(x)=yolo(x)+yL(x)+--+ynln(x)=XYili(x)其中10(X),11(x),ln(x)為以X0,Xi,???,Xn為節(jié)點(diǎn)的n次插值基函數(shù)，則Ln(x)是一次數(shù)不超過(guò)n的多項(xiàng)式,且滿足Ln(xj)=Yj,L=0,1,…，n再由插值多項(xiàng)式的唯一性，得Pn(X)=Ln(x)2建立正規(guī)方程組：￡(EXij'k)ak=LXijyi,j=0,1,-,n平方誤差：I=E(Sakxik-yi)2對(duì)給定數(shù)據(jù)點(diǎn){(Xi,Yi)}(i=0,l,…，m),在取定的函數(shù)類①中，求p(x)￡①，使誤差的平方和E2最小,E?=ZEp(Xi)-Yi]#inc1ude<stdio.h>#inc1ude<stdio.h>#include<stdlib.h>#include<conio.h>#include<a11oc.h>voiddifference(float*x,float*y>intn)(ofloat*f;intk,i；f=(f1oat*)malloc(n*sizeof(float));for(k=l;k<=n;k++)(0]=y[kl;for(i=0;i<k;i+-F)of[i+1]=(ni]-y[i])/(X[k]-x[i]);3y[k]=f[k];葉dreturn;}intmain()(inti,n；floatx[20],y[20],xx,yy;printf("請(qǐng)輸入數(shù)據(jù)個(gè)數(shù)n：”)；從幾何意義上講，就是尋求與給定點(diǎn){(Xi,Yi)}(i=0,1,…,m)的距離平方和為最小的曲線y=p(x)。函數(shù)p(x)稱為擬合函數(shù)或最小二乘解，求擬合函數(shù)p(x)的方法稱為曲線擬合的最小二乘法。得到的兩個(gè)關(guān)于a0、al為未知數(shù)的兩個(gè)方程組，解這兩個(gè)方程組得出：a0=(EYi)/m—al(LXi)/ma1=[mZXiYi-(SXiZYi)]/[mEXi2-(EXi)2)]即最終的擬合多項(xiàng)式各項(xiàng)系數(shù)重要儀器設(shè)備臺(tái)式或筆記本計(jì)算機(jī)實(shí)驗(yàn)記錄(寫出實(shí)驗(yàn)內(nèi)容中的程序代碼和運(yùn)營(yíng)結(jié)果)(可分欄或加頁(yè))oscanf("%(1"，&n)；printf(”\n");for(i=0;i<=n-l;i++)(，printf("x[%d]=",i);scanfC'%fn,&x[i]);，printfC'y[%dl=",i)；aoprintf("\nM);)叩rintf(”\nM);difference(x,(float*)y,n);printf("請(qǐng)輸入插值xx:");oscanf，％"&xx);yy=yl20]；for(i=n-1；i>=0;i--)yy=yy*(xx-x[i])+y[i]；printf(H\n近似值為：(%f)=%Ann,xx,yy);)2#include<iostrearn.h>#include<fstream.h>#defineN15doublepower(double&a,intn){doubleb=l;?for(inti=0;i<n;i+4-)?b*=a；returnb;}voidGaussQ;doubleX[N],Y[N]^umX[N],sumY[N],a[N][N],b[N],1[N][N],x[N]?intmain()(ofstrearnoutdata；ifstreamindata；doub1es；?intij,kindex;coutvv”請(qǐng)輸入已知點(diǎn)的個(gè)數(shù)n=M；,cin?n；cout<<endl;ocoutvv”請(qǐng)輸入X和Y：”vvendI；〃輸入給定數(shù)據(jù)or(i=0;i<n；i++)?cout?nX["?i?H]二";。cin?X[iJ;?osumX[l]+=X[i];?cout?,,Y[,,?i?M]=M;??cin?Y[i]；?dsumY[l]+=Y[i]；acout?endl;■}?cout?MsumX[1]=n?sumX[1]<<"\t,,?MsumY[l]="?sumY[lJ<<endi;coutvV”請(qǐng)輸入擬合次數(shù)index=M;。cin>>index;?cout<<end1;i=n;sumX[()]=i；for(i=2；iv=2*index；i++)?sumX[i]=0;?for(j=O;j<n；j++)?sumX[i]+=power(X[j],i);^cout<<MsumX[,,?i?nl="?sumX[il<<end1;?for(i=2;i<=index+1;i++)q?sumY[i]=0;?for(j=0;j<n;j++)?sumY[i]+=power(X[j],i—1)*Y[j];?coutvv"sun】Y[,,<<i?,，l=H?sumY[il?endl；。}for(i=1；i<=index+1;i++)。{//建立正規(guī)方程組??for(j=1;j<=index+l;j++)a[i]Ej]=sumX[i+j-2];。b[i]=sumY[i]；)k=l；〃用高斯消元法解方程組，do{。for(j=k+l;j<=index+l;j++)I[j][k]=a[j][k]/a[k][k];?for(i=k+1；i<=index+l;i++)(for(j=k+1；j<=index+1;j++)Mi][j]=a[i][j]-l[i][kl*a[k][j]；gblij=b[i]-l[i][k]*b[k]；aif(k==index+l)break;。k++；}while(l)；x[index+1]=b[index+1]/a[index4-1][index+1];?for(i=index；i>=l;i-){9s=0；ofor(j=i+1;j<=index+l;j++)。s=s+a[i][j]*xU]；，x[i]=(b[i]-s)/a[i][i];?cout?"擬合系數(shù)為：";//輸出擬合系數(shù)afor(i=l;i<=index+1;i++)?cout?x[il?M\t";doub1em=0;coutvVendlvv”平方誤差為:"；for(i=0;i<n；i++)doublet=x[l]+x[21*X[i]-Y[i]；gm=m+power(t,2);Ocout?m?endi;)實(shí)驗(yàn)結(jié)果和分析1代數(shù)插值x(0)-0.40xtl1-0.55y【lJ-0.57815x(2J-0.65VE21-0.6967Sxt31-0.814</[31-0.88811x(4J-0.90J-l.02652x(5J-l.05y153-1.25386詩(shī)懈人插值xx，0.596近似值為，(0-596000>-0?6319232最小二乘法請(qǐng)輸入已知點(diǎn)的個(gè)數(shù)n=7請(qǐng)輸入X和丫：X[0]=0V[0]=1X[U=0.5V[l]=1.75Xt2]=0.6￥E2]=1.96Xt3]=0.7V[31=2.19XE4]=0.8?[4]=2.44XC51-0.9￥[51=2.71X[6]=1.0VE6]=3.00半:sunX[U=4.5sun?[11=15.05請(qǐng)輸入擬合次數(shù)index=2sunX[2]=3.55sunX[3]=2.925sunX[43=2.4979sum?[23=10.975sunY[31=8.9729擬盒系婺力門11平芳誤差為：2.4979Pressanykeytocontinue實(shí)臉?lè)治觯?拉格朗日插值的插值多項(xiàng)式特別容易建立，缺陷是增長(zhǎng)節(jié)點(diǎn)是原有多項(xiàng)式不能運(yùn)用，必須重新建立，即所有基函數(shù)都要重新計(jì)算，這就導(dǎo)致計(jì)算量的浪費(fèi)。2數(shù)據(jù)擬合的具體作法是：對(duì)給定的數(shù)據(jù)(2，%)(i=0,1m),在取定的函數(shù)類中，求p(x)屬于此函數(shù)類,使誤差n=p(Xi)—yi(i=0,1,…，m)的平方和最小，即Sri2=S(Sp(xj)-yi)2=min從幾何意義上講，就是尋求與給定點(diǎn)(X,，//)(i=0,1,???,m)的距離平方和為最小的曲線y=p(x)。心得體會(huì)(碰到的問(wèn)題和解決方法)1通過(guò)做代數(shù)差值實(shí)驗(yàn)，加深了對(duì)牛頓插值和拉格朗口差值的理解，體現(xiàn)了這些算法的優(yōu)點(diǎn)，提高了自己的能力，鞏固了知識(shí)。2本次實(shí)驗(yàn)設(shè)計(jì)，很好地實(shí)現(xiàn)了最小二乘法的程序模擬。感覺(jué)收獲很大。特別是對(duì)平方誤差的計(jì)算模塊，一次成功。最后還對(duì)程序進(jìn)行優(yōu)化,刪去冗余以及標(biāo)準(zhǔn)化模塊。很好，堅(jiān)持下去便是。eIseif(s>0)(b=x；)printf("％f\t%f\n”,a,b);}printf("%AnM,x);printf(M%f\nn,s);return0；}2.割線法：#includeustdio.h#include"math.hMintmain()(floatc,a=l.0,b=2.0;whi1e(l)(c=b-(b*b51tb+4*b*b-10)*(b-a)/(b*b*b+4*b*b—(a*a*a+4*a*a));0if(fabs(b-c)<0.5*0.00001)break;b=c;printf("％f\n0,b);)printf(H%f\nM,c);}流程圖；（\「）「【j32H廠IiI運(yùn)營(yíng)結(jié)果；1二分法

■■卯楠梅，「TTTT1.2500001.250000■■卯楠梅，「TTTT1.2500001.2500001.3125001.3437501.3593751.3593751.3632811.3632811.3642581.3647461.3649901.3651121.3651731.3652041.3652191.3652271.3652271.3652291.3652301.3652300.0000011.5000001.3750001.3750001.3750001.3750001.3671881.3671881.3652341.3652341.3652341.3652341.3652341.3652341.3652341.3652341.3652341.3652311.3652311.3652311.3650381.3652701.3652221.3652321.3652302割線法實(shí)驗(yàn)結(jié)果和分析兩種方法均能求出方程的解，但割線法比二分法的收效速度更快，且程序的代碼更簡(jiǎn)潔。心得體會(huì)（碰到的問(wèn)題和解決方法）通過(guò)實(shí)驗(yàn)，加深了對(duì)方程求根方法的理解，加強(qiáng)了實(shí)踐操作能力，實(shí)現(xiàn)了理論和實(shí)踐相結(jié)合。

實(shí)驗(yàn)名稱實(shí)驗(yàn)二線性方程組的直接求解實(shí)驗(yàn)?zāi)康暮鸵?guī)定合理運(yùn)用①②③④JGauss消元法、LU分解法、追趕法求解下列方程組：-1231x,1F14012x2=8_24lJ[x3J[130.3x10-1559.143iTxJ「59.1715.291-6.130-12x2_46.7811.29529—1121lJ|_x4J[2_-4215Tx,][-287210x2_-74836x3--7126112()J[x4J|_-3_~21Tx,1r-7-21x2-5:(n=5,10,100,??)211-5_12h_|[-5]實(shí)驗(yàn)內(nèi)容高斯消元：1ik二aik/akkai尸aij—1ik*aj(k=l,2,…，n—1i=k+l,k+2,???,nj=k+l,k+2,???,n+1)由回代過(guò)程求得原方程組的解：xn=ann+1/annxk=(akn+i—Sakjxj/akk追趕法：當(dāng)矩陣A為三對(duì)角矩陣，在A的LU分解中，L取下三角陣，U取單位上三角陣，這樣求解方程組Ax=d的方法稱為追趕法.LU分解法：將系數(shù)矩陣A轉(zhuǎn)化為A=L*U,L為單位下三角矩陣,U為普通上三角矩陣,然后通過(guò)解方程組l*y=b,u*x=y,來(lái)求解xo重要儀器設(shè)備臺(tái)式或筆記本計(jì)算機(jī)實(shí)驗(yàn)記錄(寫出實(shí)驗(yàn)內(nèi)容中的程序代碼和運(yùn)營(yíng)結(jié)果)(可分欄或加頁(yè))高斯消元法#inc1ude<stdio.h>#definen3main(){inti,j,k；floata[n][n],c[n][n],b[n],d[n];for(i=0;i<n；i++){for(j=0;j<n;j++){scanf(u%fH,&a[i][j]);c[i][j}scanf(〃％f〃,&b[i]);d[i]=b[i];}for(k=0；k<n;k++)[b[k]=d[k]/c[k][k];for(i=0;i<n；i++){if(i==k)continue;c[i]Lk]=c[i][k]/cEk][k];for(j=k+1;j<n;j++){a[k][j]=c[k][j]/c[k][k]；aLi][j]=c[i][j]-c[i][k]*c[k][j];}b[i]=d[i]-c[i][k]*d[k];}for(i=0;i<n；i++){d[i]=b[i];for(j=k+1;j<n;j++)c[i][j]=a[i][j];}}for(i=0;i<n;i++)printf(z，b[%d]=%f\n〃，i,b[i]);)LU分解法：include<stdio.h>include<math,h>#dcfincL3()doublea[L][L],b[L],1[L][L],u[L][L],x[L],y[L];intmain()(intn,i,j,k,r;oprintf("請(qǐng)輸入矩陣無(wú)次:\n")；scanf("%d",&n);叩rintf("請(qǐng)輸入矩陣各