chap6.h

LDL分解法

#include<iostream.h>
#include<math.h>
#include<fstream.h>
#include<iomanip.h> 

const int N=10;
 
 class matrix
 {
 public:
	 void initialize(int nn, int mm);            //初始化
	 void input_fixed();                         //直接用老师给出的系数矩阵和右端项进行计算
	 void input_client();                        //由用户输入系数矩阵和右端项进行计算
	 void LDL();                                 //LDL分解
	 void solve(int num);                        //求解方程组
	 void output(int num);                       //显示方程组的解
     ~matrix();                                  //析构函数
  
 private :
      
	 int n,m;                                    //系数矩阵的阶数为n，右端项一共有m组


	 double *a;                                  //a存储系数矩阵
	 double *b;                                  //存储全部的右端项
	 double *x;                                  //公式中的X,即方程组的解

	 double *d;                                  //对应于课本上的对角矩阵D
	 double *aa;                                 //对应于课本上的辅助矩阵A~
	 double *l;                                  //对应于课本上的下三角矩阵L
 };

 void matrix::initialize(int nn, int mm)
 { 
	 n=nn;
	 m=mm;

	 a=new double[(n+1)*n/2];      //a存储系数矩阵
	 b=new double[n*m];            //存储全部的右端项
	 x=new double[n];              //公式中的X,即方程组的解
	 
	 d=new double[n];              //对应于课本上的对角矩阵D
	 aa=new double[(n+1)*n/2];     //对应于课本上的辅助矩阵A~
	 l=new double[(n+1)*n/2];      //对应于课本上的下三角矩阵L
 }

 matrix::~matrix()
 {
	 delete[] a;
	 delete[] b;
	 delete[] x;

	 delete[] d;
	 delete[] aa;
	 delete[] l;
 }

 void matrix::input_fixed()
 {
	a[0]=5;
	a[1]=7;   a[2]=10;
	a[3]=6;   a[4]=8;   a[5]=10;
	a[6]=5;   a[7]=7;   a[8]=9;   a[9]=10;
    
	//两组右端项
	b[0]=23;  b[1]=32;  b[2]=33;  b[3]=31;  
	b[4]=92;  b[5]=128; b[6]=132; b[7]=124;
 }

 void matrix::input_client()
 {
	int i=0,j=0;
	cout<<endl<<"请按示例输入为对称正定阵的系数矩阵:"<<endl;
	cout<<"(按行输入，由对称性，只输入矩阵的下三角部分"<<endl;
	cout<<"输入的矩阵必须是正定矩阵，否则用LDL'分解，结果会出错)"<<endl;
	cout<<"示例:"<<endl<<" 5"<<endl<<"-4  6"<<endl<<" 1 -4  6"<<endl<<endl;
	//输出相关提示信息

	cout<<"系数矩阵A如下:"<<endl;//输入系数矩阵
	for(i=0;i<(n+1)*n/2;i++)
	{
		cin>>a[i];
	}

	cout<<endl<<"请按示例输入右端项:(每行输入一组右端项)"<<endl;
	cout<<"示例:"<<endl<<"第1组右端项: 1 3 4"<<endl<<"第2组右端项: 5 1 2"<<endl<<endl;
	//输出相关提示信息

	cout<<"方程组右端项D如下:"<<endl;//输入各组右端项
	for(i=0;i<m;i++)
	{
		cout<<"第"<<i+1<<"组右端项:"<<ends;
		for(j=0;j<n;j++)
		{
			cin>>b[i*n+j];
		}
	}
}


 void matrix::LDL()
 {
	 int i=0,j=0,k=0,s=0;//,m=0,n=0;//j,k,m是用于循环的变量,n用于存放数据在数组中的位置
	 
	 //下面使用的公式参见书本P151
	 d[0]=a[0];
	 l[0]=1;
	 for(j=2;j<n+1;j++)
	 {
		 for(k=1;k<j;k++)
		 {
			 i=j*(j-1)/2+k-1;
			 aa[i]=a[i];
			 for(s=1;s<k;s++)
			 {
				 aa[i]=aa[i]-aa[j*(j-1)/2+s-1]*l[k*(k-1)/2+s-1];
			 } //6.32第一个式子
			 
			 l[i]=aa[i]/d[k-1];//6.32第二个式子
		 }
		 
		 l[j*(j+1)/2-1]=1;
		 d[j-1]=a[(j+1)*j/2-1];
		 for(s=1;s<j;s++)
		 {
			 d[j-1]=d[j-1]-aa[j*(j-1)/2+s-1]*l[j*(j-1)/2+s-1];
		 }//6.32第三个式子
	 }
 }

 void matrix::solve(int num)
 {
	 int i=0,j=0;
	 double *y=new double[n];//公式中的Y
	 //下面先求解下三角方程组LY=b
	 for(i=1;i<=n;i++)
	 {
		 y[i-1]=b[i-1+(num-1)*n];  //其中num代表第num组右端项
		 for(j=1;j<i;j++)
		 {
			 y[i-1]=y[i-1]-y[j-1]*l[i*(i-1)/2+j-1];
		 }
	 }
	 
	 for(i=1;i<=n;i++)
	 {
		 y[i-1]=y[i-1]/d[i-1];//计算公式中Z的各个分量
	 }
	 
	 //求解上三角方程组L'X=Z
	 for(i=n;i>=1;i--)
	 {
		 x[i-1]=y[i-1];
		 for(j=n;j>=i+1;j--)
		 {
			 x[i-1]=x[i-1]-x[j-1]*l[j*(j-1)/2+i-1];
		 }
	 }
	 
	 
	 delete[] y;
 }

 void matrix::output(int num)
 {
	 int i=0,j=0;
	 cout<<"对于系数矩阵为A的方程组"<<endl;
	 //将系数矩阵显示出来
	 cout<<setiosflags(ios::right);
	 cout<<"A = [ ";
	 cout<<setw(2)<<a[0]<<endl;

	 for(i=1;i<n;i++)
	 {
		 for(j=0;j<i+1;j++)
		 {
			 cout<<setw(8)<<a[(i+1)*i/2+j];
			 if(j==n-1)
				 cout<<" ]";
		 }
		 cout<<endl;
	 }
	 cout<<"(由对称性，只显示下三角的元素)"<<endl<<endl;
	 
	 
	 cout<<"当右端项 b = [ ";//将右端项现实出来
	 for(j=0;j<n;j++)
	 {
		 cout<<b[j+(num-1)*n]<<ends<<ends;
	 }
	 cout<<"]'"<<endl;
	 
	 //将方程组的解显示出来
	 cout<<"方程组解为:"<<endl<<"X = [ ";
	 for(j=0;j<n;j++)
	 {
		 cout<<x[j]<<ends<<ends;
	 }
	 cout<<"]'"<<endl<<endl;
	 
 }