lu.cpp

LU分解法解线性方程组

#include "iostream.h"
#include "math.h"

void main()
{
    //program d1r2
    //driver program for routine lubksb,ludcmp
	int i,j,l,p=0;
    int n = 3;
    double a[4][4], b[4], a1[10], x[4];
	int indx[4];
    //输入已知的方程组的系数矩阵
    a[1][1]=1; a[1][2]=2; a[1][3]=3;
    a[2][1]=2; a[2][2]=2; a[2][3]=3;
    a[3][1]=3; a[3][2]=3; a[3][3]=3;
    //输入已知的方程组的右端向量
    b[1]=1;
    b[2]=2;
    b[3]=3;
    cout<<endl;
    cout<<"已知的方程组的右端向量"<<endl;
	cout.width(10);
    cout<<b[1]<<endl;
	cout.width(10);
	cout<<b[2]<<endl;
	cout.width(10);
	cout<<b[3]<<endl;
    for (i = 1; i<=n; i++)
	{
        for (j = 1; j<=n; j++)
		{
            a1[(i-1)*n+j] = a[i][j];
        }
    }
    ludcmp(a1, n, indx, p);
    for (l = 1; l<=n; l++)
	{
        x[l] = b[l];
    }
    lubksb(a1, n, indx, x);
    cout<<endl;
    cout<<"计算出的方程组的解"<<endl;
	cout.width(10);
    cout<<x[1]<<endl;
	cout.width(10);
	cout<<x[2]<<endl;
	cout.width(10);
	cout<<x[3]<<endl;
    //将计算出的b乘以系数矩阵，以验证计算结果正确
    for (l = 1; l<=n; l++)
	{
        b[l] = 0.0;
        for (j = 1; j<=n; j++)
		{
            b[l] = b[l] + a[l][j] * x[j];
        }
    }
    cout<<endl;
    cout<<"计算出的解乘以系数矩阵的结果"<<endl;
	cout.width(10);
    cout<<b[1]<<endl;
	cout.width(10);
	cout<<b[2]<<endl;
	cout.width(10);
	cout<<b[3]<<endl;
}