math.cpp

代数方程

// Math.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "stdio.h"
#include "math.h"

#define ABS(x) (x)>0?(x):-(x)
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}

void mult(double **a,double *x,double *y,int n)
//矩阵与向量相乘
{
	int i,j;
	for(i=0;i<n;i++)
	{
		y[i]=0;
		for(j=0;j<n;j++)
			y[i]+=a[i][j]*x[j];
	}
}

void gauss(double **a,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]);
		for(i=n-1;i>=k;i--)
			a[k][i]/=a[k][k];
		for(i=k+1;i<n;i++)
		{
			for(j=n-1;j>=k;j--)
				a[i][j]-=a[i][k]*a[k][j];
		}
	}
}

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);
}

void reverse(double **a,double **r,int n)
//矩阵求逆
{
	int i,j;
	double *b,*x;
	b=new double[n];
	x=new double[n];
	doolittle_deco(a,n);
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			b[j]=0;
		b[i]=1;
		doolittle_back(a,b,x,n);
		for(j=0;j<n;j++)
			r[j][i]=x[j];
	}
	delete b,x;
}

void linear(double **x,double *y,double *beta,int n,int p)
//最小二乘法多元回归
{
	double **a,*b;
	int i,j,k;
	a=new double*[p];
	for(i=0;i<p;i++)
		a[i]=new double[p];
	for(i=0;i<p;i++)
		for(j=0;j<p;j++)
		{
			a[i][j]=0;
			for(k=0;k<n;k++)
				a[i][j]+=x[k][i]*x[k][j];
		}
	b=new double[p];
	for(i=0;i<p;i++)
	{
		b[i]=0;
		for(j=0;j<n;j++)
			b[i]+=x[j][i]*y[j];
	}
	gauss(a,b,beta,p);
	for(i=0;i<p;i++)
		delete a[i];
	delete a,b;
}

void polyfit(double *x,double *y,double *beta,int n,int p)
//一元多项式拟合
{
	int i,j;
	double **xx;
	xx=new double*[n];
	for(i=0;i<n;i++)
	{
		xx[i]=new double[p];
		for(j=0;j<p;j++)
			xx[i][j]=pow(x[i],j);
	}
	linear(xx,y,beta,n,p);
	for(i=0;i<n;i++)
		delete xx[i];
	delete xx;
}

void polyfitw(double *x,double *y,double *w,double *beta,int n,int p)
//一元多项式加权拟合
{
	int i,j;
	double **xx;
	xx=new double*[n];
	for(i=0;i<n;i++)
	{
		w[i]=sqrt(w[i]);
		xx[i]=new double[p];
		for(j=0;j<p;j++)
			xx[i][j]=pow(x[i],j)*w[i];
		y[i]*=w[i];
	}
	linear(xx,y,beta,n,p);
	for(i=0;i<n;i++)
		delete xx[i];
	delete xx;
}

void cholesky(double **a,double **L,int n)
//Cholesky分解，得到的L是下三角矩阵
{
	int k,i,m;
	for(k=0;k<n;k++)
	{
		L[k][k]=a[k][k];
		for(m=0;m<k;m++)
			L[k][k]-=L[k][m]*L[k][m];
		L[k][k]=sqrt(L[k][k]);
		for(i=k+1;i<n;i++)
		{
			L[i][k]=a[i][k];
			for(m=0;m<k;m++)
				L[i][k]-=L[i][m]*L[k][m];
			L[i][k]/=L[k][k];
		}
	}
	for(i=0;i<n;i++)
		for(m=i+1;m<n;m++)
			L[i][m]=0;
}

void reverse_L(double **L,double **R,int n)
{
	int i,k,m;
	for(k=0;k<n;k++)
	{
		R[k][k]=1/L[k][k];
		for(i=k-1;i>=0;i--)
		{
			R[k][i]=0;
			for(m=i+1;m<=k;m++)
				R[k][i]-=R[k][m]*L[m][i];
			R[k][i]/=L[i][i];
		}
	}
	for(i=0;i<n;i++)
		for(m=i+1;m<n;m++)
			R[i][m]=0;
}

void conjugate(double **a,double *b,double *x,int n,double eps)
{
	double *p,*r,alfa,beta,*ap,temp;
	p=new double[n];
	r=new double[n];
	ap=new double[n];
	int i,j,t;
	for(i=0;i<n;i++)
	{
		p[i]=b[i];
		for(j=0;j<n;j++)
			p[i]-=a[i][j]*x[j];
		r[i]=p[i];
	}
	for(t=0;t<3;t++)
	{
		for(temp=0,i=0;i<n;i++)
			temp+=r[i]>0?r[i]:-r[i];
		if(temp<eps)
		{
			printf("Times:%d\n",t);
			break;
		}
		for(i=0;i<n;i++)
		{
			ap[i]=0;
			for(j=0;j<n;j++)
				ap[i]+=a[i][j]*p[j];
		}
		for(alfa=0,i=0;i<n;i++)
			alfa+=r[i]*p[i];
		for(temp=0,i=0;i<n;i++)
			temp+=ap[i]*p[i];
		alfa/=temp;
		for(i=0;i<n;i++)
			x[i]+=alfa*p[i];
		for(i=0;i<n;i++)
		{
			r[i]=b[i];
			for(j=0;j<n;j++)
				r[i]-=a[i][j]*x[j];
		}
		for(beta=0,i=0;i<n;i++)
			beta-=r[i]*ap[i];
		beta/=temp;
		for(i=0;i<n;i++)
			p[i]=r[i]+beta*p[i];
	}
}

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++)
		x[i]=0;
	conjugate(a,b,x,n,1e-2);//共轭梯度法
	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;
}