ludecomp.cpp

此代码可有于计算方法中的各种数学计算,如解线性方程组,求积分等,很方便!

#define N  3       //方程组阶数
#define EPS 1e-16  //精度要求

#include "math.h"
#include <iostream>
using namespace std;

int main()
{
	int i,j;
	double a[N+1][N+1],b[N+1];
	void LUdecomp(double a[][N+1]);
	void Solve(double a[][N+1],double b[]);
	
	//输入系数矩阵
	cout<<"请输入系数矩阵:"<<endl;
	for(i=1;i<=N;i++)
		for(j=1;j<=N;j++)
			cin>>a[i][j];
	
	//输入右边的值
	cout<<"请输入右边的值:"<<endl;
	for(i=1;i<=N;i++) 
		cin>>b[i];
	
	LUdecomp(a);      //调用函数对方程组进行LU分解
	
	//打印LU分解后的结果
	cout<<"LU分解后的结果为:"<<endl;
	for(i=1;i<=N;i++)
	{
		for(j=1;j<=N;j++)
			cout<<" "<<a[i][j]<<" ";
		cout<<endl;
	}
	
	//对方程组求解
	cout<<"方程组的解为:"<<endl;
	if(fabs(a[N][N]>EPS))
	{
		Solve(a,b);
		cout<<endl;
		for(i=1;i<=N;i++)
			cout<<"x"<<i<<"="<<b[i]<<"  ";
	}
	return 0;
}

//对系数矩阵进行LU分解
void LUdecomp(double a[][N+1])
{
	int i,j,k,s;
	for(k=1;k<=N;k++)
	{
		for(j=k;j<=N;j++)
			for(s=1;s<=k-1;s++)
				a[k][j]-=a[k][s]*a[s][j];
			if((fabs(a[k][k])<EPS)&&(k<N))
				cout<<"LU分解失败!"<<endl;
				
			for(i=k+1;i<=N;i++)
			{
				for(s=1;s<=k-1;s++)
					a[i][k]-=a[i][s]*a[s][k];
				a[i][k]/=a[k][k];
			}
	}
}

//求解方程组
void Solve(double a[][N+1],double b[])
{
	int k,j;
	for(k=1;k<=N;k++)
		for(j=1;j<=k-1;j++)
			b[k]-=a[k][j]*b[j];
		for(k=N;k>=1;k--)
		{
			for(j=k+1;j<=N;j++)
				b[k]-=a[k][j]*b[j];
			b[k]/=a[k][k];
		}
}