d1r11.cpp

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

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

int sgn(double pa)
{
	if (pa>0.0)
	{
		return 1;
	}
	else
	{
		if (pa<0.0)
		{
			return -1;
		}
	}
	return 0;
}

void qrdcmp(double a[], int m, int n, double q[])
{
	int i,j,k;
	double s,t,h,f;
    for (i = 1; i<=m; i++)
	{
        for (j = 1; j<=m; j++)
		{
            q[(i-1)*n+j] = 0.0;
        }
        q[(i-1)*n+i] = 1.0;
    }
    for (k = 1; k<=m-1; k++)
	{
        s = 0.0;
        for (i = k; i<=m;i++)
		{
            s = s + fabs(a[(i-1)*n+k]);
        }
        if (s != 0.0)
		{
            t = 0.0;
            for (i = k; i<=m; i++)
			{
                a[(i-1)*n+k] = a[(i-1)*n+k] / s;
                t = t + a[(i-1)*n+k] * a[(i-1)*n+k];
            }
            t = -sqrt(t) * sgn(a[(k-1)*n+k]);
            a[(k-1)*n+k] = a[(k-1)*n+k] - t;
            h = -t * a[(k-1)*n+k];
            for (j = k + 1; j<=n; j++)
			{
                f = 0.0;
                for (i = k; i<=m; i++)
				{
                    f = f + a[(i-1)*n+k] * a[(i-1)*n+j];
                }
                f = f / h;
                for (i = k; i<=m; i++)
				{
                    a[(i-1)*n+j] = a[(i-1)*n+j] - a[(i-1)*n+k] * f;
                }
            }
            for (j = 1; j<=m;j++)
			{
                f = 0.0;
                for (i = k; i<=m; i++)
				{
                    f = f + a[(i-1)*n+k] * q[(i-1)*n+j];
                }
                f = f / h;
                for (i = k; i<=m; i++)
				{
                    q[(i-1)*n+j] = q[(i-1)*n+j] - a[(i-1)*n+k] * f;
                }
            }
            a[(k-1)*n+k] = t * s;
            for (i = k + 1; i<=m; i++)
			{
                a[(i-1)*n+k] = 0.0;
            }
        }
    }
}

void qrbksb(double a[], int n, double q[], double b[], double x[])
{
	int i,j;
	double sum;
    for (i = 1; i<=n; i++)
	{
        sum = 0.0;
        for (j = 1; j<=n; j++)
		{
            sum = sum + q[(i-1)*n+j] * b[j];
        }
        x[i] = sum;
    }
    for (i = n; i>=1; i--)
	{
        sum = x[i];
        for (j = i + 1; j<=n; j++)
		{
            sum = sum - a[(i-1)*n+j] * x[j];
        }
        if (a[(i-1)*n+i] == 0.0)
		{
			cout<<"a is singular matrix."<<endl;
		}
        x[i] = sum / a[(i-1)*n+i];
    }
}

void main()
{
    //program d1r11
    //driver program for routine qrbksb,qrdcmp
    int i,j,l,n = 5;
    double a[6][6], b[6], q[26], c[26], r[6], x[6];
    //输入已知的方程组的系数矩阵
    a[1][1]=1.4; a[1][2]=2.1; a[1][3]=2.1; a[1][4]=7.4; a[1][5]=9.6;
    a[2][1]=1.6; a[2][2]=1.5; a[2][3]=1.1; a[2][4]=0.7; a[2][5]=5.0;
    a[3][1]=3.8; a[3][2]=8.0; a[3][3]=9.6; a[3][4]=5.4; a[3][5]=8.8;
    a[4][1]=4.6; a[4][2]=8.2; a[4][3]=8.4; a[4][4]=0.4; a[4][5]=8.0;
    a[5][1]=2.6; a[5][2]=2.9; a[5][3]=0.1; a[5][4]=9.9; a[5][5]=7.7;
    //输入已知的方程组的右端向量b
    b[1] = 1.1;
    b[2] = 1.6;
    b[3] = 4.7;
    b[4] = 9.1;
    b[5] = 0.1;
    cout<<endl;
    cout<<"已知的方程组的右端向量"<<endl;
	cout.setf(ios::fixed|ios::right);
	cout.precision(5);
	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;
    for (i = 1; i<=n; i++)
	{
        for (j = 1; j<=n; j++)
		{
           c[(i-1)*n+j] = a[i][j];
        }
    }
    qrdcmp(c, n, n, q);
    for (i = 1; i<=n; i++)
	{
        r[i] = b[i];
    }
    qrbksb(c, n, q, r, x);
	cout<<endl;
    cout<<"计算出的方程组的解"<<endl;
	cout.width(12);    cout<<x[1]<<endl;
	cout.width(12);    cout<<x[2]<<endl;
	cout.width(12);    cout<<x[3]<<endl;
	cout.width(12);    cout<<x[4]<<endl;
	cout.width(12);    cout<<x[5]<<endl;
    //将计算出的解乘以系数矩阵，以验证计算结果正确
    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(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;
}