《数值分析》课程设计.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]*(j
y[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