d9r7.cpp

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

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

void fleg(double x, double pl[], int nl)
{
	pl[1] = 1.0;
	pl[2] = x;
	if (nl > 2)
	{
		double twox = 2.0 * x;
		double f2 = x;
		double d = 1.0;
		for (int j = 3; j<=nl; j++)
		{
			double f1 = d;
			f2 = f2 + twox;
			d = d + 1.0;
			pl[j] = (f2 * pl[j - 1] - f1 * pl[j - 2]) / d;
		}
	}
}

double plgndr(int l, int m, double x)
{
	int i;
	double pll,pmmp1,somx2,fact;
    if (m < 0 || m > l || fabs(x) > 1.0)
	{
		cout<<"bad arguments int plgndr"<<endl;
		exit(1);
	}
    double pmm = 1.0;
    if (m > 0)
	{
        somx2 = sqrt((1.0 - x) * (1.0 + x));
        fact = 1.0;
        for (i = 1; i<=m; i++)
		{
            pmm = -pmm * fact * somx2;
            fact = fact + 2.0;
        }
    }
    if (l == m)
	{
        return pmm;
	}
    else
	{
        pmmp1 = x * (2 * m + 1) * pmm;
        if (l == m + 1)
		{
            return pmmp1;
		}
        else
		{
            for (int ll = m + 2; ll<=l; ll++)
			{
                pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
                pmm = pmmp1;
                pmmp1 = pll;
            }
            return pll;
        }
    }
}

void main()
{
    //program d9r7
    //driver for routine fleg
    int nval = 5;
    double x,dx = 0.2;
    int i,j,npoly = 5;
    double afunc[6];
    cout<<"Legendre polynomials"<<endl;
    cout<<"      n=1        n=2        n=3        n=4        n=5"<<endl;
	cout.setf(ios::fixed);
	cout.precision(4);
    for (i = 1; i<=nval; i++)
	{
        x = i * dx;
        fleg(x, afunc, npoly);
        cout<<" x = ";
		cout<<x<<endl;
        for (j = 1; j<=npoly; j++)
		{
			cout.width(11);
            cout<<afunc[j];
        }
        cout<<"   routine fleg"<<endl;
        for (j = 1; j<=npoly; j++)
		{
			cout.width(11);
            cout<<plgndr(j - 1, 0, x);
        }
        cout<<"   routine plgndr";
        cout<<endl;
    }
}