solveequal.cpp

数值分析过程中

#include"solveEqual.h"
#include<iostream>
#define N 2
using namespace std;
//解下三角矩阵
void solveDowntriEqual(double L[N][N],double d[N], double x[N]);
//解上三角矩阵
void solveUptriEqual(double U[N][N],double d[N], double x[N]);
//解线性方程
void solveEqual( double A[N][N],double d[N], double x[N]) 
{
double L[N][N],U[N][N];
double y[N];
//cout<<"before triangle"<<endl;
triangle(A,L,U);
//cout<<"after trangle:"<<endl;
//print(A);
//print(L);
//print(U);
solveDowntriEqual(L,d,y);
solveUptriEqual(U,y,x);
}

//矩阵的三角分解
void triangle(double A[N][N],double L[N][N],double U[N][N])
{
	for(int i = 0; i <= N-1; i++) {
		U[0][i] = A[0][i];
		L[i][0] = A[i][0]/U[0][0];
	}
	/*
	cout<<"U:"<<endl;
	print(U[0]);
	cout<<"L:"<<endl;
	for(int i=0; i<= N-1; i++)
	cout<<"\t"<<L[i][0];
	cout<<endl;
	*/
	//先求出u的r-1行和L的r-1列, 在计算U的第r行,L的第r列
	for(int r = 1; r <= N-1; r++)
	{
		for(int i = 0 ;i <= N-1; i++)
		{
			U[r][i] = A[r][i];
			for(int k = 0; k <= r-1 ; k++)
			{
				U[r][i] = U[r][i] - L[r][k]*U[k][i];
			}
		}
		for(int i = 0; i <= N-1; i++)
		{
			L[i][r] = A[i][r];
			for(int k = 0; k <= r-1; k++)
			{
				L[i][r] -= L[i][k]*U[k][r];
			}
			L[i][r] = L[i][r]/U[r][r];
		}
	}

}

//解下三角矩阵
void solveDowntriEqual(double L[N][N],double d[N],double x[N])
{
x[0] = d[0];//x1 = d1
for(int i = 1;i <= N-1;i++)
{
	x[i] = d[i];
	for(int j = 0; j <= i-1 ;j++)
	{
		x[i] = x[i] - L[i][j]*x[j];	
	}
}
}
//解上三角矩阵
void solveUptriEqual(double U[N][N],double d[N],double x[N])
{
x[N-1] = d[N-1]/U[N-1][N-1];
for(int i = N-2; i >= 0; i--)
{
	x[i] = d[i];
	for(int j = i+1; j <= N-1;j++)
	{
		x[i] = x[i] - U[i][j]*x[j];
	}
	x[i] = x[i]/U[i][i];
}
}
//输出格式
void print(double X[N])
{
	for(int i = 0 ;i <= N-1; i++)
	{
	cout<<"\t"<<X[i];
	}
	cout<<endl;
}

void print(double X[N][N])
{
	for(int i = 0 ;i <= N-1; i++)print(X[i]);
}