d1r3.cpp

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

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

void tridag(double a[],double b[], double c[],double r[], double u[], int n)
{ 
    const int nmax = 100;
    double gam[100];
    if (b[1] == 0.0)
	{
		return;
	}
    double bet = b[1];
    u[1] = r[1] / bet;
    for (int j = 2; j<=n; j++)
	{
        gam[j] = c[j - 1] / bet;
        bet = b[j] - a[j] * gam[j];
        if (bet == 0.0)
		{
			return;
		}
        u[j] = (r[j] - a[j] * u[j - 1]) / bet;
    }
    for (j = n - 1; j>=1; j--)
	{
        u[j] = u[j] - gam[j + 1] * u[j + 1];
    }
}

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