计算数值实验报告(太原理工大学)

本科实验报告课程名称：计算数值方法实验地点：综合楼五层506室专业班级：计科1002学号：2010001414学生姓名：xxx指导教师：王峥2012年6月20太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号2010001414学生姓名xxx实验日期2012.6.5成绩课程名称计算数值方法实验题目实验一方程求根一、课题名称方程求根：熟悉使用、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。选择上述方法中的两种方法求方程：二分法f(x)=x3+4x2-10=0在[1,2]内的一个实根，且要求满足精度|x*-xn|0.5×10-5迭代法：用迭代公式x=f(x)进行迭代计算，直到满足|x*-xn|0.5×10-5为止。二分法：设f（x）在[a,b]上连续，且f（a1）*f（x1）0,记（a2，b2）=（x1，b1）带入计算式进行计算直到|x*-xn|0.5×10-5为止。二、目的和意义（1）了解非线性方程求根的常见方法，如二分法、迭代法、牛顿法、割线法。（2）加深对方程求根方法的认识，掌握算法。会进行误差分析，并能对不同方法进行比较。三、计算公式（1）迭代法1）.首先对给定的计算公式进行变形使其能够迭代或者找出相应迭代速度较快的式子。2）.带入求好的式子到循环中去比如：（2）二分法：f(x)在区间（x，y）上连续1）.先找到a、b属于区间（x，y），使f(a)，f(b)异号，说明在区间(a,b)内一定有零点，然后求f[(a+b)/2],2）.如果f[(a+b)/2]=0，该点就是零点，如果f[(a+b)/2]0,则在区间（(a+b)/2，b)内有零点，反之在（a，(a+b)/2）内有零点带入1）中继续。四、主要仪器设备Vc++9.0C-freeCodeBlocks1()(0,1,2,)kkxxk五、结构程序设计迭代法:#includestdio.h#includemath.hmain(){inti;doublexn[15],y,x1,x2,m;printf(请输入x1,x2的值：\n);scanf(%lf%lf,&x1,&x2);printf(请输入精度要求：\n);scanf(%lf,&m);printf(nxn\n);i=0;do{xn[0]=(x1+x2)/2;xn[i+1]=sqrt(10/(4+xn[i]));//迭代printf(%5d%5lf\n,i,xn[i]);y=fabs(xn[i+1]-xn[i]);i++;if(ym)break;}while(1);}二分法：#includestdio.h#includemath.hmain(){intm,n,o,p;doublea,b,l;printf(请输入x^3,x^2,x的系数和常数p：\n);scanf(%d%d%d%d,&m,&n,&o,&p);//140-10printf(请输入x1,x2:\n);scanf(%lf%lf,&a,&b);//12printf(请输入精度要求:\n);scanf(%lf,&l);//0.5x10^5printf(nanbnxnf(xn)\n);doublex,fx;inti=1;do{x=(b+a)/2;fx=m*x*x*x+n*x*x+o*x+p;printf(%5d%5f%5f%5f%5f\n,i,a,b,x,fx);i++;if(fx==0)break;if(fx0)b=x;elseif(fx0)a=x;if((b-a)l)break;//进行计算并返值}while(1);}六、结果讨论和分析二分法：迭代法：分析讨论：使用不同的方法，可以不同程度的求得方程的解，不同的方法速度不同，求得的结果也稍有区别，当然和要求精度也有关系。刚开始的时候用数组对二分法进行求解，发现循环到第二次就无法实现值的传递，于是换了另外一种方法代替了数组。流程图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号2010001414学生姓名xxx实验日期2012.6.5成绩课程名称计算数值方法实验题目实验二线性方程组的直接解法一、课题名称线性方程组的直接解法合理利用Gauss消元法、LU分解法、追赶法求解下列方程组：①13814142210321321xxx②2178.4617.5911212592.1121130.6291.51314.59103.0432115xxxx③3772201161263841027851244321xxxx④55572112112112121nnxxxx（n=5,10,100…）二、目的和意义（1）了解线性方程组常见的直接解法，如Guass消元法、LU分解法、追赶法。（2）加深对线性方程组求解方法的认识，掌握算法。（3）会进行误差分析，并能对不同方法进行比较。三、计算公式高斯分解法：⑴将原方程组化为三角形方阵的方程组：lik=aik/akkaij=aij-lik*akjk=1,2,…,n-1i=k+1,k+2,…,nj=k+1,k+2,…,n+1⑵由回代过程求得原方程组的解：xn=ann+1/annxk=(akn+1-∑akjxj)/akk(k=n-1,n-2,…,2,1)LU分解法：将系数矩阵A转化为A=L*U,L为单位下三角矩阵，U为普通上三角矩阵，然后通过解方程组l*y=b,u*x=y,来求解x.追赶法:用来求对角方程组；将系数矩阵A转化为A=L*U,L为普通下n-1对角矩阵，U为单位上n-1对角矩阵，然后通过解方程组l*y=b,u*x=y,来求解x.四、主要仪器设备Vc++9.0C-freeCodeBlocks五、结构程序设计Gauss消元法：#includeiostream#includecmathusingnamespacestd;intmain(){intn,i,j,k;doublea[100][100],b[100],o;cout输入未知数个数:endl;cinn;cout输入数列:endl;for(i=1;i=n;i++)for(j=1;j=n+1;j++)cina[i][j];for(i=1;i=n;i++)for(j=i+1;j=n;j++)if(fabs(a[j][i])1e-7){o=a[i][i]/a[j][i];for(k=i;k=n+1;k++)a[j][k]=a[j][k]*o-a[i][k];}for(i=n;i0;i--){b[i]=a[i][n+1]/a[i][i];for(j=i-1;j0;j--)a[j][n+1]=a[j][n+1]-b[i]*a[j][i];}cout解得:endl;for(i=1;i=n;i++)coutb[i]endl;//system(pause);return0;}列主元素消元法：#includeiostream#includecmath#defineN20usingnamespacestd;voidload();floata[N][N];intm;intmain(){inti,j;intc,k,n,p,r;floatx[N],l[N][N],s,d;cout下面请输入未知数的个数m=;cinm;coutendl;cout请按顺序输入增广矩阵a:endl;load();for(i=0;im;i++){for(j=i;jm;j++)c=(fabs(a[j][i])fabs(a[i][i]))?j:i;/*找列最大元素*/for(n=0;nm+1;n++){s=a[i][n];a[i][n]=a[c][n];a[c][n]=s;}/*将列最大数防在对角线上*/for(p=0;pm+1;p++)couta[i][p]\t;coutendl;for(k=i+1;km;k++){l[k][i]=a[k][i]/a[i][i];for(r=i;rm+1;r++)/*化成三角阵*/a[k][r]=a[k][r]-l[k][i]*a[i][r];}}x[m-1]=a[m-1][m]/a[m-1][m-1];for(i=m-2;i=0;i--){d=0;for(j=i+1;jm;j++)d=d+a[i][j]*x[j];x[i]=(a[i][m]-d)/a[i][i];/*求解*/}cout该方程组的解为:endl;for(i=0;im;i++)coutx[i]=x[i]\t;return0;}voidload(){inti,j;for(i=0;im;i++)for(j=0;jm+1;j++)cina[i][j];}LU分解法：#includestdio.hvoidsolve(floatl[][100],floatu[][100],floatb[],floatx[],intn){inti,j;floatt,s1,s2;floaty[100];for(i=1;i=n;i++)/*第一次回代过程开始*/{s1=0;for(j=1;ji;j++){t=-l[i][j];s1=s1+t*y[j];}y[i]=(b[i]+s1)/l[i][i];}for(i=n;i=1;i--)/*第二次回代过程开始*/{s2=0;for(j=n;ji;j--){t=-u[i][j];s2=s2+t*x[j];}x[i]=(y[i]+s2)/u[i][i];}}voidmain(){floata[100][100],l[100][100],u[100][100],x[100],b[100];inti,j,n,r,k;floats1,s2;for(i=1;i=99;i++)/*将所有的数组置零，同时将L矩阵的对角值设为1*/for(j=1;j=99;j++){l[i][j]=0,u[i][j]=0;if(j==i)l[i][j]=1;}printf(输入方程组的个数n:\n);/*输入方程组的个数*/scanf(%d,&n);printf(读取原矩阵A（x的系数）:\n);/*读取原矩阵A*/for(i=1;i=n;i++)for(j=1;j=n;j++)scanf(%f,&a[i][j]);printf(读取列矩阵B（y的值）:\n);/*读取列矩阵B*/for(i=1;i=n;i++)scanf(%f,&b[i]);for(r=1;r=n;r++)/*求解矩阵L和U*/{for(i=r;i=n;i++){s1=0;for(k=1;k=r-1;k++)s1=s1+l[r][k]*u[k][i];u[r][i]=a[r][i]-s1;}for(i=r+1;i=n;i++){s2=0;for(k=1;k=r-1;k++)s2=s2+l[i][k]*u[k][r];l[i][r]=(a[i][r]-s2)/u[r][r];}}printf(输出矩阵L:\n);//输出矩阵Lfor(i=1;i=n;i++){for(j=1;j=n;j++)printf(%7.3f,l[i][j]);printf(\n);}printf(输出矩阵U:\n);//输出矩阵Ufor(i=1;i=n;i++){for(j=1;j=n;j++)printf(%7.3f,u[i][j]);printf(\n);}solve(l,u,b,x,n);printf(解为:\n);for(i=1;i=n;i++)printf(x%d=%f\n,i,x[