《数值分析》课程设计.docx
《《数值分析》课程设计.docx》由会员分享，可在线阅读，更多相关《《数值分析》课程设计.docx（30页珍藏版）》请在冰点文库上搜索。
《数值分析》课程设计
问题的提出
3.3用SOR方法解下列方程组（去松弛因子w=1.2），要求
。
3.4设
，
计算
。
3.5用选列主元Gauss消元法求解方程组
3.6用追赶法解三对角方程组
3.7用三角分解法解方程组
3.8用选主元消元法计算下列行列式
。
一、问题分析
1.超松弛法是迭代方法的一种加速方法，其计算公式简单，但需要选择合适的松弛因子，以保证迭代过程有较快的收敛速度。
2.A的条件数计算首先要获得A的逆，而求A的逆可以转化为求n个方程组。
3.完全主元消元法在计算过程中花费了大量的时间用于寻找主元。
同时，各变量的位置在消元过程中也可能会发生变化。
而列选主元法则可消除这个弊病。
4.追赶法主要是解三对角方程组。
所谓追指消元过程，赶指回代过程。
5.Gauss消元法是通过逐步消元过程，将方程组的系数矩阵A转变为一个上三角矩阵。
三角分解法，就是把系数矩阵分解为两个三角阵。
6.将某一向量坐标同乘以某非零实数，加到另一向量上，行列式的值不变。
用选主元法将行列式矩阵变为三角阵，对角线上的数值相乘即为行列式的值。
二、编程解决
3.3Sor法c语言编程：
#include
#include
#include
#defineomega1.2//取值不合适结果可能发散
voidmain（）
{
doublea[5][5];
doubleb[5],x[5],f,t,y[5]={0,0,0,0,0};
inti,j,n,cnt=0;
printf（"阶数:
"）;
scanf（"%d",&n）;
printf（"请输入%d阶的A矩阵\n",n）;
for（i=0;ifor（j=0;jscanf（"%lf",&a[i][j]）;
printf（"请输入B矩阵\n"）;
for（i=0;iscanf（"%lf",&b[i]）;
printf（"count\t"）;
for（i=0;iprintf（"x[%d]\t\t",i）;
printf（"收敛程度\n"）;
do
{
for（i=0;ix[i]=y[i];
for（i=0;i{
t=0;
for（j=0;jt=t+a[i][j]*（jy[j]:
x[j]）;
y[i]=x[i]+omega*（b[i]-t）/a[i][i];
}
printf（"%d",cnt++）;
for（i=0;iprintf（"\t%lf",x[i]）;
f=0;
for（i=0;if+=fabs（y[i]-x[i]）;
printf（"\t%g\n",f）;
}while（f>1e-4&&cnt<100）;
}
所得结果：
3.4求逆、算条件数编程：
#include
#include
#include
#defineN5//可修改，以改变可解决的最大维数。
但会浪费更多的空间
doublea[N][N],b[N][N],c[N][N];
intn,flag=1;
voidinput（）
{
inti,j;
do
{
printf（"阶数:
"）;
scanf（"%d",&n）;
if（n>N||n<=0）//宏定义中N规定了数组的最大维数，超过则出错
printf（"溢出，请重新输入!
\n"）;
}while（n>N||n<=0）;
printf（"请输入%d阶的A矩阵\n",n）;
for（i=0;ifor（j=0;j{
if（i==j）
b[i][j]=1;
else
b[i][j]=0;
scanf（"%lf",&a[i][j]）;
c[i][j]=a[i][j];
}
}
voidwatch（）
{
inti,j;
for（i=0;i{
for（j=0;jprintf（"%4.2lf\t",a[i][j]）;
for（j=0;jprintf（"%4.2lf\t",b[i][j]）;
printf（"\n"）;
}
system（"pause"）;
}
voidswap（double*p,double*q）
{
doubletmp;
tmp=*p;
*p=*q;
*q=tmp;
}
voidswap_x（ints,intt）//S<==>T
{
intj;
for（j=0;j{
swap（&a[s][j],&a[t][j]）;
swap（&b[s][j],&b[t][j]）;
}
}
voidself_x（inti,doubles）//I=I/s
{
intj;
for（j=0;j{
a[i][j]=a[i][j]/s;
b[i][j]=b[i][j]/s;
}
}
voidinter_x（inti,intr,doublet）//I=I-R*t
{
intj;
for（j=0;j{
a[i][j]=a[i][j]-a[r][j]*t;
b[i][j]=b[i][j]-b[r][j]*t;
}
}
voidsort（intk）
{
inti,maxi=k;
doublemax;
max=fabs（a[k][k]）;//进入循环前先假定首数最大
for（i=k;iif（fabs（a[i][k]）>max）
{
maxi=i;
max=fabs（a[i][k]）;
}
/*将最大的一个数放到首数位置*/
if（maxi!
=k）
swap_x（maxi,k）;
}
voidinv（）
{
inti,k;
for（k=0;k{
sort（k）;//消元之前先排序
if（a[k][k]==0）
{
flag=0;
return;
}
self_x（k,a[k][k]）;
for（i=k+1;iinter_x（i,k,a[i][k]）;
watch（）;
}
for（k=n-1;k>0;k--）//向上消元过程
{
for（i=0;iinter_x（i,k,a[i][k]）;
watch（）;
}
}
voidoutput（）
{
inti,j;
if（flag）
{
printf（"A的逆矩阵如下:
\n"）;
for（i=0;i{
for（j=0;jprintf（"%4.2lf\t",b[i][j]）;
printf（"\n"）;
}
}
else
printf（"A矩阵不可逆\n"）;
}
doublecond（）//求条件数
{
doubleu[N],v[N],t,max1,max2;
inti,j;
for（i=0;i{
t=0;
for（j=0;j{
t+=fabs（b[i][j]）;
}
u[i]=t;
}
for（i=0;i{
t=0;
for（j=0;j{
t+=fabs（c[i][j]）;
}
v[i]=t;
}
max1=u[0];
max2=v[0];
for（i=0;i{
if（u[i]>max1）
max1=u[i];
if（v[i]>max2）
max2=v[i];
}
return（max1*max2）;
}
voidmain（）
{doublecond（）;
input（）;
inv（）;//求逆矩阵
output（）;
printf（"thecondis:
%lf",cond（））;//求条件数
printf（"\n"）;
}
输出结果：
3.5选列主元gauss消元法编程：
#include
#include
#include
#defineN5
doublea[N][N];
doubleb[N],x[N];
intn,flag=1;
voidinput（）
{
inti,j;
do
{
printf（"阶数:
"）;
scanf（"%d",&n）;
if（n>N）
printf（"溢出，请重新输入!
\n"）;
}while（n>N）;
printf（"请输入%d阶的A矩阵\n",n）;
for（i=0;ifor（j=0;jscanf（"%lf",&a[i][j]）;
printf（"请输入B矩阵\n"）;
for（i=0;iscanf（"%lf",&b[i]）;
}
voidswap（double*p,double*q）
{
doubletmp;
tmp=*p;
*p=*q;
*q=tmp;
}
voidswap_x（ints,intt）
{
inti;
for（i=t;iswap（&a[s][i],&a[t][i]）;
swap（&b[s],&b[t]）;
}
voidsort（intk）
{
inti,maxi=k;
doublemax;
max=fabs（a[k][k]）;//进入循环前先假定首数最大
for（i=k;iif（fabs（a[i][k]）>max）
{
maxi=i;
max=fabs（a[i][k]）;
}
/*将最大的一个数放到首数位置*/
if（maxi!
=k）
swap_x（maxi,k）;
}
voidgauss（）
{
inti,j,k;
doublet;
for（k=0;k{
sort（k）;//消元之前先排序
if（a[k][k]==0）//若发现排序后首数为0，终止!
{
flag=0;
break;
}
for（i=k+1;i{
t=a[i][k]/a[k][k];
a[i][k]=0;
for（j=k+1;ja[i][j]=a[i][j]-a[k][j]*t;
b[i]=b[i]-b[k]*t;
}
for（i=0;i{
for（j=0;jprintf（"%6.3lf\t",a[i][j]）;
printf（"%6.3lf\n",b[i]）;
}
system（"pause"）;
}
if（a[k][k]==0）//回代之前发现系数为0，终止!
flag=0;
else
{
x[n-1]=b[n-1]/a[n-1][n-1];//回代
for（i=n-2;i>=0;i--）
{
t=0;
for（j=i+1;jt+=a[i][j]*x[j];
x[i]=（b[i]-t）/a[i][i];
}
}
}
voidoutput（）
{
inti;
if（flag）
for（i=0;iprintf（"x[%d]=%6.3lf\t",i+1,x[i]）;
else
printf（"消元过程被强行终止!
"）;
printf（"\n"）;
}
voidmain（）
{
input（）;
gauss（）;
output（）;
}
输出结果：
3.6追赶法编程：
#include
#include
#defineN5
intn;
doublea[N],b[N],c[N],d[N],u[N],q[N],x[N];
voidinput（）
{
inti;
do
{
printf（"阶数:
"）;
scanf（"%d",&n）;
if（n>N）
printf（"溢出，请重新输入!
\n"）;
}while（n>N）;
printf（"a="）;
for（i=1;iscanf（"%lf",&a[i]）;
printf（"b="）;
for（i=0;iscanf（"%lf",&b[i]）;
printf（"c="）;
for（i=0;iscanf（"%lf",&c[i]）;
printf（"d="）;
for（i=0;iscanf（"%lf",&d[i]）;
}
voidchase（）
{
inti,j;
//变换成上三角矩阵
u[0]=c[0]/b[0];
q[0]=d[0]/b[0];
for（i=1;i{
u[i]=c[i]/（b[i]-u[i-1]*a[i]）;
q[i]=（d[i]-q[i-1]*a[i]）/（b[i]-u[i-1]*a[i]）;
}
for（i=0;i{
for（j=0;j{
if（ji+1）
printf（"0\t"）;
elseif（i==j）
printf（"1\t"）;
else
printf（"%4.2lf\t",u[i]）;
}
printf（"%lf\n",q[i]）;
}
//回代过程
x[n-1]=q[n-1];
for（i=n-2;i>=0;i--）
x[i]=q[i]-u[i]*x[i+1];
}
voidoutput（）
{
inti;
for（i=0;iprintf（"x[%d]=%6.3lf\t",i,x[i]）;
}
voidmain（）
{
input（）;
chase（）;
output（）;
}
输出结果：
3.7三角分解法解方程组
程序：
#include
#include
#defineN6//定义一个维数N，以应用于6维以下
doublea[N][N],b[N],X[N],Y[N],u[N][N],v[N][N],c[N][N];
intn;
voidinput（）//输入a系数矩阵，b常数矩阵
{
inti,j;
printf（"pleaseinputthedimension（lessthan7）:
\n"）;//输入维数
scanf（"%d",&n）;
printf（"pleaseinputthematrixa[n][n]:
\n"）;
for（i=0;ifor（j=0;jscanf（"%lf",&a[i][j]）;
printf（"pleaseinputthematrixb[n]:
\n"）;
for（i=0;iscanf（"%lf",&b[i]）;
}
voidsolveUandV（）//求解u和v三角矩阵
{
inti,j,k;
doubletemp1,temp2;
for（i=n-1;i>=0;i--）
{
for（j=i;j>=0;j--）
{
temp1=0;
for（k=i+1;k{
temp1+=u[j][k]*v[k][i];
}
u[j][i]=a[j][i]-temp1;
temp2=0;
for（k=i+1;k{
temp2+=u[i][k]*v[k][j];
}
v[i][j]=（a[i][j]-temp2）/u[i][i];
}
}
printf（"上三角u:
\n"）;
for（i=0;i{
for（j=0;jprintf（"%lf\t",u[i][j]）;
printf（"\n"）;
}
printf（"下三角v:
\n"）;
for（i=0;i{
for（j=0;jprintf（"%lf\t",v[i][j]）;
printf（"\n"）;
}
}
voidsubstituting（）//回代
{
inti,j;
doubletemp1,temp2;
printf（"the中间解Yis:
\n"）;
for（i=n-1;i>=0;i--）//回代求出Y矩阵的值
{
temp1=0;
for（j=i+1;jtemp1+=Y[j]*u[i][j];
Y[i]=（b[i]-temp1）/u[i][i];
printf（"%lf\t\t",Y[i]）;
}
printf（"\n"）;
for（i=0;i{
temp2=0;
for（j=0;j{
temp2+=X[j]*v[i][j];
}
X[i]=Y[i]-temp2;
}
}
voidverification（）//验证拆分的正确性
{
inti,j,k;
doublet;
for（i=0;i{
for（j=0;j{
t=0;
for（k=0;kt+=u[i][k]*v[k][j];
c[i][j]=t;
}
}
printf（"the验证的a1is:
\n"）;//输出验证得到的矩阵
for（i=0;i{
for（j=0;jprintf（"%lf\t",c[i][j]）;
printf（"\n"）;
}
}
voidoutput（）//输出解x
{
inti;
printf（"thevalueoftheseequationsis:
\n"）;
for（i=0;iprintf（"x[%d]=%lf\n",i,X[i]）;
}
voidmain（）
{
input（）;//输入
solveUandV（）;//拆分a为u、v上三角和下三角矩阵
substituting（）;//回代求出y（中间解）和x
verification（）;//验证拆分的正确性
output（）;//输出解
}
结果：
3.8选主元消元法编程：
#include
#include
#include
#defineN5//可修改，以改变可解决的最大维数。
但会浪费更多的空间
doublea[N][N];
intn,flag=1;
voidinput（）
{
inti,j;
do
{
printf（"阶数:
"）;
scanf（"%d",&n）;
if（n>N）//宏定义中N规定了数组的最大维数，超过则出错
printf（"溢出，请重新输入!
\n"）;
}while（n>N）;
printf（"请输入%d阶的A矩阵\n",n）;
for（i=0;ifor（j=0;jscanf（"%lf",&a[i][j]）;
}
voidswap（double*p,double*q）
{
doubletmp;
tmp=*p;
*p=*q;
*q=tmp;
}
voidswap_x（ints,intt）
{
inti;
for（i=t;iswap（&a[s][i],&a[t][i]）;
}
voidswap_y（ints,intt）
{
inti,tmp;
for（i=0;iswap（&a[i][s],&a[i][t]）;
}
voidsort（intk）
{
inti,j,maxi=k,maxj=k;
doublemax;
max=fabs（a[k][k]）;//进入循环前先假定首数最大
for（i=k;ifor（j=k;jif（fabs（a[i][j]）>max）
{
maxi=i;
maxj=j;
max=fabs（a[i][j]）;
}
/*将最大的一个数放到首数位置*/
if（maxi!
=k）
swap_x（maxi,k）;
if（maxj!
=k）
swap_y（maxj,k）;
}
voidgauss（）
{
inti,j,k;
doublet;
for（k=0;k
升级会员 