计算数值实验报告太原理工大学.docx

计算数值实验报告太原理工大学.docx

《计算数值实验报告太原理工大学.docx》由会员分享，可在线阅读，更多相关《计算数值实验报告太原理工大学.docx（32页珍藏版）》请在冰点文库上搜索。

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

本科实验报告

课程名称：

计算数值方法

实验地点：

综合楼五层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-freeCodeBlocks

五、结构程序设计

迭代法:

#include

#include

main（）

{

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（y

}while

（1）;

}

二分法：

#include

#include

main（）

{

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-10

printf（"请输入x1,x2:

\n"）;

scanf（"%lf%lf",&a,&b）;//12

printf（"请输入精度要求:

\n"）;

scanf（"%lf",&l）;//0.5x10^5

printf（"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（fx>0）b=x;

elseif（fx<0）a=x;

if（（b-a）

break;//进行计算并返值

}while

（1）;

}

六、结果讨论和分析

二分法：

迭代法：

分析讨论：

使用不同的方法，可以不同程度的求得方程的解，不同的方法速度不同，求得的结果也稍有区别，当然和要求精度也有关系。

刚开始的时候用数组对二分法进行求解，发现循环到第二次就无法实现值的传递，于是换了另外一种方法代替了数组。

流程图：

实验地点

综合楼五层506室

指导教师

王峥

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

学院名称

计算机科学与技术

专业班级

计科1002

学号

2010001414

学生姓名

xxx

实验日期

2012.6.5

成绩

课程名称

计算数值方法

实验题目

实验二线性方程组的直接解法

一、课题名称

线性方程组的直接解法

合理利用Gauss消元法、LU分解法、追赶法求解下列方程组：

①

②

③

④

（n=5,10,100…）

二、目的和意义

（1）了解线性方程组常见的直接解法，如Guass消元法、LU分解法、追赶法。

（2）加深对线性方程组求解方法的认识，掌握算法。

（3）会进行误差分析，并能对不同方法进行比较。

三、计算公式

高斯分解法：

将原方程组化为三角形方阵的方程组：

lik=aik/akk

aij=aij-lik*akjk=1,2,…,n-1

i=k+1,k+2,…,nj=k+1,k+2,…,n+1

由回代过程求得原方程组的解：

xn=ann+1/ann

xk=（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消元法：

#include

#include

usingnamespacestd;

intmain（）

{

intn,i,j,k;

doublea[100][100],b[100],o;

cout<<"输入未知数个数:

"<

cin>>n;

cout<<"输入数列:

"<

for（i=1;i<=n;i++）

for（j=1;j<=n+1;j++）

cin>>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<<"解得:

"<

for（i=1;i<=n;i++）cout<

//system（"pause"）;

return0;

}

列主元素消元法：

#include

#include

#defineN20

usingnamespacestd;

voidload（）;

floata[N][N];

intm;

intmain（）{

inti,j;

intc,k,n,p,r;

floatx[N],l[N][N],s,d;

cout<<"下面请输入未知数的个数m=";

cin>>m;

cout<

cout<<"请按顺序输入增广矩阵a:

"<

load（）;

for（i=0;i

{

for（j=i;j

c=（fabs（a[j][i]）>fabs（a[i][i]））?

j:

i;/*

找列最大元素

*/

for（n=0;n

{s=a[i][n];a[i][n]=a[c][n];a[c][n]=s;}/*

将列最大数防在对角线上

*/

for（p=0;p

cout<

cout<

for（k=i+1;k

{

l[k][i]=a[k][i]/a[i][i];

for（r=i;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;j

d=d+a[i][j]*x[j];

x[i]=（a[i][m]-d）/a[i][i];/*

求解

*/

}

cout<<"该方程组的解为:

"<

for（i=0;i

cout<<"x["<

return0;

}

voidload（）

{

inti,j;

for（i=0;i

for（j=0;j

cin>>a[i][j];

}

LU分解法：

#include

voidsolve（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;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;j>i;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"）;//输出矩阵L

for（i=1;i<=n;i++）

{

for（j=1;j<=n;j++）

printf（"%7.3f",l[i][j]）;

printf（"\n"）;

}

printf（"输出矩阵U:

\n"）;//输出矩阵U

for（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[i]）;

}

追赶法：

#include

#defineN3

main（）

{

doubleA[3][3],b[3];

printf（"请按顺序输入x的系数:

\n"）;

inta,c;

{for（a=0;a<3;a++）

for（c=0;c<3;c++）

scanf（"%lf",&A[a][c]）;}

printf（"请按顺序输入y的值:

\n"）;

intk;

for（k=0;k<3;k++）

scanf（"%lf",&b[k]）;

inti;

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",i,d[i]）;

}

六、结果讨论和分析

Gauss消元法：

列主元素消元法：

LU分解法：

追赶法：

分析讨论

从消元过程可以看出，对于n阶线性方程组，只要各步主元素不为零，经过n-1步消元，就可以得到一个等价的系数矩阵为上三角形阵的方程组，然后再利用回代过程可求得原方程组的解.由于列主元素法相似且优于完全主元素法所以省略了后者。

消元过程相当于分解A为单位下三角阵L与上三角阵U的乘积，解方程组Ly=b回代过程就是解方程组Ux=y。

其中的L为n阶单位下三角阵、U为上三角阵.在A的LU分解中,L取下三角阵,U取单位上三角阵,这样求解方程组Ax=d的方法称为追赶法。

另外是追赶法和其他方法求同一方程结果不一样，我多次修改源程序，也不知道原因。

再就是追赶法有很大的局限性还待改良。

流程图：

实验地点

综合楼五层506室

指导教师

王峥

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

学院名称

计算机科学与技术

专业班级

计科1002

学号

2010001414

学生姓名

xxx

实验日期

2012.6.5

成绩

课程名称

计算数值方法

实验题目

实验三线性方程组的迭代解法

一、课题名称

线性方程组的迭代解法

使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。

二、目的和意义

学习使用雅可比迭代法或高斯-赛德尔迭代法

三、计算公式

雅克比迭代法：

设线性方程组

Ax=b

的系数矩阵A可逆且主对角元素a11,a22,…,ann均不为零，令

D=diag（a11,a22,…,ann）

并将A分解成

A=（A-D）+D

从而线性方程组可写成

Dx=（D-A）x+b

则有迭代公式

x（k+1）=B1x（k）+f1

其中，B1=I-D-1A,f1=D-1b。

四、主要仪器设备

Vc++9.0C-freeCodeBlocks

五、结构程序设计

雅克比迭代法：

#include

#include

main（）

{

inti;

doublex1[20],x2[20],x3[20];

doublex10,x20,x30;

printf（"请输入x1,x2,x3的初值：

\n"）;

scanf（"%lf%lf%lf",&x10,&x20,&x30）;

printf（"nx1[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",i,x1[i],x2[i],x3[i]）;

}}

六、实验结果与分析：

雅克比迭代法：

分析讨论：

其实，这两个迭代法是之前迭代法的升级，多了几个迭代式子而已，而且两者相差不大比较简单，所以选择了雅克比迭代法进行求解，但是没有与另一种方法高斯赛德尔迭代法进行实质性的比较。

流程图：

实验地点

综合楼五层506室

指导教师

王峥

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

学院名称

计算机科学与技术

专业班级

计科1002

学号

2010001414

学生姓名

xxx

实验日期

2012.6.10

成绩

课程名称

计算数值方法

实验题目

实验四矩阵特征值与特征向量问题

一、课题名称

使用幂法求A模为最大的特征值及其相应的特征向量。

二、目的和意义

（1）了解矩阵特征值与特征向量问题解法，掌握幂法。

（2）加深对矩阵特征值与特征向量问题求解方法的认识，掌握算法。

三、计算公式

幂法：

由已知的非零向量x0和矩阵A的乘幂构造向量序列{xn}以计算矩阵A的按模最大特征值及其特征向量的方法，称为幂法。

迭代公式：

结果可取

四、主要仪器设备

Vc++9.0C-freeCodeBlocks

五、结构程序设计

源代码：

#include

#include

#defineN3

#defineeps1e-6

#defineKM30

floatMaxValue（floatx[],intn）

{

floatMax=x[0];

inti;

for（i=1;i

if（fabs（x[i]）>fabs（Max））Max=x[i];

returnMax;

}

voidPowerMethod（float*A）

{

floatU[N],V[N],r1,r2,temp;

inti,j,k=0;

while（k

{

k++;

for（i=0;i

{

temp=0;

for（j=0;j

V[i]=temp;

}

for（i=0;i

if（k==1）r1=MaxValue（V,N）;

elser2=MaxValue（V,N）;

if（fabs（r2-r1）

r1=r2;

}

printf（"r=%f\n",r2）;

for（i=0;i

}

voidmain（）

{

floatA[N][N]={{2,-1,0},{-1,2,-1},{0,-1,2}};floatU[N];//A的值

U[0]=1;

U[1]=1;

U[2]=1;//x0的值

PowerMethod（A[0]）;

}

六、结果讨论和分析

分析讨论

由于该程序将A矩阵和x0的值编写在程序中，所以要想修改成其他矩阵时比较麻烦，所以也有一定的局限性。

幂法是一种求任意矩阵A的按模最大特征值及其对应特征向量的迭代算法。

该方法的最大优点是计算简单，容易在计算机上实现，对稀疏矩阵较为适合，但有时收敛速度很慢。

流程图：

实验地点

综合楼五层506室

指导教师

王峥

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

学院名称

计算机科学与技术

专业班级

计科1002

学号

2010001414

学生姓名

xxx

实验日期

2012.6.10

成绩

课程名称

计算数值方法

实验题目

实验五代数插值

一、课题名称

（使用拉格朗日插值法或牛顿插值法求解：

已知f（x）在6个点的函数值如下表所示，运用插值方法，求f（0.596）的近似值。

x

0.40

0.55

0.65

0.80

0.90

1.05

f（x）

0.41075

0.57815

0.69675

0.88811

1.02652

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）为以