太原理工大学计算机数值方法实验报告.docx
《太原理工大学计算机数值方法实验报告.docx》由会员分享，可在线阅读，更多相关《太原理工大学计算机数值方法实验报告.docx（41页珍藏版）》请在冰点文库上搜索。
太原理工大学计算机数值方法实验报告
本科实验报告
课程名称：
计算机数值方法
实验地点：
专业班级：
学号：
学生姓名：
指导教师：
成绩：
年月日
太原理工大学学生实验报告
学院名称
计算机科学与技术学院
专业班级
学号
学生姓名
实验日期
成绩
课程名称
计算机数值方法
实验题目
实验一方程求根
一、实验目的和要求：
（1）实验目的：
熟悉使用二分法、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。
求方程：
f（x）=x*x*x+4*x*x-10=0在[1,2]内的一个实根，且要求满足精度|x*-xn|<0.001
（2）实验要求：
1.应用C,C++或JAVA编出通用程序，源程序要有详细的注释和说明；
2.比较计算结果，对不同方法进行比较分析；
3.实验完成，提交实验结果并写出报告，分析计算结果是否符合问题的要求，找出计算成功的原因或计算失败的教训。
二、实验内容和原理：
（1）增值寻根法：
基本思想为：
从初始值x0开始，按规定的一个初始步长h来增值。
令x（n+1）=x（n）+h,（n=0,1,2...）,同时计算f（x（n+1））.在增值过程中会遇到三种情况：
1.f（x（n+1））=0,此时x（n+1）即为方程根。
2.f（x（n））和f（x（n+1））同号，说明区间内无根。
3.f（x（n））和f（x（n+1））同号，说明区间内有根，则把步长缩小，直至满足精度要求为止，x（n）或x（n+1）就是满足精度的近似根。
（2）二分法：
基本思想为：
设f（x）在[a,b]内连续，且f（a）*f（b）<0，则方程f（x）=0在（a,b）内有实根x*.然后逐步对分区间[a,b]，通过判断两端点函数值乘积的符号，进一步缩小有根区间从而求出满足精度要求的近似值。
（3）牛顿迭代法：
基本思想为给定一个初始值由牛顿法的迭代公式：
x（n+1）=x（n）-f（x（n））/f’（x（n））（n=0,1,2...）逐步求出x（n），直至（x（n+1）-x（n））的绝对值小于给定精度，则x（n+1）即可作为近似值。
（4）双点割线法：
由给出的两个初始近似值，再根据双点割线法迭代公式：
x（n+1）=x（n）-（f（x（n））/（f（x（n））-f（x（n-1）））*（x（n）-x（n-1））,（n=1,2,3...）逐步求出x（n），直至x（n+1）-x（n）的绝对值满足精度，则x（n+1）即可作为近似值。
（5）单点割线法：
由给出的两个初始近似值，再根据双点割线法迭代公式：
x（n+1）=x（n）-（f（x（n））/（f（x（n））-f（x（0）））*（x（n）-x（0））,（n=1,2,3...）逐步求出x（n），直至x（n+1）-x（n）的绝对值满足精度，则x（n+1）即可作为近似值。
三、主要仪器设备：
笔记本电脑
四、操作方法：
源代码：
（1）增值寻根法：
#include
doublefun（doublex）{//原函数
return（x*x*x+4*x*x-10）;//求解方程f（x）=x*x*x+4*x*x-10=0的根，精度为10-3.
}
intmain（）{
doublea=1.25,h=1,x=a;
printf（"初始近似值为：
%lf\n",a）;
do{
if（fun（x）==0）{printf（"根为：
%f",x）;return0;}/*如果初始值函数值为0，则初始值即为根*/
elseif（fun（x）*fun（x+h）>0）/*如果fun（x）和fun（x+h）同号则使X加h并跳出本次循环执行下一次*/
{x=x+h;continue;}
elseif（fun（x）*fun（x+h）<0）//若异号则说明在X和X+h之间存在函数根，则缩//短步长继续寻根
{h=h/10.0;}
}while（h>0.001）;//当不满足精度要求时继续执行循环体
printf（"根为：
%f\n",x）;//跳出循环说明满足精度要求则x可近似作为方程根
return0;
}
（2）二分法：
#include
#include
#defineesp1e-3//精度
doublef（doublex）//原函数
{
return（x*x*x+4*x*x-10）;
}
doubleroot（double（*fun）（double）,doubleleft,doubleright,doubledeviation）//用二分法求方程根
{//其中形参*fun为指向原函数的指针
doublex,y;
while（fabs（right-left）>deviation）//当不满足精度要求继续执行循环体
{
x=（right+left）/2;//求出中点的横坐标
y=fun（x）;//计算中点的函数值
if（y==0）returnx;//如果中点函数值等于0则中点即为所求根
if（y>0）right=x;//若中点函数值大于0则令其为区间右终点值
elseleft=x;//否则令其为左端点值
}
return（right+left）/2;//以中点值作为最终近似值
}
intmain（）
{doublea,b;
printf（"请输入区间左右端点值：
\n"）;
scanf（"%lf%lf",&a,&b）;
printf（"近似根为%lf\n",root（f,1,2,esp））;
return0;
}
（3）牛顿迭代法：
#include
#include
doublef1（doublex）//原函数
{
returnx*x*x+4*x*x-10;
}
doublef2（doublex）//导函数
{
return3*x*x+8*x;
}
doublenewton（doublex0,doublee）
{
doublex1;
do{
x1=x0;//按牛顿迭代公式Xn+1=Xn-（f1（Xn）/f2（Xn））求解
x0=x1-f1（x1）/f2（x1）;
}while（fabs（x0-x1）>e）;//当不满足精度要求时继续执行循环体
returnx0;//满足精度要求时返回Xn+1的值
}
intmain（）
{
doublex0=1.5;//初始近似值
doublee=pow（10,-3）;//精度
printf（"初始近似值为：
%lf\n",x0）;
printf（"近似根为：
%lf\n",newton（x0,e））;
return0;
}
（4）双点割线法：
#include
#include
#defineesp10e-3//精度
doublef（doublex）//原函数
{
return（x*x*x+4*x*x-10）;
}
doublefun（doublex0,doublex1）//双点割线法求近似根
{
doublex;
inti=0;
while（fabs（x1-x0）>esp）//当不满足精度要求继续执行循环体
{
x=x1-f（x1）*（x1-x0）/（f（x1）-f（x0））;
x0=x1;
x1=x;
i++;
}//达到精度要求后跳出循环
printf（"执行循环次数为i=%d\n",i）;
return（x1）;//返回最终近似根
}
main（）
{
doublem,n,result;
printf（"请输入两个双精度初始近似值mn:
\n"）;//输入两个双精度初始近似值
scanf（"%lf%lf",&m,&n）;
result=fun（m,n）;
printf（"result=%lf\n",result）;//输出最后结果，保留四位小数
}
（5）单点割线法：
#include
#include
#defineesp5*（10e-3）//精度
doublef（doublex）//原函数
{
return（x*x*x-3*x-1）;
}
doublefun（doublex0,doublex1）//单点割线法求近似根
{
doublex;
inti=0;
while（fabs（x1-x0）>esp）//当不满足精度要求继续执行循环体
{
x=x1-f（x1）*（x1-x0）/（f（x1）-f（x0））;
x1=x;
i++;
}//达到精度要求后跳出循环
printf（"执行循环次数为i=%d\n",i）;
return（x1）;//返回最终近似根
}
main（）
{
doublem,n,result;
printf（"pleaseinputmn:
\n"）;//输入两个双精度初始近似值
scanf（"%lf%lf",&m,&n）;
result=fun（m,n）;
printf（"result=%.4f\n",result）;//输出最后结果，保留四位小数
}
五、实验结果：
程序运行结果：
（1）增值寻根法：
初始近似值为：
1.250000
根为：
1.360000
Pressanykeytocontinue
（2）二分法：
请输入区间左右端点值：
12
近似根为1.364746
Pressanykeytocontinue
（3）牛顿迭代法：
初始近似值为：
1.500000
近似根为：
1.365230
Pressanykeytocontinue
（4）双点割线法：
请输入两个双精度初始近似值mn:
12
执行循环次数为i=4
result=1.365212
Pressanykeytocontinue
（5）单点割线法：
pleaseinputmn:
12
执行循环次数为i=15
result=1.365200
Processreturned14（0xE）executiontime:
4.411s
Pressanykeytocontinue.
六、结果分析：
（1）通过对不同方法比较可知，收敛速度：
牛顿迭代法>双点割线法>单点割线法,二分法和增值寻根法的收敛速度就相对较慢，且结果相对前三种方法较不精确。
（2）但二分法和增值寻根法的优点是简单，对f（x）只要求连续，局限性是只能用于求实根，不能求复根及偶数重根。
（3）牛顿法则收敛很快，而且可求复根，缺点是对重根收敛较慢，要求函数一阶导数存在。
（4）割线法省去了牛顿法求函数导数过程，简化了计算步骤，且收敛速度也较快，其中双点割线法比单点割线法收敛速度快。
七、心得体会：
本次实验内容较为简单，最主要的是要学会将数学方法转变为程序上机实现，同时要学会对不同方法进行比较，确定时间复杂度最小的算法。
实验地点
指导教师
太原理工大学学生实验报告
学院名称
计算机科学与技术学院
专业班级
学号
学生姓名
实验日期
成绩
课程名称
计算机数值方法
实验题目
实验二线性方程组的直接求解法
一、实验目的和要求:
（1）实验目的：
合理利用Gauss消元法、LU分解法求解下列方程组：
0.001X1+2.000X2+3.000X3=1.000
-1.000X1+3.712X2+4.623X3=2.000
-2.000X1+1.072X2+5.643X3=3.000
（2）实验要求：
1.应用C,C++或JAVA编出通用程序，源程序要有详细的注释和说明；
2.比较计算结果，对不同方法进行比较分析；
3.实验完成，提交实验结果并写出报告，分析计算结果是否符合问题的要求，找出计算成功的原因或计算失败的教训。
二、实验内容和原理:
（1）Gauss消元法：
基本思想为：
对于n阶线性方程组，只要各步主元素不为0，经过n-1步消元，就可以得到一个等价的的系数矩阵为上三角形矩阵的方程组，然后再利用回代过程即可求得原方程的解。
时间复杂度约为O（n3）。
（2）Gauss列主元素消元法：
基本思想：
在用高斯消元法求解方程组时，用作除法的小主元素可能使舍入误差增加，因此需要考虑依次按列选主元素，然后换行使之变到主元素位置上，再进行消元计算。
（3）Gauss完全主元素消元法：
基本思想：
首先在系数矩阵A中选取绝对值最大的元素作为主元素，然后交换相应行和列，进行高斯消元，其次在A
（1）的2~n行及2~n列选取主元素进行消元，依次进行下去。
（4）LU分解法：
当系数矩阵A满足顺序主子式不为0时，可将A分解为为一个单位下三角矩阵L和一个上三角矩阵U的乘积，且分解唯一，然后方程式变为Ly=b,Ux=y,接着先求y，再求出x。
三、主要仪器设备：
笔记本电脑
四、操作方法：
（1）Gauss消元法:
源代码：
/*矩阵A用于存放线性方程组的增广矩阵，向量X表示线性方程组的解,例题为P43例2*/
#include
#include
intmain（）
{
doublem,p,A[10][10],X[10];
intn,i,j,k,q;
charc;
printf（"请输入方程的阶数（小于等于8）：
\n"）;
scanf（"%d",&n）;
for（i=1;i<=n;i++）//方便起见从A[1][1]开始存入数据
{
printf（"请输入增广矩阵第%d行：
\n",i）;
for（j=1;j<=n+1;j++）
scanf（"%lf",&A[i][j]）;
}//for
for（i=1;i{
m=fabs（A[i][i]）;j=i;//m表示矩阵第i列中的最大值，j用来标记最大值所在的行数
for（k=i+1;k<=n;k++）
{
if（fabs（A[k][i]）>m）
{m=fabs（A[k][i]）;j=k;}//易错
}//for
for（q=i;q<=n+1;q++）//交换第i行和第j行
{
p=A[i][q];
A[i][q]=A[j][q];
A[j][q]=p;
}//for
//*************************消去过程*********
for（k=i+1;k<=n;k++）
{
m=-A[k][i]/A[i][i];
for（j=i;j<=n+1;j++）
A[k][j]=A[k][j]+A[i][j]*m;//易错
}//for
}//for
//*************回代过程*************
X[n]=A[n][n+1]/A[n][n];
for（i=n-1;i>=1;i--）
{
p=0.0;
for（j=i+1;j<=n;j++）
p=p+A[i][j]*X[j];
X[i]=（A[i][n+1]-p）/A[i][i];
}//for
printf（"线性方程组的解为：
\n"）;
for（i=1;i<=n;i++）
printf（"x[%d]=%lf\n",i,X[i]）;
c=getchar（）;
return0;
}
（2）Gauss列主元素消元法：
源代码：
/*矩阵A用于存放线性方程组的增广矩阵，向量X表示线性方程组的解,例题为P43例2*/
#include
#include
intmain（）
{
doublem,p,A[10][10],X[10];
intn,i,j,k,q;
charc;
printf（"请输入方程的阶数（小于等于8）：
\n"）;
scanf（"%d",&n）;
for（i=1;i<=n;i++）//方便起见从A[1][1]开始存入数据
{
printf（"请输入增广矩阵第%d行：
\n",i）;
for（j=1;j<=n+1;j++）
scanf（"%lf",&A[i][j]）;
}//for
for（i=1;i{
m=fabs（A[i][i]）;j=i;//m表示矩阵第i列中的最大值，j用来标记最大值所在的行数
for（k=i+1;k<=n;k++）
{
if（fabs（A[k][i]）>m）
{m=fabs（A[k][i]）;j=k;}//易错
}//for
for（q=i;q<=n+1;q++）//交换第i行和第j行
{
p=A[i][q];
A[i][q]=A[j][q];
A[j][q]=p;
}//for
//*************************消去过程*********
for（k=i+1;k<=n;k++）
{
m=-A[k][i]/A[i][i];
for（j=i;j<=n+1;j++）
A[k][j]=A[k][j]+A[i][j]*m;//易错
}//for
}//for
//*************回代过程*************
X[n]=A[n][n+1]/A[n][n];
for（i=n-1;i>=1;i--）
{
p=0.0;
for（j=i+1;j<=n;j++）
p=p+A[i][j]*X[j];
X[i]=（A[i][n+1]-p）/A[i][i];
}//for
printf（"线性方程组的解为：
\n"）;
for（i=1;i<=n;i++）
printf（"x[%d]=%.3lf\n",i,X[i]）;
c=getchar（）;
return0;
}
（3）Gauss完全主元素消元法：
源代码：
/*矩阵A用于存放线性方程组的增广矩阵，向量X表示线性方程组的解*/
#include
#include
intmain（）
{
doublem,p,A[10][10],X[10];
intn,i,j,k,q,l,c,w;
chara;
printf（"请输入方程的阶数（小于等于8）：
\n"）;
scanf（"%d",&n）;
for（i=1;i<=n;i++）//方便起见从A[1][1]开始存入数据
{
printf（"请输入增广矩阵第%d行：
\n",i）;
for（j=1;j<=n+1;j++）
scanf（"%lf",&A[i][j]）;
}//for
for（i=1;i{
m=fabs（A[i][i]）;j=i;l=i;//m表示矩阵第i列中的最大值，j用来标记最大值所在的行数,l用来标记最大值所在的列数
for（k=i;k<=n;k++）
for（c=i;c<=n;c++）
if（fabs（A[k][c]）>m）
{m=fabs（A[k][c]）;j=k;l=c;}//易错，忘加fabs
for（q=i;q<=n+1;q++）//交换第i行和第j行
{
p=A[i][q];
A[i][q]=A[j][q];
A[j][q]=p;
}//for
for（w=i;w<=n;w++）//交换第i列与第c列
{
p=A[w][i];
A[w][i]=A[w][c];
A[w][c]=p;
}
//*************************消去过程*********
for（k=i+1;k<=n;k++）
{
m=-A[k][i]/A[i][i];
for（j=i;j<=n+1;j++）
A[k][j]=A[k][j]+A[i][j]*m;//易错
}//for
}//for
//*************回代过程*************
X[n]=A[n][n+1]/A[n][n];
for（i=n-1;i>=1;i--）
{
p=0.0;
for（j=i+1;j<=n;j++）
p=p+A[i][j]*X[j];
X[i]=（A[i][n+1]-p）/A[i][i];
}//for
printf（"线性方程组的解为：
\n"）;
for（i=1;i<=n;i++）
printf（"x[%d]=%lf\n",i,X[i]）;
a=getchar（）;
return0;
}
（4）LU分解法：
源代码：
#include
#include
intmain（）
{
doublep,A[12][12],x[12],L[12][12],R[12][12],b[12],y[12];
intn,i,j,k,m;
printf（"请输入方程组的阶数（小于等于10）：
\n"）;
scanf（"%d",&n）;
for（i=1;i<=n;i++）
{
printf（"请输入系数矩阵的第%d行：
\n",i）;
for（j=1;j<=n;j++）
scanf（"%lf",&A[i][j]）;
}//for
printf（"请输入右端向量b:
\n"）;
for（i=1;i<=n;i++）
scanf（"%lf",&b[i]）;
for（j=1;j<=n;j++）
R[1][j]=A[1][j];//上三角矩阵R的第一行为A的第一行
for（i=2;i<=n;i++）
L[i][1]=A[i][1]/R[1][1];//求出L的第一列
for（k=2;k<=n-1;k++）
{
for（j=k;j<=n;j++）
{p=0.0;
for（m=1;m<=k-1;m++）
p=p+L[k][m]*R[m][j];
R[k][j]=A[k][j]-p;
}//for
for（i=k+1;i<=n;i++）
{
p=0.0;
for（m=1;m<=k-1;m++）
p=p+L[i][m]*R[m][k];
L[i][k]=（A[i][k]-p）/R[k][k];
}//for
}//for
p=0.0;
for（j=1;j<=n-1;j++）//求R[n][n]
p=p+L[n][j]*R[j][n];
R[n][n]=A[n][n]-p;
y[1]=b[1];//回代过程求y
for（k=2;k<=n;k++）
{
p=0.0
升级会员 