yang.cpp

数学

// suanfa.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"


#include <stdio.h>
#include <math.h>

//BEGIN_A放置初始距阵，B放b。然后程序依次输出列选高斯，全选高斯，HOUSEHOLDER变换，QR求特征值

double begin_a[5][5]={{10,1,2,3,4},{1,9,-1,2,-3},{2,-1,7,3,-5},{3,2,3,12,-1},{4,-3,-5,-1,15}};
double a[5][5];
double q[5][5];
double b[5] = {0,0,18,24,31};

int n=5;

void init_a()
{
	int i,j;
	for (i=0;i<n;i++)
		for (j=0;j<n;j++)
		a[i][j] = begin_a[i][j];
}

int all_lu(int u[],int v[],int n)
{
	int k,i,j,p,q;
	double max;
	for (k=0;k<=n-2;k++)
	{
		max = 0.0; p = k; q = k;

		for (i=k;i<=n-1;i++)
			for (j=k;j<=n-1;j++)
			{
				if (fabsl(a[i][j])>fabsl(max)) 
				{
					max = a[i][j];
					p = i;
					q = j;
				}
			};

		u[k] = p;
		v[k] = q;		
		if (p>k)
			for (j=0;j<n;j++)
			{
				max = a[k][j];
			    a[k][j] = a[p][j];
			    a[p][j] = max;
		    };
		if (q>k)
			for (i=0;i<n;i++)
			{
				max = a[i][k];
			    a[i][k] = a[i][q];
			    a[i][q] = max;
		    };		
		if (a[k][k] == 0.0)
		{
			printf("error\n");
			return (1);
		}		

		for (i=k+1;i<n;i++)
		{
			a[i][k] = a[i][k]/a[k][k];
			for(j=k+1;j<n;j++)
				a[i][j]=a[i][j]-a[i][k]*a[k][j];
		};
	};
	return (0);
}

int row_lu(int u[],int n) 

{
	int i,k,j;
	int p;
	double max;

	for(k=0;k<n-1;k++)
	{
		max=0; p = k;
		
		for(i=k;i<n;i++)
			if(fabsl(max)<fabsl(a[i][k]))
			{
				max=a[i][k];
				p=i;
			};

		u[k]=p;
		if(p!=k)
			for(i=0;i<n;i++)
			{
				max=a[k][i];
				a[k][i]=a[p][i];
				a[p][i]=max;
			};		
		if (a[k][k] == 0.0)
		{
			printf("error\n");
			return (1);
		};

		for(i=k+1;i<n;i++)
		{
			a[i][k]=a[i][k]/a[k][k];
			for(j=k+1;j<n;j++)
				a[i][j]=a[i][j]-a[i][k]*a[k][j];
		};
	};
return 0;
}


void Left(double b[],int n) 
{
	int i,j;
	for(i=0;i<n;i++)
		for(j=i+1;j<n;j++)
			b[j]-=b[i]*a[j][i];
}

void Right(double b[],int n) 
{
	int i,j;
	for(i=n-1;i>0;i--)
	{
		b[i]/=a[i][i];
		for(j=0;j<i;j++)
			b[j]-=b[i]*a[j][i];
	}
	b[0]/=a[0][0];
}

void answer_lu_row()
{
	int r,i,j;
	double t;
	int u[20];

	r = row_lu(u,n); 
	if(r==0)
	{
		for (i=0;i<n;i++)
		{
			for (j=0;j<n;j++)
				printf("%10.5e",a[i][j]);
			printf("\n");
		};
		
		for (i=0;i<n-1;i++) 
		{
			if(u[i]!=i)
			{				
				t=b[i];
				b[i]=b[u[i]];
				b[u[i]]=t;
			}
		};
		Left(b,n); 
		Right(b,n);
		for (i=0;i<n;i++)
			printf("x[%d]=%f, ",i+1,b[i]);
		printf("\n");
	}
	printf("\n");
}

void answer_lu_all()
{
	int r,i,j;
	double t;
	int u[20];
	int v[20];

	r = all_lu(u,v,n);
	if(r==0)
	{
		for (i=0;i<n;i++)
		{
			for (j=0;j<n;j++)
				printf("%10.5e",a[i][j]);
			printf("\n");
		}
		
		for (i=0;i<n-1;i++) 
		{
			if(u[i]!=i)
			{				
				t = b[i];
				b[i] = b[u[i]];
				b[u[i]] = t;
			}
		};
		Left(b,n); 
		Right(b,n);
		for (i=0;i<n-1;i++)
		{
			if(v[i]!=i)
			{
				t = b[i];
				b[i] = b[v[i]];
				b[v[i]] = t;
			}
		};
		for (i=0;i<n;i++)
			printf("x[%d]=%f, ",i+1,b[i]);
		printf("\n");
	}
	printf("\n");
}



int house(double b[],double c[])
{
	int i,j,k;

	double h,f,g,h2;

	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			q[i][j] = a[i][j];
	for(i=n-1;i>=1;i--)
	{
		h = 0.0;
		if (i>1)
			for (k=0;k<i;k++) h = h+q[i][k]*q[i][k];
		if (h == 0)
		{
			c[i] = 0.0;
			if (i==1) c[i] = q[i][i-1];
			b[i] = 0.0;
		}
		else
		{
			c[i] = sqrtl(h);
			if (q[i][i-1]>0.0) c[i] = -c[i];
			h = h-q[i][i-1]*c[i];
			q[i][i-1] = q[i][i-1]-c[i];
			f = 0.0;
			for (j=0;j<i;j++)
			{
				q[j][i] = q[i][j]/h;
				g = 0.0;
				for(k=0;k<=j;k++) g = g+q[i][k]*q[j][k];
				if (j+1<=i-1)
					for (k=j+1;k<i;k++) g = g+q[k][j]*q[i][k];
				c[j] = g/h;
				f = f+g*q[j][i];
			}
			h2 = f/(h+h); 
			for (j=0;j<i;j++)
			{
				f = q[i][j];
				g = c[j]-h2*f;
				c[j] = g;
				for(k=0;k<=j;k++) q[j][k] = q[j][k] - f*c[k] - g*q[i][k];
			}
			b[i] = h;
		}
	};
	for(i=0;i<n-1;i++) c[i] = c[i+1];
    c[n-1] = 0.0;
	
	b[0] = 0.0;
	for (i=0;i<=n-1;i++)
	{
		if ((b[i]!=0.0)&&(i-1>=0))
			for (j=0;j<i;j++)
			{
				g = 0.0;
				for (k=0;k<i;k++) g = g+q[i][k]*q[k][j];
				for (k=0;k<i;k++) q[k][j] = q[k][j]-g*q[k][i];
			}
		b[i] = q[i][i]; 
		q[i][i] = 1.0;
		if (i-1>=0)
			for(j=0;j<i;j++)
			{
				q[i][j] = 0;
				q[j][i] = 0;
			};
	};
	
	for (i=0;i<n;i++)
		for (j=0;j<n;j++)
			a[i][j] = 0;
	for (i=0;i<n;i++) {a[i][i] = b[i]; a[i][i+1] = c[i];};
	for (i=1;i<n;i++) a[i][i-1] = c[i-1];
	return 0;

}

int qr(int l,double eps,double b[],double c[])
{
	int i,j,k,m,it;
	double d,f,h,g,p,r,e,s;
	c[n-1] = 0;
	d = 0;
	f = 0;
	for (j=0;j<n;j++)
	{
		it = 0;
		h = eps*(fabsl(b[j])+fabsl(c[j]));
		if (h>d) d = h;
		m = j;
		while ((m<=n-1)&&(fabsl(c[m])>d)) m++;
		if (m!=j)
		{
			do
			{
				if (it==l)
				{
					printf("error\n");
					return -1;
				}
				it++;
				g = b[j];
				p = (b[j+1]-g)/(2.0*c[j]);
				r = sqrtl(p*p+1.0);
				if (p>=0.0) b[j] = c[j]/(p+r);
				else b[j] = c[j]/(p-r);
				h = g-b[j];
				for (i=j+1;i<=n-1;i++) b[i] = b[i]-h;
				f = f+h; p =b[m]; e = 1.0; s =0.0;
				for (i=m-1;i>=j;i--)
				{
					g = e*c[i]; h = e*p;
					if (fabsl(p)>=fabsl(c[i]))
					{
						e = c[i]/p; r =sqrtl(e*e+1.0);
						c[i+1] = s*p*r; s = e/r; e =1.0/r;
					}
					else
					{
						e = p/c[i];
						r = sqrtl(e*e+1.0);
						c[i+1] = s*c[i]*r;
						s = 1.0/r; e = e/r;
					}
					p = e*b[i]-s*g;
					b[i+1] = h+s*(e*g+s*b[i]);
					for (k=0; k<n;k++)
					{				
						h = q[k][i+1]; q[k][i+1] = s*q[k][i]+e*h;
						q[k][i] = e*q[k][i] - s*h;
					}
				}
				c[j] = s*p;
				b[j]=e*p;
			}
			while (fabsl(c[j])>d);
		}
		b[j]+=f;
	};
	for (i=0;i<n;i++)
	{
		k = i;
		p = b[i];
		if (i+1<n)
		{
			j = i+1;
			while ((j<=n-1)&&(b[j]<=p))
			{
				k = j; p = b[j]; j = j+1;
			}
		};
		if (k!=i)
		{
			b[k] = b[i];
			b[i] =p;
			for (j=0;j<n;j++)
			{
				p = q[j][i];
				q[j][i] = q[j][k];
				q[j][k] = p;
			}
		}
	}
	return 1;
}





int _tmain(int argc, _TCHAR* argv[])
{
	
	int i,j;
	double zhu[10],ci[10];
    init_a();
	answer_lu_all();
	init_a();
	answer_lu_row();
	init_a();
	house(zhu,ci);
	for (i=0;i<n;i++)
	{
		for (j=0;j<n;j++) printf("%10.5e ",a[i][j]);
		printf("\n");
	};
	printf("\n");
	qr(60,0.000001,zhu,ci);
	for (i=0;i<n;i++)
	{
		for (j=0;j<n;j++) printf("%10.5e ",q[i][j]);
		printf("\n");
	};
	printf("\n");
	for (i=0;i<n;i++)
		printf("%10.5e ",zhu[i]);
	printf("\n");

	return 0;	
}