数值分析 课程设计.docx
《数值分析 课程设计.docx》由会员分享，可在线阅读，更多相关《数值分析 课程设计.docx（17页珍藏版）》请在冰豆网上搜索。
数值分析课程设计
数值分析课程设计
求解线性方程组
作者姓名：
孙怡丰
学号：
200930980221
指导教师：
张昕副教授
学院名称
理学院
专业名称
统计学
提交日期
2011年6月10日
一、问题的提出
分别用SOR方法和高斯消元的LU分解算法（lii=1,i=1,…,n）求解给定的线性方程组AX=B,以感受迭代法和直接法的不同特点。
二、实验内容
1.自定义函数SOR（A,B,w,MAXN,TOL），以实现SOR方法求解线性方程组AX=B，其中
A——系数矩阵；
B——常数列向量；
w——松弛因子；
MAXN——迭代的最大次数
TOL——
达到的精度上限
返回值有以下四种可能：
a）-2：
SOR方法不收敛；（不收敛的依据为
的某个分量值超出区间[-108,108]。
）
b）-1：
矩阵有一列全为0；
c）0：
算法经过MAXN次迭代还未收敛；
d）k：
SOR方法经k次迭代收敛，求得方程组的解向量X记录下来.
2.自定义函数Direct（A,B），以实现高斯LU分解的方法求解线性方程组AX=B，其中
A——系数矩阵；
B——常数列向量；
返回值有两种可能：
a）“LUdecompsitionfailed.”：
分解过程中U的对角线元素至少一个为0；
b）X：
分解过程中
3.分别使用步骤1中定义的函数SOR（A,B,w,MAXN,TOL）和步骤2中定义的函数Direct（A,B）进行测试，记录返回值及X值（算法收敛或有效的情形,保留4位小数）:
（1）测试1：
MAXN=1000，TOL=10-9，w分别取1，1.05，1.1，1.2，1.3，1.6，1.95；
（2）测试2：
MAXN=1000，TOL=10-9，w=1；
（3）测试3：
MAXN=1000，TOL=10-9，w=1.2；
（4）测试4：
MAXN=1000，TOL=10-9，w=1，1.1，1.3，1.8；
（5）测试5：
:
n阶Hilbert矩阵定义为
取n=3,
MAXN=1000，TOL=10-9，w=1,1.3,1.6,1.9;
测试6：
A为4阶Hilbert矩阵，
MAXN=10000，TOL=10-6，w=1,1.3,1.6,1.8,1.9.
三、实验结果及分析
SOR方法
测试1：
W=1
W=1.05
W=1.1
W=1.2
W=1.3
W=1.6
W=1.95
测试2：
测试3：
测试4：
w=1
w=1.1
w=1.3
w=1.8
测试5：
w=1
w=1.3
w=1.6
w=1.9
测试6：
w=1
w=1.3
w=1.6
w=1.8
w=1.9
高斯LU方法
测试1：
测试2：
测试3：
测试4：
测试5：
测试6：
四、关于本设计的体会
虽然在多数实验中，两种方法的结果是一样的，但是在通过松弛度w的调整中，SOR可能会得出结果，但是在精确度上一般会差点。
高斯在数字的精确度上好点但是，其使用起来并不能全部算出结果，适用的范围相对较小。
所以，两种方法各有优点，在做题目时应结合两种方法，对结果比较后才可以做出正确的判断。
五、参考文献
《数值分析》（第三版）北京理工大学出版社
六、附录
原程序：
SOR:
#include"stdlib.h"
#include"stdio.h"
#include"conio.h"
#include"string.h"
#include"math.h"
#defineN100//初始的一个大矩阵//
floatcre_sch（intn,float*w,floata[N][N],floatb[N],floatlw[N][N],floatfw[N]）
{
inti,j,k;
floattmp1[N][N],tmp2[N][N],rev[N][N],tmp;
for（i=0;ifor（j=0;j{
if（j==i）
{
tmp1[i][j]=a[i][j];
tmp2[i][j]=（1-*w）*a[i][j];
rev[i][j]=1;
}
elseif（j{
tmp1[i][j]=*w*a[i][j];
tmp2[i][j]=0;
rev[i][j]=0;
}
else
{
tmp1[i][j]=0;
tmp2[i][j]=-*w*a[i][j];
rev[i][j]=0;
}
}
for（j=0;jfor（i=0;ifor（k=0;k<=j;k++）
{
if（iif（i==j）
{
rev[i][k]*=1/tmp1[j][j];
continue;
}
rev[i][k]+=rev[j][k]*（-tmp1[i][j]）;
}
for（i=0;ifor（j=0;j{
tmp=0.0;
for（k=0;ktmp+=rev[i][k]*tmp2[k][j];
lw[i][j]=tmp;
}
for（i=0;i{
tmp=0.0;
for（k=0;ktmp+=*w*rev[i][k]*b[k];
fw[i]=tmp;
}
}
floatTable（intn,floata[N][N],floatb[N],float*w）
{
inti,j,k;
floatlw[N][N],fw[N];
printf（"PleaseinputthematrixAbyrow!
\n"）;
for（i=0;i{
printf（"Row%d:
",i+1）;
k=0;
for（j=0;j{scanf（"%f",&a[i][j]）;
if（a[i][j]==0）
k=k+1;}
if（k==n）
{printf（"thematrixAhasalineof0"）;
exit（0）;}
}
printf（"Pleaseinputthevectorb:
"）;
for（i=0;iscanf（"%f",&b[i]）;
printf（"Inputw:
"）;
scanf（"%f",w）;
cre_sch（n,w,a,b,lw,fw）;
printf（"\nThematrixAandvectorb:
\n"）;
for（i=0;i{
for（j=0;jprintf（"%10.5f",a[i][j]）;
printf（"%10.5f",b[i]）;
printf（"\n"）;
}
printf（"\nTheSORiterativescheme（matrixLw&vectorfw）:
\n"）;
for（i=0;i{
for（j=0;jprintf（"%10.5f",lw[i][j]）;
printf（"%10.5f",fw[i]）;
printf（"\n"）;
}
}
floatinit_vec（intn,floatx[N]）
{
inti;
printf（"\nPleaseinputtheinitialiterationvectorx:
"）;
for（i=0;iscanf（"%f",&x[i]）;
printf（"\nTheinitialiterationvectorx:
\n"）;
for（i=0;iprintf（"%10.5f",x[i]）;
printf（"\n"）;
}
floatsor（intn,floata[N][N],floatb[N],floatx[N],floatw）
{
inti,j,k;
floattmp1,tmp2,x2[N];
for（k=0;;k++）
{
for（i=0;ix2[i]=x[i];
for（i=0;i{
tmp1=0.0;
tmp2=0.0;
for（j=0;jtmp1+=a[i][j]*x[j];
for（j=i;jtmp2+=a[i][j]*x2[j];
x[i]=x2[i]+w*（b[i]-tmp1-tmp2）/a[i][i];
}
for（i=0,j=0;iif（fabs（x2[i]-x[i]）<0.000001）j++;//精度要求,可改变//
if（j==n）
{
printf（"\nThisSORiterativeschemeisconvergent!
\n"）;
printf（"Numberofiterations:
%d",k+1）;
break;
}
if（k==10000）//最大迭代次数，可改变//
{
printf（"ThisSORiterativeschememaybenotconvergent!
"）;
break;
}
}
printf（"\nTheresults:
\n"）;
for（i=0;iprintf（"%12.7f",x[i]）;
}
intmain（）
{
intn;
floatx[N],a[N][N],b[N],w;
printf（"Inputn:
"）;
scanf（"%d",&n）;
Table（n,a,b,&w）;
init_vec（n,x）;
sor（n,a,b,x,w）;
}
高斯LU：
/****************LU分解法*******************/
#defineAL100//局阵最大行
#defineROW100//局阵最大列
#include
#include
#include
#include
voidmain（）{
intn;//阶
doubleA[AL][ROW];//系数矩阵
doubleb[AL];//右端项
inti,j;
charflag='n';//选择标志
//.......文件读入.........
/*cout<<"是否从文件读入（Y/N）:
";
cin>>flag;
if（flag=='y'||flag=='Y'）{
cout<<"请输入文件名:
（file.txt）";
charfile[10];
//已存在文件:
file1.txt,file2.txt
cin>>file;
cout<ifstreamfin（file）;
cout<<"输入阶:
"<fin>>n;
cout<cout<<"输入系数矩阵A:
"<for（i=0;ifor（j=0;jfin>>A[i][j];
cout<}
cout<}
cout<<"输入右端项b:
"<for（i=0;ifin>>b[i];
cout<}
cout<cout<}*/
//........手工输入............
//else{
cout<<"输入阶:
"<cin>>n;
cout<<"输入系数矩阵A:
"<for（i=0;ifor（j=0;jcin>>A[i][j];
cout<<"输入右端项b:
"<for（i=0;icin>>b[i];
cout<//}
//.........计算L矩阵和U矩阵（都放到A矩阵里）..........
intk;
//intm;
doublesum_1=0;
for（k=0;k//L矩阵算法
for（i=k+1;isum_1=0;
for（intm=0;msum_1+=A[i][m]*A[m][k];
A[i][k]=（A[i][k]-sum_1）/A[k][k];
}
//U矩阵算法
for（i=k+1;isum_1=0;
for（intm=0;msum_1+=A[k+1][m]*A[m][i];
A[k+1][i]-=sum_1;
}
}
//.........判断矩阵对角线上是否为0..........
intwr;
for（intnum=0;num{if（A[num][num]==0）
wr=1+wr;
}if（wr==n）
{cout<<"\nLUdecompsitionfailed.";
exit（0）;}
else
cout<<"\n分解过程中对角线上没有0";
//..............LUX=b......................
//第一步:
UX=X_*,算LX_*=b,得到Y（即b[i]）
for（i=0;isum_1=0;
for（j=i;j>0;j--）
sum_1+=A[i][j-1]*b[j-1];
b[i]=b[i]-sum_1;
}
//第二步:
UY=b,得出答案b[i]
for（i=n-1;i>=0;i--）{
sum_1=0;
for（j=i+1;jsum_1=sum_1+A[i][j]*b[j];
}
b[i]=（b[i]-sum_1）/A[i][i];
}
cout<<"\n答案输出:
";
for（i=0;icout<cout<//L矩阵和U矩阵都仍旧存放在A矩阵里
/*for（i=0;ifor（j=0;jcout<cout<}
七、教师评价
升级会员 