计算方法.docx
《计算方法.docx》由会员分享，可在线阅读，更多相关《计算方法.docx（25页珍藏版）》请在冰点文库上搜索。
计算方法
Cramer算法解方程组
Gauss列主元解方程组
Gauss全主元解方程组
用Doolittle算法解方程组
//解线性方程组
#include
#include
#include
//----------------------------------------------全局变量定义区
constintNumber=15;//方程最大个数
doublea[Number][Number],b[Number],copy_a[Number][Number],copy_b[Number];//系数行列式
intA_y[Number];//a[][]中随着横坐标增加列坐标的排列顺序,如a[0][0],a[1][2],a[2][1]...则A_y[]={0,2,1...};
intlenth,copy_lenth;//方程的个数
doublea_sum;//计算行列式的值
char*x;//未知量a,b,c的载体
//----------------------------------------------函数声明区
voidinput（）;//输入方程组
voidprint_menu（）;//打印主菜单
intchoose（）;//输入选择
voidcramer（）;//Cramer算法解方程组
voidgauss_row（）;//Gauss列主元解方程组
voidguass_all（）;//Gauss全主元解方程组
voidDoolittle（）;//用Doolittle算法解方程组
intDoolittle_check（doublea[][Number],doubleb[Number]）;//判断是否行列式>0,若是,调整为顺序主子式全>0
voidxiaoqu_u_l（）;//将行列式Doolittle分解
voidcalculate_u_l（）;//计算Doolittle结果
double&calculate_A（intn,intm）;//计算行列式
doublequanpailie_A（）;//根据列坐标的排列计算的值,如A_y[]={0,2,1},得sum=a[0][A_y[0]]*a[1][A_y[1]]*a[2][A_y[2]]=a[0][0]*a[1][2]*a[2][1];
voidexchange（intm,inti）;//交换A_y[m],A_y[i]
voidexchange_lie（intj）;//交换a[][j]与b[];
voidexchange_hang（intm,intn）;//分别交换a[][]和b[]中的m与n两行
voidgauss_row_xiaoqu（）;//Gauss列主元消去法
voidgauss_all_xiaoqu（）;//Gauss全主元消去法
voidgauss_calculate（）;//根据Gauss消去法结果计算未知量的值
voidexchange_a_lie（intm,intn）;//交换a[][]中的m和n列
voidexchange_x（intm,intn）;//交换x[]中的x[m]和x[n]
voidrecovery（）;//恢复数据
//主函数
voidmain（）
{
intflag=1;
input（）;//输入方程
while（flag）
{
print_menu（）;//打印主菜单
flag=choose（）;//选择解答方式
}
}
//函数定义区
voidprint_menu（）
{
system（"cls"）;
cout<<"------------方程系数和常数矩阵表示如下:
\n";
for（intj=0;jcout<<"系数"<cout<<"\t常数";
cout<for（inti=0;i{
for（j=0;jcout<:
left）<cout<<"\t"<}
cout<<"-----------请选择方程解答的方案----------";
cout<<"\n1.克拉默（Cramer）法则";
cout<<"\n2.Gauss列主元消去法";
cout<<"\n3.Gauss全主元消去法";
cout<<"\n4.Doolittle分解法";
cout<<"\n5.退出";
cout<<"\n输入你的选择:
";
}
voidinput（）
{inti,j;
cout<<"方程的个数:
";
cin>>lenth;
if（lenth>Number）
{
cout<<"Itistoobig.\n";
return;
}
x=newchar[lenth];
for（i=0;ix[i]='a'+i;
//输入方程矩阵
//提示如何输入
cout<<"====================================================\n";
cout<<"请在每个方程里输入"<\n";
cout<<"例：
\n方程:
a";
for（i=1;i{
cout<<"+"<}
cout<<"=10\n";
cout<<"应输入:
";
for（i=0;icout<cout<<"10\n";
cout<<"==============================\n";
//输入每个方程
for（i=0;i{
cout<<"输入方程"<";
for（j=0;jcin>>a[i][j];
cin>>b[i];
}
//备份数据
for（i=0;ifor（j=0;jcopy_a[i][j]=a[i][j];
for（i=0;icopy_b[i]=b[i];
copy_lenth=lenth;
}
//输入选择
intchoose（）
{
intchoice;charch;
cin>>choice;
switch（choice）
{
case1:
cramer（）;break;
case2:
gauss_row（）;break;
case3:
guass_all（）;break;
case4:
Doolittle（）;break;
case5:
return0;
default:
cout<<"输入错误,请重新输入:
";
choose（）;
break;
}
cout<<"\n是否换种方法求解（Y/N）:
";
cin>>ch;
if（ch=='n'||ch=='N'）return0;
recovery（）;
cout<<"\n\n\n";
return1;
}
//用克拉默法则求解方程.
voidcramer（）
{
inti,j;doublesum,sum_x;charch;
//令第i行的列坐标为i
cout<<"用克拉默（Cramer）法则结果如下:
\n";
for（i=0;iA_y[i]=i;
sum=calculate_A（lenth,0）;
if（sum!
=0）
{
cout<<"系数行列式不为零,方程有唯一的解:
";
for（i=0;i{ch='a'+i;
a_sum=0;
for（j=0;jA_y[j]=j;
exchange_lie（i）;
sum_x=calculate_A（lenth,0）;
cout<exchange_lie（i）;
}
}
else
{
cout<<"系数行列式等于零,方程没有唯一的解.";
}
cout<<"\n";
}
double&calculate_A（intn,intm）//计算行列式
{inti;
if（n==1）{
a_sum+=quanpailie_A（）;
}
else{for（i=0;i{exchange（m,m+i）;
calculate_A（n-1,m+1）;
exchange（m,m+i）;
}
}
returna_sum;
}
doublequanpailie_A（）//计算行列式中一种全排列的值
{
inti,j,l;
doublesum=0,p;
for（i=0,l=0;ifor（j=0;A_y[j]!
=i&&jif（A_y[j]>i）l++;
for（p=1,i=0;ip*=a[i][A_y[i]];
sum+=p*（（l%2==0）?
（1）:
（-1））;
returnsum;
}
//高斯列主元排列求解方程
voidgauss_row（）
{
inti,j;
gauss_row_xiaoqu（）;//用高斯列主元消区法将系数矩阵变成一个上三角矩阵
for（i=0;i{
for（j=0;jcout<cout<}
if（a[lenth-1][lenth-1]!
=0）
{
cout<<"系数行列式不为零,方程有唯一的解：
\n";
gauss_calculate（）;
for（i=0;i{
cout<}
}
else
cout<<"系数行列式等于零,方程没有唯一的解.\n";
}
voidgauss_row_xiaoqu（）//高斯列主元消去法
{
inti,j,k,maxi;doublelik;
cout<<"用Gauss列主元消去法结果如下:
\n";
for（k=0;k{
j=k;
for（maxi=i=k;iif（a[i][j]>a[maxi][j]）maxi=i;
if（maxi!
=k）
exchange_hang（k,maxi）;//
for（i=k+1;i{
lik=a[i][k]/a[k][k];
for（j=k;ja[i][j]=a[i][j]-a[k][j]*lik;
b[i]=b[i]-b[k]*lik;
}
}
}
//高斯全主元排列求解方程
voidguass_all（）
{
inti,j;
gauss_all_xiaoqu（）;
for（i=0;i{
for（j=0;jcout<cout<}
if（a[lenth-1][lenth-1]!
=0）
{
cout<<"系数行列式不为零,方程有唯一的解：
\n";
gauss_calculate（）;
for（i=0;i{
for（j=0;x[j]!
='a'+i&&jcout<}
}
else
cout<<"系数行列式等于零,方程没有唯一的解.\n";
}
voidgauss_all_xiaoqu（）//Gauss全主元消去法
{
inti,j,k,maxi,maxj;doublelik;
cout<<"用Gauss全主元消去法结果如下:
\n";
for（k=0;k{
for（maxi=maxj=i=k;i{
for（j=k;jif（a[i][j]>a[maxi][maxj]）
{maxi=i;
maxj=j;
}
}
if（maxi!
=k）
exchange_hang（k,maxi）;
if（maxj!
=k）
{
exchange_a_lie（maxj,k）;//交换两列
exchange_x（maxj,k）;
}
for（i=k+1;i{
lik=a[i][k]/a[k][k];
for（j=k;ja[i][j]=a[i][j]-a[k][j]*lik;
b[i]=b[i]-b[k]*lik;
}
}
}
voidgauss_calculate（）//高斯消去法以后计算未知量的结果
{
inti,j;doublesum_ax;
b[lenth-1]=b[lenth-1]/a[lenth-1][lenth-1];
for（i=lenth-2;i>=0;i--）
{
for（j=i+1,sum_ax=0;jsum_ax+=a[i][j]*b[j];
b[i]=（b[i]-sum_ax）/a[i][i];
}
}
voidDoolittle（）//Doolittle消去法计算方程组
{
doubletemp_a[Number][Number],temp_b[Number];inti,j,flag;
for（i=0;ifor（j=0;jtemp_a[i][j]=a[i][j];
flag=Doolittle_check（temp_a,temp_b）;
if（flag==0）cout<<"\n行列式为零.无法用Doolittle求解.";
xiaoqu_u_l（）;
calculate_u_l（）;
cout<<"用Doolittle方法求得结果如下:
\n";
for（i=0;i{
for（j=0;x[j]!
='a'+i&&jcout<}
}
voidcalculate_u_l（）//计算Doolittle结果
{inti,j;doublesum_ax=0;
for（i=0;i{
for（j=0,sum_ax=0;jsum_ax+=a[i][j]*b[j];
b[i]=b[i]-sum_ax;
}
for（i=lenth-1;i>=0;i--）
{
for（j=i+1,sum_ax=0;jsum_ax+=a[i][j]*b[j];
b[i]=（b[i]-sum_ax）/a[i][i];
}
}
voidxiaoqu_u_l（）//将行列式按Doolittle分解
{inti,j,n,k;doubletemp;
for（i=1,j=0;ia[i][j]=a[i][j]/a[0][0];
for（n=1;n{//求第n+1层的上三角矩阵部分即U
for（j=n;j{for（k=0,temp=0;ktemp+=a[n][k]*a[k][j];
a[n][j]-=temp;
}
for（i=n+1;i{for（k=0,temp=0;ktemp+=a[i][k]*a[k][n];
a[i][n]=（a[i][n]-temp）/a[n][n];
}
}
}
intDoolittle_check（doubletemp_a[][Number],doubletemp_b[Number]）//若行列式不为零,将系数矩阵调整为顺序主子式大于零
{
inti,j,k,maxi;doublelik,temp;
for（k=0;k{
j=k;
for（maxi=i=k;iif（temp_a[i][j]>temp_a[maxi][j]）maxi=i;
if（maxi!
=k）
{exchange_hang（k,maxi）;
for（j=0;j{temp=temp_a[k][j];
temp_a[k][j]=temp_a[maxi][j];
temp_a[maxi][j]=temp;
}
}
for（i=k+1;i{
lik=temp_a[i][k]/temp_a[k][k];
for（j=k;jtemp_a[i][j]=temp_a[i][j]-temp_a[k][j]*lik;
temp_b[i]=temp_b[i]-temp_b[k]*lik;
}
}
if（temp_a[lenth-1][lenth-1]==0）return0;
return1;
}
voidexchange_hang（intm,intn）//交换a[][]中和b[]两行
{
intj;doubletemp;
for（j=0;j{temp=a[m][j];
a[m][j]=a[n][j];
a[n][j]=temp;
}
temp=b[m];
b[m]=b[n];
b[n]=temp;
}
voidexchange（intm,inti）//交换A_y[m],A_y[i]
{inttemp;
temp=A_y[m];
A_y[m]=A_y[i];
A_y[i]=temp;
}
voidexchange_lie（intj）//交换未知量b[]和第i列
{doubletemp;inti;
for（i=0;i{temp=a[i][j];
a[i][j]=b[i];
b[i]=temp;
}
}
voidexchange_a_lie（intm,intn）//交换a[]中的两列
{doubletemp;inti;
for（i=0;i{temp=a[i][m];
a[i][m]=a[i][n];
a[i][n]=temp;
}
}
voidexchange_x（intm,intn）//交换未知量x[m]与x[n]
{chartemp;
temp=x[m];
x[m]=x[n];
x[n]=temp;
}
voidrecovery（）//用其中一种方法求解后恢复数据以便用其他方法求解
{
for（inti=0;ifor（intj=0;ja[i][j]=copy_a[i][j];
for（i=0;ib[i]=copy_b[i];
for（i=0;ix[i]='a'+i;
a_sum=0;
lenth=copy_lenth;
}
胡哲光
（台州职业技术学院计算机工程系，浙江 台州 318000）
【摘 要】高斯消去法可以在没有舍入误差影响的条件下经过有限步的四则运算求得线性代数方程的精确解，是日前计算机上常用于求解低阶稠密矩阵方程组的有效方法。
文章主要就在C＋＋中实现带主元选择的高斯消去法求解N阶线性代数方程进行了讨论。
【关键词】高斯消去法；线性方程；矩阵C＋＋
【中图分类号】TP316 【文献标识码】A 【文章编号】1008-1151（2006）08-0044-02
解线性代数方程是科学研究和工程计算中经常遇到的问题。
对于n阶线性代数方程组：
（1）
假定该n阶线性方程有唯一解，则如果把它转换成三角形方程组，例如：
则可以采取逆推方式进行，依次求得xn、xn-1、……、x1。
该求解n阶线性方程方法就是高斯消去法，共分成两个过程：
将线性方程组转换成三角形方程组的过程称为消元过程；逆反求解的过程称为回代过程。
一、高斯消去法算法
若将线性方程组
（1）用矩阵表示则可写为：
AX=b。
其中A＝[aij]n×n，是n阶非奇异的系数矩阵，X=（x1，x2，……，xn）T是末知向量，b＝（b1，b2，……，bn）T是右端向量。
将n阶矩阵A和右端向量b合成一个矩阵
1.对该矩阵进行消元计算，公式如下：
对k＝1，2，……，n－1依次计算
消元后原矩阵成为三角矩阵：
2.逆反矩阵
当aii≠0（I=1，2，……，n）时方程组有唯一解，求解过程可以从最后一个方程入手，采用逆推方式进行，公式如下：
二、选主元技巧
以上高斯消去法是按照方程及末知数的给定顺序依次进行消元，但如果在消元
升级会员 