d15r2.cpp

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

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

double c2, factr,dx;
int m, n;
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);
}

void load1(double x1, double v1[], double y[])
{
    y[3] = v1[1];
    y[2] = -(y[3] - c2) * factr / 2.0 / (m + 1.0);
    y[1] = factr + y[2] * dx;
}

void load2(double x2, double v2[], double y[])
{
    y[3] = v2[2];
    y[2] = (y[3] - c2) * y[1] / 2.0 / (m + 1.0);
    y[1] = v2[1];
}

void score(double xf, double y[], double f[])
{
    for (int i = 1; i<=3; i++)
	{
        f[i] = y[i];
    }
}

void main()
{
    //program d15r2
    //driver for routine shootf
    int nvar = 3; int n1 = 2; int n2 = 1;
	double delta = 0.001; double eps = 0.000001;
    double dxx = 0.0001; int kmax = 100;
    double v1[2],delv1[2],v2[3],delv2[3],dv1[2],dv2[3],f[4];
    do
	{
		cout<<"input m,n,c-squared (999 to end)"<<endl;
		m = 2;
		n = 2;
		c2 = 0.1;
		if (c2 == 999)
		{
			exit(1);
		}
		if ((n < m) || (m < 0))
		{
			cout<<"improper arguments"<<endl;
		}
    }while(n < m || m < 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 i,q1;
    if (m != 0)
	{
        q1 = n;
        for (i = 1; i<=m; i++)
		{
            factr = -0.5*factr*(n + i)*((double)q1/(double)i);
            q1 = q1 - 1;
        }
    }
    dx = dxx;
    v1[1] = n * (n + 1) - m * (m + 1) + c2 / 2.0;
    if (((n - m) % 2) == 0)
	{
        v2[1] = factr;
	}
    else
	{
        v2[1] = -factr;
    }
    v2[2] = v1[1] + 1;
    delv1[1] = delta * v1[1];
    delv2[1] = delta * factr;
    delv2[2] = delv1[1];
    double h1 = 0.1;
    double hmin = 0.0;
    double x1 = -1 + dx;
    double x2 = 1 - dx;
    double xf = 0.0;
    cout<<"        mu(-1)      y(1-dx)     mu(+1)"<<endl;
    do
	{
		shootf(nvar, v1, v2, delv1, delv2, n1, n2, x1, x2,
			   xf, eps, h1, hmin, f, dv1, dv2);
		cout<<"  v   ";
		cout.width(12);		cout<<v1[1];
		cout.width(12);		cout<<v2[1];
		cout.width(12);		cout<<v2[2]<<endl;
		cout<<"  dv  ";
		cout.width(12);		cout<<dv1[1];
		cout.width(12);		cout<<dv2[1];
		cout.width(12);		cout<<dv2[2]<<endl;
		cout<<endl;
    }while (fabs(dv1[1]) > fabs(eps * v1[1]));
}