ludecompose.cpp

这是我编写的一些关于数值分析算法

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

#include "stdafx.h"
#include <math.h>
#include <iostream.h>
#define n 4

void Decompose(double a[n][n],double l[n][n],double u[n][n]);
void Gauss(double a[n][n],double b[n]);
void Exchange(double a[n][n],double b[n],int i);

int main(int argc, char* argv[])
{
	double a[n][n]={{1,2,3,4},{1,4,9,16},{1,8,27,64},{1,16,81,256}};
	double b[n]={2,10,44,190};
	double l[n][n]={0};
	double u[n][n]={0};
	Decompose(a,l,u);
	for (int i=0;i<4;i++)
	{
		for (int j=0;j<4;j++)
		{
			cout<<u[i][j]<<"    ";
		}
		cout<<endl;
	}
	Gauss(l,b);
	Gauss(u,b);
	for( i=0;i<n;i++)
		cout<<b[i]<<endl;
	return 0;
}
void Decompose(double a[n][n],double l[n][n],double u[n][n])
{
	double temp;
	for (int i=0;i<n;i++)
	{
		u[0][i]=a[0][i];
		l[i][0]=a[i][0]/u[0][0];
	}
	l[0][0]=1;

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

	}
}

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

void Exchange(double a[n][n],double b[],int i)
{
	int k=i;
	double temp;
	double max=fabs(a[i][i]);
	for (int j=i+1;j<n;j++)
	{
		if(fabs(a[j][i])>max)
		{
			max=fabs(a[j][i]);
			k=j;
		}
	}
	temp=b[i];
	b[i]=b[k];
	b[k]=temp;
	for (j=i;j<n;j++)
	{
		temp=a[i][j];
		a[i][j]=a[k][j];
		a[k][j]=temp;
	}
}