d15r1.cpp

VC++常用数据算法集,里面包含了大量的常用算法程序,很有用的哟!

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

void derivs(double x, double y[], double dydx[])
{
    dydx[1] = y[2];
    dydx[3] = 0.0;
    dydx[2] = (2.0 * x * (m + 1.0) * y[2]-(y[3] - c2 * x * x)*
		       y[1]) / (1.0 - x * x);
    dydx[4] = y[5];
    dydx[6] = 0.0;
    dydx[5] = (2.0 * x * (m + 1.0) * y[5]-(y[6] + c2 * x * x)*
		       y[4]) / (1.0 - x * x);
}

double c2, factr;
int m, n;
void main()
{
    //program d15r1
    //driver for routine shoot
    //Solves for eigenvalues of spheroidal harmonics. Both
    //Prolate and Oblate case are handled simultaneously,
    //leading to six first-order equations. Unknown to
    //SHOOT, these are actually two independent sets of
    //three coupled equations, one set with c^2 positive
    //and the other with c^2 negative.
    int nvar = 6; int n2 = 2;
	double delta = 0.001; double eps = 0.000001;
    double v[3], delv[3], f[3], dv[3];
    double dx = 0.0001;
    do
	{
		cout<<"input m,n,c-squared (999 to end)"<<endl;
		m = 2;
		n = 2;
		c2 = 0.1;
		if (c2 == 999)
		{
			exit(1);
		}
	}while ((n < m) || (m < 0) || (n < 0));
	cout.setf(ios::fixed|ios::right);
	cout.precision(6);
	cout.width(1);
    cout<<m;
	cout.width(6);
	cout<<n;
	cout.width(12);
	cout<<c2<<endl;
    factr = 1;
	int q1;
    if (m != 0)
	{
        q1 = n;
        for (int i = 1; i<=m; i++)
		{
            factr = -0.5*factr*(n + i)*((double)q1/(double)i);
            q1 = q1 - 1;
        }
    }
    v[1] = n * (n + 1) - m * (m + 1) + c2 / 2.0;
    v[2] = n * (n + 1) - m * (m + 1) - c2 / 2.0;
    delv[1] = delta * v[1];
    delv[2] = delv[1];
    double h1 = 0.1;
    double hmin = 0.0;
    double x1 = -1.0 + dx;
    double x2 = 0.0;
	cout<<endl;
    cout<<"            prolate                   oblate"<<endl;
    cout<<"  mu(m,n)   error est.   mu(m,n)    error est."<<endl;
    do
	{
		shoot(nvar, v, delv, n2, x1, x2, eps, h1, hmin, f, dv);
		cout.width(13);
		cout<<v[1];
		cout.width(13);
		cout<<dv[1];
		cout.width(13);	
		cout<<v[2];
		cout.width(13);
		cout<<dv[2]<<endl;
    }while (fabs(dv[1]) > fabs(eps * v[1]) || fabs(dv[2]) >
		    fabs(eps * v[2]));
}