sfroid.cpp

此代码为VC++常用数值算法这本书中附赠的光D部分的源代码

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

double x[42], h, c2, anorm;
int mm, n;
void main()
{
    //program sfroid
    int ne = 3; int m = 41; int nb = 1; int nci = 3;
	int ncj = 3; int nck = 42; int nsi = 3; int nsj = 7; 
	int nyj = 3; int nyk = 41;
    double scalv[4]; double y[124], c[379], s[22];
	int indexv[4]; int itmax = 100;
    double conv = 0.000005; double slowc = 1.0;

    h = 1.0 / (m - 1);
    c2 = 0.0;
    cout<<"enter m,n"<<endl;
    mm = 2;
    n = 2;
    if ((n + mm % 2) == 1)
	{
        indexv[1] = 1;
        indexv[2] = 2;
        indexv[3] = 3;
	}
    else
	{
        indexv[1] = 2;
        indexv[2] = 1;
        indexv[3] = 3;
    }
    anorm = 1.0;
	int i,k;
	double q1;
    if (mm != 0)
	{
        q1 = (double)n;
        for (i = 1; i<=mm; i++)
		{
            anorm = -0.5 * anorm * (n + i)*(q1 / (double)i);
            q1 = q1 - 1.0;
        }
    }
	double fac1,fac2,deriv;
    for (k = 1; k<=m - 1; k++)
	{
        x[k] = (k - 1) * h;
        fac1 = 1.0 - x[k] * x[k];
        fac2 = pow(fac1 , (-(double)mm / 2.0));
        y[k] = plgndr(n, mm, x[k]) * fac2;
        deriv = -((n - mm + 1) * plgndr(n + 1, mm, x[k]) -
			    (n + 1) * x[k] * plgndr(n, mm, x[k])) / fac1;
        y[nyk+k] = mm * x[k] * y[k] / fac1 + deriv * fac2;
        y[2*nyk+k] = double(n * (n + 1) - mm * (mm + 1));
    }
    x[m] = 1.0;
    y[m] = anorm;
    y[2*nyk+m] = n * (n + 1) - mm * (mm + 1);
    y[nyk+m] = (y[2*nyk+m] - c2) * y[m] / (2.0 * (mm + 1.0));
    scalv[1] = fabs(anorm);
    if (y[nyk+m] > fabs(anorm))
	{
        scalv[2] = y[nyk+m];
	}
    else
	{
        scalv[2] = fabs(anorm);
    }
    if (y[2*nyk+m] > 1)
	{
        scalv[3] = y[2*nyk+m];
	}
    else
	{
        scalv[3] = 1;
    }
	cout.setf(ios::fixed|ios::right);
	cout.precision(6);
    do
	{
		cout<<"enter c^2 or  999 to end"<<endl;
		cin>>c2;
		if (c2 == 999)
		{
			break;
		}
		solvde(itmax, conv, slowc, scalv, indexv, ne, nb, m, y,
			   nyj, nyk, c, nci, ncj, nck, s, nsi, nsj);
		cout<<"m = ";
		cout.width(1);		cout<<mm;
		cout<<"   n = ";
		cout.width(1);		cout<<n;
		cout<<"   c^2 = ";
		cout.width(6);		cout<<c2;
		cout<<"   lambda = ";
		cout.width(8);
		cout<<y[2*nyk+1] + mm * (mm + 1)<<endl;
    }while(1);
}