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

本科实验报告

课程名称： 计算数值方法

实验地点： 综合楼五层506室

专业班级：计科1002 学号： 2010001414

学生姓名： xxx

指导教师： 王峥

2012 年 6 月 20

太原理工大学学生实验报告

学院名称 学生姓名 课程名称 计算机科学与技术 xxx 计算数值方法 专业班级 实验日期 实验题目 计科 1002 2012.6.5 实验一 方程求根 学号 成绩 2010001414 一、课题名称 方程求根：熟悉使用、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。选择上述方法中的两种方法求方程：二分法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）.首先对给定的计算公式进行变形使其能够迭代或者找出相应迭代速度较快的式子。 xk (k0,1,2,)1(xk) 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.0 C-free CodeBlocks 五、结构程序设计 迭代法: #include #include main() { int i; double xn[15],y,x1,x2,m ; printf(\"请输入x1,x2的值：\\n\" ); scanf(\"%lf%lf\ printf(\"请输入精度要求：\\n\" ); scanf(\"%lf\ printf(\" n xn\\n\"); i=0; do{ xn[0]=(x1+x2)/2 ; xn[i+1]= sqrt(10/(4+xn[i])); //迭代 printf(\"%5d %5lf\\n\ y= fabs(xn[i+1]-xn[i]) ; i++; if(y #include main() { int m,n,o,p; double a,b,l; printf(\"请输入x^3, x^2, x的系数和常数p：\\n\"); scanf(\"%d%d%d%d\ 4 0 -10 printf(\"请输入x1,x2:\\n\"); scanf(\"%lf%lf\ printf(\"请输入精度要求:\\n\"); scanf(\"%lf\ printf(\" n an bn xn f(xn)\\n\"); double x,fx; int i=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++; if(fx==0) break; if(fx>0) b=x; else if(fx<0) a=x; if((b-a)指导教师 王峥

太原理工大学学生实验报告

学院名称 学生姓名 课程名称 计算机科学与技术 xxx 计算数值方法 专业班级 实验日期 实验题目 计科 1002 2012.6.5 学号 成绩 2010001414 实验二 线性方程组的直接解法 一、课题名称 线性方程组的直接解法 合理利用Gauss消元法、LU分解法、追赶法求解下列方程组： 0.31015123x114x8 ② 5.291012① 211.2241x133148③41227861x159.176.13012x246.78 952x31x2211459.14321x1715x12121x52210x27 ④（n=5,10,100…） 36x37x121n151120x43125xn二、目的和意义 （1）了解线性方程组常见的直接解法，如Guass消元法、LU分解法、追赶法。 （2）加深对线性方程组求解方法的认识，掌握算法。 （3）会进行误差分析，并能对不同方法进行比较。 三、计算公式 高斯分解法： ⑴将原方程组化为三角形方阵的方程组： lik=aik/akk aij= aij- lik* akj k=1,2,…,n-1 i=k+1,k+2, …,n j=k+1,k+2, …,n+1 ⑵由回代过程求得原方程组的解： xn= ann+1/ ann xk=( akn+1-∑akj xj)/ 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.0 C-free CodeBlocks 五、结构程序设计 Gauss消元法： #include #include using namespace std; int main() { int n,i,j,k; double a[100][100],b[100],o; cout<<\"输入未知数个数:\"<>n; cout<<\"输入数列:\"<>a[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;i>0;i--) { b[i]=a[i][n+1]/a[i][i]; for (j=i-1;j>0;j--) a[j][n+1]=a[j][n+1]-b[i]*a[j][i]; } cout<<\"解得:\"< #include #define N 20 using namespace std; void load(); float a[N][N]; int m; int main(){ int i,j; int c,k,n,p,r; float x[N],l[N][N],s,d; cout<<\"下面请输入未知数的个数m=\"; cin>>m; cout<fabs(a[i][i]))?j:i; /* 找列最大元素 */ for(n=0;n=0;i--) { d=0; for(j=i+1;j>a[i][j]; } LU分解法： #include void solve(float l[][100],float u[][100],float b[],float x[],int n) {int i,j; float t,s1,s2; float y[100]; for(i=1;i<=n;i++) /* 第一次回代过程开始 */ {s1=0; for(j=1;j=1;i--) /* 第二次回代过程开始*/ { s2=0; for(j=n;j>i;j--) { t=-u[i][j]; s2=s2+t*x[j]; } x[i]=(y[i]+s2)/u[i][i]; } } void main() {float a[100][100],l[100][100],u[100][100],x[100],b[100]; int i,j,n,r,k; float s1,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\ printf (\"读取原矩阵 A（x的系数）:\\n\");/*读取原矩阵A*/ for(i=1;i<=n;i++) for(j=1;j<=n;j++) scanf(\"%f\ printf (\"读取列矩阵 B（y的值）:\\n\");/*读取列矩阵B*/ for(i=1;i<=n;i++) scanf(\"%f\ 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\");//输出矩阵L for(i=1;i<=n;i++) { for(j=1;j<=n;j++) printf(\"%7.3f \ printf(\"\\n\"); } printf(\"输出矩阵 U:\\n\");//输出矩阵U for(i=1;i<=n;i++) { for(j=1;j<=n;j++) printf(\"%7.3f \ printf(\"\\n\"); } solve(l,u,b,x,n); printf(\"解为:\\n\"); for(i=1;i<=n;i++) printf(\"x%d=%f\\n\ } 追赶法： #include #define N 3 main() { double A[3][3],b[3]; printf(\"请按顺序输入x的系数:\\n\"); int a,c; {for(a=0;a<3;a++) for(c=0;c<3;c++) scanf(\"%lf\printf(\"请按顺序输入y的值:\\n\"); int k; for(k=0;k<3;k++) scanf(\"%lf\ int i; A[0][1]=A[0][1]/A[0][0]; for(i=1;i<2;i++) A[i][i+1]=A[i][i+1]/(A[i][i]-A[i][i-1]*A[i-1][i]); for(i=1;i<3;i++) A[i][i]=A[i][i]-A[i][i-1]*A[i-1][i]; b[0]=b[0]/A[0][0]; for(i=1;i<3;i++) b[i]=(b[i]-A[i][i-1]*b[i-1])/A[i][i]; for(i=1;i>=0;i--) b[i]=b[i]-A[i][i+1]*b[i+1]; for(i=0;i<3;i++) printf(\"x[%d]=%.6lf\\n\} 六、结果讨论和分析 Gauss消元法： 列主元素消元法： LU分解法： 追赶法： 分析讨论 从消元过程可以看出，对于n阶线性方程组，只要各步主元素不为零，经过n-1步消元，就可以得到一个等价的系数矩阵为上三角形阵的方程组，然后再利用回代过程可求得原方程组的解. 由于列主元素法相似且优于完全主元素法 所以省略了后者。消元过程相当于分解 A为单位下三角阵L与上三角阵U的乘积，解方程组Ly=b回代过程就是解方程组Ux=y。其中的L为n阶单位下三角阵、U为上三角阵. 在 A 的LU 分解中, L取下三角阵, U 取单位上三角阵,这样求解方程组Ax=d 的方法称为追赶法。另外是追赶法和其他方法求同一方程结果不一样，我多次修改源程序，也不知道原因。再就是追赶法有很大的局限性 还待改良。 流程图： 实验地点 综合楼五层506室 指导教师 王峥

太原理工大学学生实验报告

学院名称 学生姓名 课程名称 计算机科学与技术 xxx 计算数值方法 专业班级 实验日期 实验题目 计科 1002 2012.6.5 学号 成绩 2010001414 实验三 线性方程组的迭代解法 一、课题名称 线性方程组的迭代解法 使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。 10x1x22x37.2x110x22x38.3xx5x4.2231 二、目的和意义 学习使用雅可比迭代法或高斯-赛德尔迭代法 三、计算公式 雅克比迭代法： 设线性方程组 Ax=b 的系数矩阵A可逆且主对角元素a11,a22,…,ann均不为零，令 D=diag(a11,a22,…,ann) 并将A分解成 A=(A-D)+D 从而线性方程组可写成 Dx=(D-A)x+b 则有迭代公式 x其中，B1=I-D -1(k+1)=B1x+f1 (k)A,f1=Db。 -1 四、主要仪器设备 Vc++ 9.0 C-free CodeBlocks 五、结构程序设计 雅克比迭代法： #include #include main() { int i; double x1[20] ,x2[20],x3[20]; double x10, x20, x30; printf(\"请输入x1,x2,x3的初值：\\n\"); scanf(\"%lf%lf%lf\printf(\" n x1[n] x2[n] x3[n] \\n\"); for(i=0;i<18;i++) { x1[0]=x10; x2[0]=x20; x3[0]=x30; x1[i+1]=0.1*x2[i]+0.2*x3[i]+0.72; x2[i+1]=0.1*x1[i]+0.2*x3[i]+0.83; x3[i+1]=0.2*x1[i]+0.2*x2[i]+0.84; printf(\"%5d %5lf %5lf %5lf\\n\ } } 六、实验结果与分析： 雅克比迭代法： 分析讨论： 其实，这两个迭代法是之前迭代法的升级，多了几个迭代式子而已，而且两者相差不大比较简单，所以选择了雅克比迭代法进行求解，但是没有与另一种方法 高斯赛德尔迭代法进行实质性的比较。 流程图： 实验地点 综合楼五层506室

指导教师 王峥 太原理工大学学生实验报告

学院名称 学生姓名 课程名称 计算机科学与技术 xxx 计算数值方法 专业班级 实验日期 实验题目 计科 1002 2012.6.10 学号 成绩 2010001414 实验四 矩阵特征值与特征向量问题 一、课题名称 使用幂法求A模为最大的特征值及其相应的特征向量。 210 A121012二、目的和意义 （1）了解矩阵特征值与特征向量问题解法，掌握幂法。 （2）加深对矩阵特征值与特征向量问题求解方法的认识，掌握算法。 三、计算公式 幂法：由已知的非零向量x0和矩阵A的乘幂构造向量序列{xn}以计算矩阵A的按模最大特征值及其特征向量的方法，称为幂法。 迭代公式： yAxkk1mkmax(yk),k1,2,...yxkkmk 结果可取 1mkyk或1xk 1 四、主要仪器设备 Vc++ 9.0 C-free CodeBlocks 五、结构程序设计 源代码： #include #include #define N 3 #define eps 1e-6 #define KM 30 float MaxValue(float x[],int n) { float Max=x[0]; int i; for (i=1;ifabs(Max))Max=x[i]; return Max; } void PowerMethod(float *A) { float U[N],V[N],r1,r2,temp; int i,j,k=0; while(k太原理工大学学生实验报告

学院名称 学生姓名 课程名称 计算机科学与技术 xxx 计算数值方法 专业班级 实验日期 实验题目 计科 1002 2012.6.10 实验五 代数插值 学号 成绩 2010001414 一、课题名称 （使用拉格朗日插值法或牛顿插值法求解：已知f(x)在6个点的函数值如下表所示，运用插值方法，求f(0.596)的近似值。 x f(x) 0.40 0.41075 0.55 0.57815 0.65 0.69675 0.80 0.88811 0.90 1.02652 1.05 1.25386 二、目的和意义 学习使用拉格朗日插值法或牛顿插值法求解 三、计算公式 设函数在区间[a,b]上n+1互异节点x0,x1,…,xn上的函数值分别为y0,y1,…,yn,求n次插值多项式Pn(x),满足条件 Pn(xj)=yj, j=0,1,…,n 令 Ln(x)=y0l0(x)+y1l1(x)+…+ynln(x)= ∑yili(x) 其中l0(x),l1(x),…, ln(x) 为以x0,x1,…,xn为节点的n次插值基函数，则Ln(x)是一次数不超过n的多项式，且满足 Ln(xj)=yj, L=0,1,…,n 再由插值多项式的唯一性，得 Pn(x)≡Ln(x) 四、主要仪器设备 Vc++ 9.0 C-free CodeBlocks 五、结构程序设计 #include #include main() { int I; char L; double M[100][100]; double x[100],y[100]; double X=1,xx=0,w=1,N=0,P,R=1; int n=5; //cout<<\"请输入所求均差阶数：\"; //求所有阶差 //cin>>n; //for(int i=0;i<=n;i++) //{ /*cout<<\"请输入x\"<>x[i]; cout<<\"请输入y\"<>y[i]; M[i][0]=x[i]; M[i][1]=y[i]; *///用二维保存所有数据 M[0][0]=0.40;M[0][1]=0.41075; M[1][0]=0.55;M[1][1]=0.57815; M[2][0]=0.65;M[2][1]=0.69675; M[3][0]=0.80;M[3][1]=0.88811; M[4][0]=0.90;M[4][1]=1.02652; M[5][0]=1.05;M[5][1]=1.25386; for( int j=2;j<=n+1;j++) { for(int i=1;i<=n;i++) { M[i][j]=(M[i][j-1]-M[i-1][j-1])/(M[i][0]-M[i-j+1][0]); } } for(int i=1;i<=n;i++) { cout<<\"其\"<>xx; for(int i=0;i指导教师 王峥

太原理工大学学生实验报告

学院名称 学生姓名 课程名称 计算机科学与技术 xxx 计算数值方法 专业班级 实验日期 实验题目 计科 1002 2012.6.10 学号 成绩 2010001414 实验六 最小二乘法拟合多项式 一、课题名称 给定数据点（xi ,yi），用最小二乘法拟合数据的多项式，并求平方误差。 xi yi 0 1 0.5 1.75 0.6 1.96 0.7 2.19 0.8 2.44 0.9 2.71 1.0 3.00 二、目的和意义 1．熟练运用已学计算方法求解方程组 2．加深对计算方法技巧，选择正确的计算方法来求解各种方程组 3．培养使用电子计算机进行科学计算和解决问题的能力 三、计算公式 建立正规方程组： ∑（∑xi）ak=∑xi 平方误差： j+kjyi ,j=0,1,…,n I=∑（∑akxik-yi）2 四、主要仪器设备 Vc++ 9.0 C-free CodeBlocks 五、结构程序设计 源代码： #include #include #define N 15 double power(double &a,int n) { double b=1; for(int i=0;i>n; n=7; cout<>X[i]; } sumX[1]+=X[i]; //cout<<\"Y[\"<>Y[i]; sumY[1]+=Y[i]; //cout<>index; cout<=1;i--) { } cout<<\"拟合系数为:\"; //输出拟合系数 for(i=1;i<=index+1;i++) cout<综合楼五层506室 指导教师 王峥