d1r2.cpp

数值计算c++源代码,包括各种算法。很有用的。

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

void ludcmp(double a[], int n, int indx[], int d)
{
	const int nmax=100;
	const double tiny=1e-20;
	double vv[100];
	double sum,dum,aamax;
	int i,j,k,imax;
	d=1;
	for (i=1; i<=n; i++)
	{
		aamax=0.0;
		for (j=1; j<=n; j++)
		{
			if (fabs(a[(i-1)*n+j])>aamax)
			{
				aamax=fabs(a[(i-1)*n+j]);
			}
		}
		if (aamax==0.0)
		{
			cout<<"singular matrix"<<endl;
		}
		vv[i]=1.0/aamax;
	}
	for (j=1; j<=n; j++)
	{
        if (j > 1)
		{
            for (i=1; i<=j-1; i++)
			{
                sum=a[(i-1)*n+j];
                if (i>1)
				{
                    for (int k=1; k<=i-1; k++)
					{
                        sum=sum-a[(i-1)*n+k]*a[(k-1)*n+j];
                    }
                    a[(i-1)*n+j]=sum;
                }
            }
        }
        aamax=0.0;
        for (i=j; i<=n; i++)
		{
            sum=a[(i-1)*n+j];
            if (j>1)
			{
                for (k=1; k<=j-1; k++)
				{
                    sum=sum-a[(i-1)*n+k]*a[(k-1)*n+j];
				}
                a[(i-1)*n+j]=sum;
            }
            dum=vv[i]*fabs(sum);
            if (dum>=aamax)
			{
                imax=i;
                aamax=dum;
            }
        }
        if (j!=imax)
		{
            for (k=1; k<=n; k++)
			{
                dum=a[(imax-1)*n+k];
                a[(imax-1)*n+k]=a[(j-1)*n+k];
                a[(j-1)*n+k]=dum;
            }
            d=-d;
            vv[imax]=vv[j];
        }
        indx[j]=imax;
        if (j!=n)
		{
            if (a[(j-1)*n+j]==0.0)
			{
				a[(j-1)*n+j]=tiny;
			}
            dum=1.0/a[(j-1)*n+j];
            for (i=j+1; i<=n; i++)
			{
                a[(i-1)*n+j]=a[(i-1)*n+j]*dum;
            }
		}
    }
    if (a[(n-1)*n+n]==0.0)
	{
		a[(n-1)*n+n]=tiny;
	}
}

void lubksb(double a[], int n, int indx[], double b[])
{
	int i,j,ll;
	double sum;
    int ii = 0;
    for (i = 1; i<=n; i++)
	{
        ll = indx[i];
        sum = b[ll];
        b[ll] = b[i];
        if (ii != 0 )
		{
            for (j = ii; j<=i-1; j++)
			{
                sum = sum - a[(i-1)*n+j] * b[j];
            }
		}
        else
		{
			if (sum != 0.0)
			{
				ii = i;
			}
        }
        b[i] = sum;
    }
    for (i = n; i>=1; i--)
	{
        sum = b[i];
        if (i < n)
		{
            for (j = i + 1; j<=n; j++)
			{
                sum = sum - a[(i-1)*n+j] * b[j];
            }
        }
        b[i] = sum / a[(i-1)*n+i];
    }
}

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;
}