矩阵与数值分析上机实验报告

2/44大连理工大学矩阵与数值分析上机实验报告任课教师：孟兆良教学班号：06院系：电信计算机应用学生姓名：邓会杰学生学号：21306151目录实验一设，分别编制从小到大和从大到小的顺序程序分别计算 41.程序代码 42.实验结果 43.实验分析 5实验二解线性方程组，分别利用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组 5（一）用Jacobi迭代法解线性方程组 61.程序代码 62.实验结果 93.实验分析 10（二）用Gauss-Seidel迭代法解线性方程组 101.程序代码 102.实验结果 143.实验分析 14（三）高斯列主元消去法求方程组 151.程序代码 152.实验结果 173.实验分析 17（四）QR方法方程组 181.程序代码 182.实验结果 223.实验分析 22实验三非线性方程的迭代解法 23（一）求方程的根 231. 程序代码 232.实验结果 243.实验分析 24（二）利用Newton迭代法求多项式 241.程序代码 242.实验结果 263.实验分析 26实验四数值积分 26（一）基于高斯-勒让德积分 261.程序代码 262.输出结果 313.实验分析 31实验五插值与逼近 32一拉格朗日插值 321. 程序代码 322.输出结果 343.实验分析 36二切比雪夫插值节点 361.程序代码 362.输出结果 393.实验分析 40三数据拟合的最小二乘法 401. 程序代码 402.实验结果 443.实验分析 45实验一设，分别编制从小到大和从大到小的顺序程序分别计算并指出两种方法计算结果的有效位数。1.程序代码从小到大（从大到小）顺序求解#include"stdio.h"#include"math.h"#defineN10000 //S1000000时#defineN1000000voidmain(){ floatj; floatSn=0.0;//doubleSn=0.0 for(j=2;j<=N;j++)//从大到小的顺序时for(j=N;j>=2;j--) { Sn=Sn+1000000.0/(j*j-1); } printf("Sn=%f",Sn);}2.实验结果Sn为double类型时（精度16位）：从小到大输出结果是S10000=749900.004999 S1000000=749999.000001从大到小输出结果是S10000=749900.005000 S1000000=749999.000001Sn为float类型时（精度为7位）：从小到大输出结果是S10000=749849.13500 S1000000=749849.125001从大到小输出结果是S10000=749899.989583 S1000000=749998.9895833.实验分析1）当使用精度较高的数据类型计算时，从大到小与从小到大计算没有太大区别2）当使用精度较低时，计算的先后顺序对结果影响很大，当N=10000时，两个计算结果差为50.854538；当N=1000000时，两个计算结果差为149.874582.3）造成是因为从小到大计算过程中，大的数加小的数时，把小数的舍入结果为0，产生了大数吃小数的现象。舍入误差较大。而从大到小计算过程，先安排小的数参加运算，获得了较高的精度。4）两种运算的精度分别是：以double数据类型从大到小的计算结果作为精确值，分别计算float类型的有效位数从小到大计算：S10000的有效位数为：3位 S1000000有效位数是：3位从大到小计算：S10000的有效数位为：7位 S1000000有效位数是：7位实验二解线性方程组，分别利用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组，其中常向量为维随机生成的列向量，系数矩阵具有如下形式，其中为阶矩阵，为阶单位矩阵，迭代法计算停止的条件为：，给出时的不同迭代步数．（一）用Jacobi迭代法解线性方程组1.程序代码#include"stdafx.h"#include"stdio.h"#include<stdlib.h>#include"math.h"#include"time.h"#include"stdlib.h"//usingnamespacestd;#defineN20voidmain(){ doubleA[N*N][N*N]={0}; doubleB[N*N]={0}; doubleX1[N*N]={0}; doubleX2[N*N]={0}; doubletemp=0.0; doublemaxErr=0.0; doubleerr; boolflag=1; inti; intj; //构造矩阵A for(i=0;i<N*N;i=i+N) { for(j=i;j<i+N;j++) { A[j][j]=4; if(j<i+N-1) { A[j+1][j]=-1; A[j][j+1]=-1; } if(j+N<N*N) { A[j+N][j]=-1; A[j][j+N]=-1; } } } //构造常数项B矩阵 for(i=0;i<N*N;i++) { B[i]=1+(int)(10.0*rand()/(RAND_MAX+1.0));//0~10之间的随机数 printf("%f\n",B[i]); } //构造迭代常数矩阵B for(i=0;i<N*N;i++) { B[i]=B[i]/A[i][i]; } //构造迭代矩阵A for(i=0;i<N*N;i=i+N) { for(j=i;j<i+N;j++) { if(j<i+N-1) { A[j+1][j]=-A[j+1][j]/A[j][j]; A[j][j+1]=-A[j][j+1]/A[j][j]; } if(j+N<N*N) { A[j+N][j]=-A[j+N][j]/A[j][j]; A[j][j+N]=-A[j][j+N]/A[j][j]; } A[j][j]=0; } } //计算 intcount=0;//记录迭代次数 while(true) { count++; //迭代X的值，如果flag为1，则把结果写入X2中；如果flag为0，把结果写入X1中 if(flag==1) { flag=0;//flag=0表示结果写入X2 for(i=0;i<N*N;i++) { temp=0; for(j=0;j<N*N;j++) { temp+=A[i][j]*X1[j]; } X2[i]=temp+B[i]; } } else { flag=1;//表明结果写入X1中 for(i=0;i<N*N;i++) { temp=0; for(j=0;j<N*N;j++) { temp+=A[i][j]*X2[j]; } X1[i]=temp+B[i]; } } //计算误差，判断循环结束的条件 for(i=0;i<N*N;i++) { maxErr=0; err=fabs(X1[i]*X1[i]-X2[i]*X2[i]); if(err>maxErr) maxErr=err; } if(maxErr<0.0000000001) { printf("循环了%d次\n",count); break; } printf("循环了%d次\n",count); } //输出线性方程组的解 if(flag==0) { for(i=0;i<N*N;i++) printf("%f,",X2[i]); } else { for(i=0;i<N*N;i++) printf("%f,",X1[i]); } system("PAUSE");}2.实验结果N=10时，迭代次数为582N=15时，迭代次数为1208N=20时，迭代次数为2033N=30时，迭代次数为4317次3.实验分析用雅可比迭代法求解线性方程组，随着矩阵规模的增大，迭代次数和运行时间都增长迅速。（二）用Gauss-Seidel迭代法解线性方程组1.程序代码#include"stdafx.h"#include"stdio.h"#include<stdlib.h>#include"math.h"#include"time.h"#include"stdlib.h"//usingnamespacestd;#defineN30voidmain(){ doubleA[N*N][N*N]={0}; doubleB[N*N]={0}; doubleX1[N*N]={0}; doubleX2[N*N]={0}; doubletemp1=0.0; doubletemp2=0.0; doublemaxErr=0.0; doubleerr; boolflag=1; inti; intj; intk; //构造矩阵A for(i=0;i<N*N;i=i+N) { for(j=i;j<i+N;j++) { A[j][j]=4; if(j<i+N-1) { A[j+1][j]=-1; A[j][j+1]=-1; } if(j+N<N*N) { A[j+N][j]=-1; A[j][j+N]=-1; } } } //构造常数项B矩阵 for(i=0;i<N*N;i++) { B[i]=1+(int)(10.0*rand()/(RAND_MAX+1.0));//0~10之间的随机数 printf("%f\n",B[i]); } //构造迭代常数矩阵B for(i=0;i<N*N;i++) { B[i]=B[i]/A[i][i]; } //构造迭代矩阵A for(i=0;i<N*N;i=i+N) { for(j=i;j<i+N;j++) { if(j<i+N-1) { A[j+1][j]=-A[j+1][j]/A[j][j]; A[j][j+1]=-A[j][j+1]/A[j][j]; } if(j+N<N*N) { A[j+N][j]=-A[j+N][j]/A[j][j]; A[j][j+N]=-A[j][j+N]/A[j][j]; } A[j][j]=0; } } //计算 intcount=0;//记录迭代次数 while(true) { count++; if(flag==1)//flag==1时用X1中的值进行迭代，把新计算的值放入X2中 { for(i=0;i<N*N;i++) { temp1=0.0; temp2=0.0; for(j=i;j<N*N;j++) { temp1+=A[i][j]*X1[j];//旧值 } for(k=0;k<i;k++) { temp2+=A[i][k]*X2[k];//新值 } X2[i]=temp1+temp2+B[i]; } flag=0; } elseif(flag==0)//flag==2时用X2中的值进行迭代，把新计算的值放入X1中。 { for(i=0;i<N*N;i++) { temp1=0.0; temp2=0.0; for(j=i;j<N*N;j++) { temp1+=A[i][j]*X2[j];//旧值 } for(k=0;k<i;k++) { temp2+=A[i][k]*X1[k];//新值 } X1[i]=temp1+temp2+B[i]; } flag=1; } //计算误差，判断循环结束的条件 for(i=0;i<N*N;i++) { maxErr=0; err=fabs(X1[i]*X1[i]-X2[i]*X2[i]); if(err>maxErr) maxErr=err; } if(maxErr<0.0000000001) { printf("循环了%d次\n",count); break; } printf("循环了%d次\n",count); } //输出线性方程组的解 if(flag==0) { for(i=0;i<N*N;i++) printf("%f,",X2[i]); } else { for(i=0;i<N*N;i++) printf("%f,",X1[i]); } system("PAUSE");}2.实验结果N=10时，迭代次数为295次N=15时，迭代次数为661次N=20时，迭代次数为1040次N=30时，迭代次数为2213次3.实验分析1）用高斯消去解线性方程组，随着矩阵规模的增大，迭代次数也迅速增加。2）高斯方法与雅可比相比，在相同矩阵，雅可比迭代次数更多，雅可比迭代次数大约是高斯迭代次数的一倍。二用Gauss列主元消去法、QR方法求解如下方程组：（三）高斯列主元消去法求方程组1.程序代码#include"stdio.h"#include<stdlib.h>#include"math.h"#include"time.h"#include"stdlib.h"voidmain(){ floatA[4][5]={1,2,-1,1,0,2,5,0,3,4,1,7,9,2,12,8,-1,-2,1,-8}; floatX[4]={0}; floatmax=0; intmaxindex=0;//记录最大值元素的行下标。 inti; intj; intk; floattemp; for(i=0;i<3;i++) { //扫描第i列最大的元素的行号， max=0; maxindex=i; for(j=i;j<4;j++) { if(fabs(A[j][i])>fabs(max)) { max=A[j][i]; maxindex=j; } } //找到第i列最大元素的行号，与第i行交换位置； if(maxindex!=i) { for(j=0;j<5;j++) { temp=A[i][j]; A[i][j]=A[maxindex][j]; A[maxindex][j]=temp; } } //高斯消去 for(j=i+1;j<4;j++) { temp=A[j][i];//暂存 for(k=i;k<5;k++) { A[j][k]=-A[i][k]/A[i][i]*temp+A[j][k]; } } printf("第%d迭代结果：\n",i+1); for(k=0;k<4;k++) { for(j=0;j<5;j++) { printf("%f,",A[k][j]); } printf("\n"); } printf("\n"); } //把结果写入X矩阵中 for(i=3;i>=0;i--) { floatsum=0.0; for(j=i+1;j<4;j++) { sum+=A[i][j]*X[j]; } X[i]=(A[i][4]-sum)/A[i][i]; } //输出X[] for(i=0;i<4;i++) { printf("X[%d]=%f\n",i,X[i]); } system("PAUSE");}2.实验结果3.实验分析高斯列主元消去法，避免了小主元做除数，在Gauss消去法中增加选主元的过程，在第K步消元时，首先在第K列主对角元以下元素中挑选绝对值最大的数，并通过初等行变换，使得该数位于主对角上，让后继续消元。（四）QR方法方程组1.程序代码#include"stdio.h"#include"math.h"#defineN4voidmain(){ inti; intj; intk; intl; doubleX[N]; doublesum=0.0; doublesumW=0.0; doubleA[N][N]={1,2,-1,1,2,5,0,3,1,7,9,2,8,-1,-2,1}; doubleB[N]={0,4,12,-8}; doubletempB[N]={0}; doubleW[N]={0}; doubletemp1Q[N][N]={0}; doubletemp2Q[N][N]={0}; for(i=0;i<N-1;i++)//i代表列 { //计算W列向量 sum=0; for(k=i;k<N;k++) { sum+=A[k][i]*A[k][i];//第i列中从第k个元素到第N-1个元素的内积， } for(j=i;j<N;j++) { if(i==j) W[j]=A[j][i]-sqrt(sum); else { W[j]=A[j][i]; } printf("W[%d]=%f\n",j,W[j]); } //计算W的内积 sumW=0; for(j=i;j<N;j++) { sumW+=W[j]*W[j]; } //初始化矩阵 for(j=0;j<N;j++) { for(k=0;k<N;k++) { temp1Q[j][k]=0; temp2Q[j][k]=0; } tempB[j]=0; } //计算temp1Q矩阵 for(j=i;j<N;j++) { for(k=i;k<N;k++) { if(j==k) { temp1Q[j][k]=1-2.0/sumW*W[j]*W[k]; } else { temp1Q[j][k]=-2.0/sumW*W[j]*W[k]; } } } //把temp1Q[][]*A[][]暂时存放在tempQ2中;temp1Q[][]*B[]暂存到tempB[]中 for(j=i;j<N;j++) { for(k=i;k<N;k++) { for(l=i;l<N;l++) { temp2Q[j][k]+=temp1Q[j][l]*A[l][k]; } tempB[j]+=temp1Q[j][k]*B[k]; } } //把temp2Q放入A中tempB放入B中 for(j=i;j<N;j++) { for(k=i;k<N;k++) { A[j][k]=temp2Q[j][k]; } B[j]=tempB[j]; } //输出A printf("A:\n"); for(j=0;j<N;j++) { for(k=0;k<N;k++) { printf("%f,",A[j][k]); } printf("%f\n",B[j]); } } //解出X[] for(i=N-1;i>=0;i--) { floatsum=0.0; for(j=i+1;j<=N-1;j++) { sum+=X[j]*A[i][j]; } X[i]=(B[i]-sum)/A[i][i]; printf("X[%d]=%f\n",i,X[i]); } for(i=0;i<N;i++) printf("X[%d]=%f\n",i,X[i]);}2.实验结果X[0]=-1.000000X[1]=0.000000X[2]=1.000000X[3]=2.0000003.实验分析利用Householder变换对A进行QR分解，Q3Q2Q1A=R；同样Q3Q2Q1B=B’（其中Q3Q2Q1为Q的转置矩阵）Ax=BRx=B’，R为上三角矩阵，可直接解得X。实验三非线性方程的迭代解法（一）求方程的根，迭代停止的条件为：；程序代码#include"stdafx.h"#include"stdio.h"#include<stdlib.h>#include"math.h"#include"time.h"#include"stdlib.h"voidmain(){ doublex=1.0; doubletemp=1.0; intcount=0; while(true) { count++; temp=sqrt((16+log(x)-2*sin(x)-exp(x))/2); printf("%f\n",temp); if(fabs(temp-x)<0.0000000001) { x=temp; printf("x=%f\n",x); floatres=exp(x)+2*x*x+2*sin(x)-

log(x)-16; printf("迭代次数=%d\nres=%

f",count,res); break; } else { x=temp; } } system("PAUSE");}2.实验结果x=1.962923迭代次数为803.实验分析让x的平方反解整个式子，满足收敛条件，然后进行迭代。（二）利用Newton迭代法求多项式的所有实零点，注意重根的问题。1.程序代码#include"stdio.h"#include<stdlib.h>#include"math.h"#include"time.h"#include"stdlib.h"voidmain(){ inti; floattemp=0.0; intf[40]; //判断方程的根在什么位置 for(i=-20;i<20;i++) { f[i]=i*i*i*i-3*i*i*i-3*i*i+11*i-6; printf("f[%d]=%d\n",i,f[i]); } //令x1=-1.5;x2=0.5;x3=3.5 floatx[3]={-2.5,0.5,4}; intcount=0; for(i=0;i<3;i++) while(true) { count++; temp=x[i]-(x[i]*x[i]*x[i]*x[i]-3*x[i]*x[i]*x[i]-3*x[i]*x[i]+11*x[i]-6)/(4*x[i]*x[i]*x[i]-9*x[i]*x[i]-6*x[i]+11); if(fabs(temp-x[i])<0.000001) { x[i]=temp; printf("count=%d,x1=%f\n",count,x[i]); break; } else { x[i]=temp; } } system("PAUSE");}2.实验结果count=5,x1=-2.000000count=24,x1=0.999999count=30,x1=3.0000003.实验分析首先，通过画图大致预测方程的根所属范围，然后分别从各范围中任意找一个点，用牛顿法进行迭代，最后迭代收敛到方程的根。实验四数值积分分别用三点Gauss型求积公式计算积分（一）基于高斯-勒让德积分1.程序代码#include"stdio.h"#include"math.h"#definea0#defineb1doublexjifen(inti)//在[-1,1]区间上，x的i次幂的积分{ if(i==0) { return2; } elseif(i==1) { return0; } elseif(i==2) { return2.0/3; } elseif(i==3) { return0; } else { return-1; }}voidmain(){ doubleArr[3][4]; doubleX[3]={sqrt(3.0/5),0,-sqrt(3.0/5)}; doubleA[3]={0}; inti; intj; intk; doublemax; intmaxindex; doubletemp; //初始化矩阵 for(i=0;i<3;i++) { for(j=0;j<4;j++) { if(j!=3) { Arr[i][j]=pow(X[j],i); } else { Arr[i][j]=xjifen(i); } } } printf("输出Arr矩阵"); for(i=0;i<3;i++) { for(j=0;j<4;j++) { printf("%f,",Arr[i][j]); } printf("\n"); } //用列主元高斯消去法解线性方程组 for(i=0;i<3;i++) { //扫描第i列最大的元素的行号， max=0; maxindex=i; for(j=i;j<3;j++) { if(fabs(Arr[j][i])>fabs(max)) { max=Arr[j][i]; maxindex=j; } } //找到第i列最大元素的行号，与第i行交换位置； if(maxindex!=i) { for(j=0;j<4;j++) { temp=Arr[i][j]; Arr[i][j]=Arr[maxindex][j]; Arr[maxindex][j]=temp; } } //高斯消去 for(j=i+1;j<3;j++) { temp=Arr[j][i];//暂存 for(k=i;k<4;k++) { Arr[j][k]=-Arr[i][k]/Arr[i][i]*temp+Arr[j][k]; } } printf("第%d迭代结果：\n",i+1); for(k=0;k<3;k++) { for(j=0;j<4;j++) { printf("%f,",Arr[k][j]); } printf("\n"); } printf("\n"); } //把结果写入A矩阵中 for(i=2;i>=0;i--) { floatsum=0.0; for(j=i+1;j<3;j++) { sum+=Arr[i][j]*A[j]; } A[i]=(Arr[i][3]-sum)/Arr[i][i]; } //输出A[] for(i=0;i<3;i++) { printf("A[%d]=%f\n",i,A[i]); } //计算积分值 doublesum=0; for(i=0;i<3;i++) { temp=(double)(b-a)*X[i]/2+(double)(b+a)/2; sum+=(b-a)*exp(-temp)*A[i]/2; } printf("sum=%f\n",sum);}2.输出结果3.实验分析1.用勒让德积分法积分区间为[-1,1]，需要积分区间转换，X范围[a,b]，转换到t属于[-1,1]过程是x=(b-a)*t/2+(b+a)*/22.用勒让德积分公式，有三个节点时，正交多项式为：5.0/2*x*x*x+3.0/2*x=0，因此高斯点为x0=sqrt(3.0/5)x1=0x3=-sqrt(3.0/5)1 1 11 1 1X0 x1 x2X0*x0 x1*x1 x2*x2U0U1U2A0A1A2Ui为x的i次方在[-1,1]上的积分解出A0A1A24.求积分：积分为：A0*f(x)+A1*f(x)+A2*f(x)实验五插值与逼近一拉格朗日插值程序代码#include"stdafx.h"#include"stdio.h"#include"math.h"#include"stdlib.h"#defineN8intn=8;inta=-1;intb=1;floath=(float)(b-a)/(n-1);floaty[N];floatx[N];//构造Lk(X)函数floatLfunc(floatx,intk){ inti; floatsum=1; for(i=0;i<n;i++) { if(i!=k) { sum=sum*(x-(a+h*i))/((a+h*k)-

(a+h*i)); } } returnsum;}//构造差值节点voidYfunc(){ inti; //计算差值节点 for(i=0;i<n;i++) { y[i]=1.0/(1+(a+h*i)*(a+h*i)); x[i]=a+h*i; }}//差值函数floatLagrange(floatx){ inti; floatsum=0.0; for(i=0;i<n;i++) { sum+=y[i]*Lfunc(x,i); } returnsum;}voidmain(){ Yfunc(); printf("拉格朗日：%d个插值节点\n",n); for(inti=0;i<n;i++) { printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y

[i]); } floatx; printf("请输入-1到1之间的数\n"); while(1) { printf("x="); scanf("%f",&x); if(x>1||x<-1) { printf("输入错误"); break; } else printf("result：%f\n",Lagrange

(x)); } system("PAUSE");}2.输出结果3.实验分析拉格朗日插值，选取等距节点，随着节点数目的增加，曲线越平滑，越接近真实值二切比雪夫插值节点 1.程序代码#include"stdafx.h"#include"stdio.h"#include"math.h"#include"stdlib.h"#defineN10#definePi3.1415926floatX[N];floatY[N];voidinitX(){ inti; floattemp; for(i=0;i<N;i++) { temp=(2*i+1)*Pi; temp=temp/2; temp=temp/N; X[i]=cos(temp); }}voidinitY(){ inti; for(i=0;i<N;i++) { Y[i]=1.0/(1+X[i]*X[i]); }}floatLfunc(floatx,intk){ inti; floatsum=1.0; for(i=0;i<N;i++) { if(i!=k) { sum=sum*(x-X[i])/(X[k]-X[i]); } } returnsum;}floatfunInsert(floatx){ inti; floatsum=0; for(i=0;i<N;i++) { sum=sum+Y[i]*Lfunc(x,i); } returnsum;}voidmain(){ initX(); initY(); printf("切比雪夫插值节点\n"); for(inti=0;i<N;i++) { printf("X[%d]=%f",i,X[i]); printf("Y[%d]=%f\n",i,Y[i]); } floatx; printf("请输入-1~1之间的值\n"); while(1) { printf("x="); scanf("%f",&x); { if(x<=1&&x>=-1) { floatresult=funInsert(x); printf("result=%f\n",result); } else { printf("请正确输入！\n"); break; } } } system("PAUSE");}2.输出结果3.实验分析 用切比雪夫插值节点插值，与等距节点插值相比，相同数目的节点，切比雪夫节点插值更平滑。三数据拟合的最小二乘法观察物体的直线运动，得到如下数据时刻t00.91.93.03.95.0位移s010305180111求运动学方程。程序代码#include"stdio.h"#include"math.h"#defineN2voidmain(){ in