☉ 设为首页 | 加入收藏 | 广告服务.
网站首页 C语言 C++语言 Windows编程 接口&驱动 汇编语言 编程资料 软 件 其它文章 等级考试辅导专栏
当前位置:网站首页>>C++语言>>C++基础>>学习园地 双击自动滚屏
矩阵直接三角分解法的c++语言实现
发表日期:2003年10月28日 出处:原创 作者:萧何 已经有2045位读者读过此文
矩阵的计算是非常有用的,也有很多的算法,其中直接三角分解法是比较高效而且误差相对较小的一种算法,我用面向对象的设计思想,将直接三角分解法用C++实现,我的程序还有很多的不足,希望各位给予指教。(直接三角分解法的具体算法 清参照《计算机数值算法》)
#include<iostream.h>
class GS2
{
public:
GS2 (int i,int j,double e) //构造矩阵
{
x=i;
y=j;
eps=e;
a=new double[x*(x+y)];
}
void build();
void return_x(); //回代求 x[i]
bool judge (int); //判断是否|P|<=EPS
bool Doolittle (); //doolittle直接三角分解法的主体程序
void changcemain (int); //选主元
void show(); //打印x
void print();
~GS2()
{
delete []a;
}
protected:
double *a,eps;
int x,y;
};
void GS2::build()
{
for(int i=0;i<x;i++)
for(int j=0;j<x+y;j++)
{
cout<<"x["<<i+1<<"]["<<j+1<<"]=";
cin>>a[j+(x+y)*i];
}
}
//
void GS2::print()
{
for(int i=0;i<x;i++)
for(int j=0;j<x+y;j++)
{
cout<<"x["<<i+1<<"]["<<j+1<<"]="
<<a[j+(x+y)*i]<<endl;
}
}
//
bool GS2::Doolittle()
{
double temp;
for(int r=0;r<x-1;r++)
{
if(r==0)
{
changcemain(r);
if (!judge(r))
return false;
else
for(int i=r+1;i<x;i++)
a[r+(x+y)*i]/=a[r+(x+y)*r];
}
else
{
for(int i=r;i<x;i++)
{
temp=0;
for(int k=0;k<r;k++)
temp+=a[k+(x+y)*i]*a[r+(x+y)*k];
a[r+(x+y)*i]-=temp;
}
changcemain(r);
if (!judge(r))
return false;
for( i=r+1;i<x;i++)
{
a[r+(x+y)*i]/=a[r+(x+y)*r];
}
for(i=r+1;i<x+y;i++)
{
temp=0;
for(int k=0;k<r;k++)
{
temp+=a[k+(x+y)*r]*a[i+(x+y)*k];
}
a[i+(x+y)*r]-=temp;
}
} //end else
}
for(int i=x-1;i<x+y;i++)
{
temp=0;
for(int k=0;k<x-1;k++)
{
temp+=a[k+(x+y)*(x-1)]*a[i+(x+y)*k];
}
a[i+(x+y)*(x-1)]-=temp;
}
if (!judge(x-1))
return false;
return true;
}
void GS2::return_x()
{
double temp;
for(int j=0;j<y;j++)
{
for(int i=x-1;i>=0;i--)
{
temp=0;
for(int k=x-1;k>i;k--)
{
temp+=a[k+(x+y)*i]*a[(x+j)+(x+y)*k];
}
a[(x+j)+(x+y)*i]-=temp;
a[(x+j)+(x+y)*i]/=a[i+(x+y)*i];
}
}
}
void GS2::show()
{
for(int j=0;j<y;j++)
{
cout<<endl;
for(int i=0;i<x;i++)
{
cout<<"x["<<j+1<<"]["<<i+1<<"]="<<a[(x+j)+(x+y)*i]<<'\t';
}
}
}
void GS2::changcemain(int r)
{
double temp=a[r+(x+y)*r];
int t=r;
for(int i=r+1;i<x;i++)
{
if( a[r+(x+y)*i] > temp )
{
temp=a[r+(x+y)*i];
t=i;
}
}
if (t!=r)
for(i=0;i<x+y;i++)
{
temp=a[i+(x+y)*r];
a[i+(x+y)*r]=a[i+(x+y)*t];
a[i+(x+y)*t]=temp;
}
}
bool GS2::judge(int r)
{
if (a[r+(x+y)*r]>0)
return ((a[r+(x+y)*r]>eps )?true:false);
else
return (((-a[r+(x+y)*r])>eps )?true:false);
}
int main()
{
int n,m;
double e;
cout<<"please input n , m and EPS :";
cin>>n>>m>>e;
GS2 gs2(n,m,e);
gs2.build();
if(!gs2.Doolittle())
{
cout<<"此方程无解!"<<endl;
return 1;
}
gs2.print();
gs2.return_x();
gs2.show();
return 0;
}
几个常用的算法,包括消元,迭代,三角分解等,都是解普通线性方程组的
void gauss(double **a,double *b,double *x,int n)
//高斯列主元消去法求解线性方程组
{
int i,j,k,ik;
double mik,temp;
for(k=0;k<n;k++)
{
mik=-1;
for(i=k;i<n;i++)
if(ABS(a[i][k])>mik)
{
mik=ABS(a[i][k]);
ik=i;
}
for(j=k;j<n;j++)
SWAP(a[ik][j],a[k][j]);
SWAP(b[k],b[ik]);
b[k]/=a[k][k];
for(i=n-1;i>=k;i--)
a[k][i]/=a[k][k];
for(i=k+1;i<n;i++)
{
b[i]-=a[i][k]*b[k];
for(j=n-1;j>=k;j--)
a[i][j]-=a[i][k]*a[k][j];
}
}
for(i=n-1;i>=0;i--)
{
x[i]=b[i];
for(j=i+1;j<n;j++)
x[i]-=a[i][j]*x[j];
}
}
void jacobi(double **a,double *b,double *x,int n,int Times)
//雅可比迭代求线性方程组
{
int i,j,t;
double *xx=new double[n];
for(t=0;t<Times;t++)
{
for(i=0;i<n;i++)
{
xx[i]=b[i];
for(j=0;j<n;j++)
if(j!=i)
xx[i]-=a[i][j]*x[j];
xx[i]/=a[i][i];
}
for(i=0;i<n;i++)
x[i]=xx[i];
}
delete xx;
}
void gs(double **a,double *b,double *x,int n,int Times)
//高斯-赛德尔迭代法求解线性方程组
{
int i,j,t;
for(t=0;t<Times;t++)
{
for(i=0;i<n;i++)
{
x[i]=b[i];
for(j=0;j<n;j++)
if(j!=i)
x[i]-=a[i][j]*x[j];
x[i]/=a[i][i];
}
}
}
void doolittle_deco(double **a,int n)
//Doolittle三角分解,分解后的矩阵保存在a中
{
int i,j,k,m;
for(k=0;k<n;k++)
{
for(j=k;j<n;j++)
for(m=0;m<k;m++)
a[k][j]-=a[k][m]*a[m][j];
for(i=k+1;i<n;i++)
{
for(m=0;m<k;m++)
a[i][k]-=a[i][m]*a[m][k];
a[i][k]/=a[k][k];
}
}
}
void doolittle_back(double **a,double *b,double *x,int n)
//Doolittle三角分解法求解线性方程组的回代过程
{
int i,j;
double *y;
y=new double[n];
for(i=0;i<n;i++)
{
y[i]=b[i];
for(j=0;j<i;j++)
y[i]-=a[i][j]*y[j];
}
for(i=n-1;i>=0;i--)
{
x[i]=y[i];
for(j=i+1;j<n;j++)
x[i]-=a[i][j]*x[j];
x[i]/=a[i][i];
}
delete y;
}
void doolittle(double **a,double *b,double *x,int n)
//Doolittle三角分解法求解线性方程组
{
doolittle_deco(a,n);
doolittle_back(a,b,x,n);
}
Top
11 楼saint001(saint001)回复于 2003-12-17 21:34:37 得分 14前面加上
#include "stdio.h"
#include "math.h"
#define ABS(x) (x)>0?(x):-(x)
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
对它们的调用
int main(int argc, char* argv[])
{
int i,j,n;
double **a,*b,*x;
FILE *fp=fopen("d:\\data.txt","r");
fscanf(fp,"%d",&n);
a=new double*[n];
b=new double[n];
x=new double[n];
for(i=0;i<n;i++)
{
a[i]=new double[n];
for(j=0;j<n;j++)
fscanf(fp,"%lf",a[i]+j);
fscanf(fp,"%lf",b+i);
}
fclose(fp);
gauss(a,b,x,n);//Gauss主元消去法
for(i=0;i<n;i++)
printf("%f\t",x[i]);
printf("\n");
for(i=0;i<n;i++)
x[i]=1;
jacobi(a,b,x,3,10);//Jacobi迭代法
for(i=0;i<n;i++)
printf("%f\t",x[i]);
printf("\n");
for(i=0;i<n;i++)
x[i]=1;
gs(a,b,x,n,10);//GS迭代法
for(i=0;i<n;i++)
printf("%f\t",x[i]);
printf("\n");
doolittle(a,b,x,n);//Doolittle三角分解
for(i=0;i<n;i++)
printf("%f\t",x[i]);
printf("\n");
for(i=0;i<n;i++)
delete a[i];
delete a,b,x;
return 0;
}