2023年计算机数值方法实验报告

TAIYUANUNIVERSITYOFTECHNOLOGY本科实验报告课程名称：计算机数值方法实9佥项目：方程求根线性方程组的直接解法线性方程组的迭代解法

代数插值和最小二乘法拟合多项式实验地点：逸夫302专业班级：软件学号:学生姓名：指导教师:田华2023年4月24日

实验结果和分析12314012824113b[01=1.000000b[11=2.00000011高斯消元法1高斯消元法b[2]=3.1高斯消元法请输入矩阵阶数n：蓍输入矩阵各项：42158721048361261120请输入方程组的常数项: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实脸分析：高斯消元法,是先消元,再回带的过程。由程序段可以发现,始终消去对角线下方的元素。从消元过程可以看出，对于n阶线性方程组，只要各步主元素不为零，通过n-1步消元，就可以得到一个等价的系数矩阵为上三角形阵的方程组，然后再运用回代过程可求得原方程组的解。LU分解法，分解矩阵为单位下三角阵L与上三角阵U的乘积，然后解方程组Ly=b,回代，解方程组Ux=y。其中的L为n阶单位下三角阵、U为上三角阵.对于追赶法，追赶法是合用于三角矩阵的线性方程组的求解的方法，并不合用于其他类型矩阵。心得体会（碰到的问题和解决方法）本证实验难度比较式在编译时经常出现各种错误「程序代码也比较繁琐「谏深感觉到自己的上机操作能力有限，应加强自己的编程能力，以后•要继续努力。

实验名称实验三线性方程组的迭代求解实验目的和规定使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。jlO不一“一2巧=72(一巧+10巧-2巧=83°［一天一巧+54=42实验内容设线性方程组Ax=b的系数矩阵A可逆,且主对角元素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重要仪器设备台式或笔记本计算机实验记录(写出实验内容中的程序代码和运营结果)(可分栏或加页)雅可比迭代法#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);)高斯一赛德尔迭代法#includeiostream”#include"iomanip”usingnamespacestd;intmain(){intiJ,k=(),m,n;doublee1,e2=0.0；coutvv”请输入精度e：H；，cin»el;coutvv”请输入系数矩阵行数：:dcin»m;coutvv”请输入系数矩阵列数：”；in»n；cout<<endl；doubIe(**a)=newdouble*[m]；〃生成二维动态数组or(i=0；i<=m;i++)(ga[i]=newdouble[n];)，double(*b)=newdouble[m]；ddouble(*x)=newdouble[nJ;cout«"请输入系数矩阵:"vVendl;cout«”——--"————一—・…n«endl;for(intnum1=0;numl<m;numl++)。for(intnum2=0;num2<n；num2++)(g，cin>>a[num1][num2];gcoutvvend1;cout<<”输入的系数矩阵为了VVendl；dfor(intnum3=0;num3<m；num3++){gfor(intnum4=0;num4<n;num4++)(®®dcout«a[num3][num4]«H";®cout<<endl;cout-------------------M«endl;ocoutvv”请输入矩阵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;}实验结果和分析1雅克比迭代k»[0]=l.099811&[U=1.199811&[23=1.299777信按任意键继续一.2高斯一一赛德尔迭代请输入系数矩阵:10-1-2-110-2-1-15请输入矩阵b,7.28.34.2值值值值值值值值fi:值值值次KK-弋代弋代弋M-・8啾欹■欲漱欹滋的第电r¥2ly.3l¥4ly5说曾士值值值值fi:值值值次KK-弋代弋代弋值值值值fi:值值值次KK-弋代弋代弋M-・8啾欹■欲漱欹滋的第电r¥2ly.3l¥4ly5说曾士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心得体会(碰到的问题和解决方法)本次实验，让我对这两种方法更加理解，在编程操作上也更加纯熟，此后继续努力，不断丰富自己的知识，增强操作能力。实验名称实验四代数插值和最小二乘法拟合实验目的和规定1使用拉格朗日插值法或牛顿插值法求解:已知f(X)在6个点的函数值如下表所示，运用插值方法，求f(().596)的近似值。X0.400.550.650.800.901.05f(X)0.410.5780.696().888111.02651.25386学生姓名实验成绩实验名称实验一方程求根学生姓名实验目的和规定(必填)熟悉使用、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。选择上述方法中的两种方法求方程:二分法f(x)=x3+4x2-10=0在[1,2]内的一个实根，且规定满足精度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；)函数f(x)在区间(x,y)上连续，先在区间(x,y)拟定a与b,若f(a),f(b)异号,说明在区间(a,b)内存在零点，然后求f[(a+b)/2]。假设F(a)V0,F(b)>0,a<b,①假如f[(a+b)/2]=0,该点即为零点；②假如f[(a+b)/2]<0,则区间((a+b)/2,b)内存在零点，(a+b)/2Na；③假如f[(a+b)/2]>0,则区间(a,(a+b)/2)内存在零点,(a+b)/2Wb；返回①重新循环，不断接近零点。通过每次把f(x)的零点所在区间收缩一半的方法，使区间内的两个端点逐步逼近函数零点，最终求得零点近似值。重要仪器设备台式或笔记本计算机实验记录(写出实验内容中的程序代码和运营结果)(可分栏或加页)075157522给定数据点(H,y),用最小二乘法拟合数据的多项式，并求平方误差。Xi0().50.60.70.80.91.011.751.962.192.442.713.00实验内容1设函数在区间［a,b］上n+1互异节点Xo,X1,•••,Xn上的函数值分别为yo,yi,・・・，y2求n次插值多项式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为节点的n次插值基函数，则Ln(x)是一次数不超过n的多项式,且满足Ln(xj)=Yj,L=0,1,…，n再由插值多项式的唯一性，得Pn(X)=Ln(x)2建立正规方程组：£(EXij'k)ak=LXijyi,j=0,1,-,n平方误差：I=E(Sakxik-yi)2对给定数据点{(Xi,Yi)}(i=0,l,…，m),在取定的函数类①中，求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("请输入数据个数n：”)；从几何意义上讲，就是寻求与给定点{(Xi,Yi)}(i=0,1,…,m)的距离平方和为最小的曲线y=p(x)。函数p(x)称为拟合函数或最小二乘解，求拟合函数p(x)的方法称为曲线拟合的最小二乘法。得到的两个关于a0、al为未知数的两个方程组，解这两个方程组得出：a0=(EYi)/m—al(LXi)/ma1=[mZXiYi-(SXiZYi)]/[mEXi2-(EXi)2)]即最终的拟合多项式各项系数重要仪器设备台式或笔记本计算机实验记录(写出实验内容中的程序代码和运营结果)(可分栏或加页)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("请输入插值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”请输入已知点的个数n=M；,cin»n；cout<<endl;ocoutvv”请输入X和Y：”vvendI；〃输入给定数据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”请输入拟合次数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++)。{//建立正规方程组«®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«"拟合系数为：";//输出拟合系数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;)实验结果和分析1代数插值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诗懈人插值xx，0.596近似值为，(0-596000>-0・6319232最小二乘法请输入已知点的个数n=7请输入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请输入拟合次数index=2sunX[2]=3.55sunX[3]=2.925sunX[43=2.4979sum?[23=10.975sunY[31=8.9729拟盒系婺力门11平芳误差为：2.4979Pressanykeytocontinue实脸分析：1拉格朗日插值的插值多项式特别容易建立，缺陷是增长节点是原有多项式不能运用，必须重新建立，即所有基函数都要重新计算，这就导致计算量的浪费。2数据拟合的具体作法是：对给定的数据(2，%)(i=0,1m),在取定的函数类中，求p(x)属于此函数类,使误差n=p(Xi)—yi(i=0,1,…，m)的平方和最小，即Sri2=S(Sp(xj)-yi)2=min从几何意义上讲，就是寻求与给定点(X,，//)(i=0,1,•••,m)的距离平方和为最小的曲线y=p(x)。心得体会(碰到的问题和解决方法)1通过做代数差值实验，加深了对牛顿插值和拉格朗口差值的理解，体现了这些算法的优点，提高了自己的能力，巩固了知识。2本次实验设计，很好地实现了最小二乘法的程序模拟。感觉收获很大。特别是对平方误差的计算模块，一次成功。最后还对程序进行优化,删去冗余以及标准化模块。很好，坚持下去便是。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运营结果；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割线法实验结果和分析两种方法均能求出方程的解，但割线法比二分法的收效速度更快，且程序的代码更简洁。心得体会（碰到的问题和解决方法）通过实验，加深了对方程求根方法的理解，加强了实践操作能力，实现了理论和实践相结合。

实验名称实验二线性方程组的直接求解实验目的和规定合理运用①②③④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]实验内容高斯消元：1ik二aik/akkai尸aij—1ik*aj(k=l,2,…，n—1i=k+l,k+2,•••,nj=k+l,k+2,•••,n+1)由回代过程求得原方程组的解：xn=ann+1/annxk=(akn+i—Sakjxj/akk追赶法：当矩阵A为三对角矩阵，在A的LU分解中，L取下三角阵，U取单位上三角阵，这样求解方程组Ax=d的方法称为追赶法.LU分解法：将系数矩阵A转化为A=L*U,L为单位下三角矩阵,U为普通上三角矩阵,然后通过解方程组l*y=b,u*x=y,来求解xo重要仪器设备台式或笔记本计算机实验记录(写出实验内容中的程序代码和运营结果)(可分栏或加页)高斯消元法#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("请输入矩阵无次:\n")；scanf("%d",&n);叩rintf("请输入矩阵各