d1r5.cpp

vc++常用数值算法集合

#include "iostream.h"
#include "math.h"
#include "stdlib.h"
#include "iomanip.h"
#include "randnumber.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 mprove(double a[], double alud[], int n, int indx[], double b[], double x[])
{
	int i, j;
    const int nmax = 100;
    double r[100];
	double sdp;
    for (i = 1; i<=n; i++)
	{
        sdp = -b[i];
        for (j = 1; j<=n; j++)
		{
            sdp = sdp + a[(i-1)*n+j] * x[j];
        }
        r[i] = sdp;
    }
    lubksb(alud, n, indx, r);
    for (i = 1; i<=n; i++)
	{
        x[i] = x[i] - r[i];
    }
}

void main()
{
    //program d1r5
    //driver program for routine mprove
    int i,j,n = 5;
	int d=0;
    double a[26], a1[6][6], a2[26], b[6], b1[6];
	int indx[6];
	CRandNumber rand;
    //输入已知的方程组的系数矩阵
    a1[1][1]=1; a1[1][2]= 2; a1[1][3] = 3; a1[1][4] = 4; a1[1][5] = 5;
    a1[2][1]=2; a1[2][2]= 3; a1[2][3] = 4; a1[2][4] = 5; a1[2][5] = 1;
    a1[3][1]=1; a1[3][2]= 1; a1[3][3] = 1; a1[3][4] = 1; a1[3][5] = 1;
    a1[4][1]=4; a1[4][2]= 5; a1[4][3] = 1; a1[4][4] = 2; a1[4][5] = 3;
    a1[5][1]=5; a1[5][2]= 1; a1[5][3] = 2; a1[5][4] = 3; a1[5][5] = 4;
    //输入已知的方程组的右端向量b
    b[1] = 1;
    b[2] = 1;
    b[3] = 1;
    b[4] = 1;
    b[5] = 1;
    for (i = 1; i<=n; i++)
	{
        for (j = 1; j<=n;j++)
		{
            a[(i-1)*n+j] = a1[i][j];
			a2[(i-1)*n+j] = a1[i][j];
        }
    }
    for (i = 1; i<=n; i++)
	{
        b1[i] = b[i];
    }
    ludcmp(a, n, indx, d);
    lubksb(a, n, indx, b);
    //输出方程组的解b
    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;
	cout.width(10);    cout<<b[4]<<endl;
	cout.width(10);    cout<<b[5]<<endl;
    for (i = 1; i<=n; i++)
	{
        b[i] = b[i] * (1.0 + 0.2 * rand.fRandom());
    }
    //输出干扰后的解b
    cout<<endl;
    cout<<"输出干扰后的解"<<endl;
	cout.width(12);    cout<<b[1]<<endl;
	cout.width(12);    cout<<b[2]<<endl;
	cout.width(12);    cout<<b[3]<<endl;
	cout.width(12);    cout<<b[4]<<endl;
	cout.width(12);    cout<<b[5]<<endl;
    mprove(a2, a, n, indx, b1, b);
    //输出改善后的解b
    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;
	cout.width(10);    cout<<b[4]<<endl;
	cout.width(10);    cout<<b[5]<<endl;
}