d15r2.cpp

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

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

//int kmax, kount, dxsav;
double c2, factr,dx;
int m, n;

int sgn(double f)
{
	if (f>0)
	{
		return 1;
	}
	if (f<0)
	{
		return -1;
	}
	return 0;
}

void lubksb(double a[], int n, int an,int indx[], double b[])
{
	int i,j,ll;
	double sum;
    int ii = 0;
    for (i = 1; i<=n; i++)
	{
        ll = indx[i];
        sum = b[ll];
        b[ll] = b[i];
        if (ii != 0 )
		{
            for (j = ii; j<=i-1; j++)
			{
                sum = sum - a[(i-1)*an+j] * b[j];
            }
		}
        else
		{
			if (sum != 0.0)
			{
				ii = i;
			}
        }
        b[i] = sum;
    }
    for (i = n; i>=1; i--)
	{
        sum = b[i];
        if (i < n)
		{
            for (j = i + 1; j<=n; j++)
			{
                sum = sum - a[(i-1)*an+j] * b[j];
            }
        }
        b[i] = sum / a[(i-1)*an+i];
    }
}

void ludcmp(double a[], int n, int an,int indx[], int& d)
{
	const int nmax=100;
	const double tiny=1e-20;
	double vv[100];
	double sum,dum,aamax;
	int i,j,k,imax;
	d=1;
	for (i=1; i<=n; i++)
	{
		aamax=0.0;
		for (j=1; j<=n;j++)
		{
			if (fabs(a[(i-1)*an+j])>aamax)
			{
				aamax=fabs(a[(i-1)*an+j]);
			}
		}
		if (aamax==0.0)
		{
			cout<<"singular matrix"<<endl;
		}
		vv[i]=1.0/aamax;
	}
	for (j=1; j<=n;j++)
	{
        if (j > 1)
		{
            for (i=1; i<=j-1;i++)
			{
                sum=a[(i-1)*an+j];
                if (i>1)
				{
                    for (int k=1; k<=i-1; k++)
					{
                        sum=sum-a[(i-1)*an+k]*a[(k-1)*an+j];
                    }
                    a[(i-1)*an+j]=sum;
                }
            }
        }
        aamax=0.0;
        for (i=j; i<=n; i++)
		{
            sum=a[(i-1)*an+j];
            if (j>1)
			{
                for (k=1; k<=j-1; k++)
				{
                    sum=sum-a[(i-1)*an+k]*a[(k-1)*an+j];
				}
                a[(i-1)*an+j]=sum;
            }
            dum=vv[i]*fabs(sum);
            if (dum>=aamax)
			{
                imax=i;
                aamax=dum;
            }
        }
        if (j!=imax)
		{
            for (k=1; k<=n; k++)
			{
                dum=a[(imax-1)*an+k];
                a[(imax-1)*an+k]=a[(j-1)*an+k];
                a[(j-1)*an+k]=dum;
            }
            d=-d;
            vv[imax]=vv[j];
        }
        indx[j]=imax;
        if (j!=n)
		{
            if (a[(j-1)*n+j]==0.0)
			{
				a[(j-1)*n+j]=tiny;
			}
            dum=1.0/a[(j-1)*an+j];
            for (i=j+1; i<=n; i++)
			{
                a[(i-1)*an+j]=a[(i-1)*an+j]*dum;
            }
		}
    }
    if (a[(n-1)*an+n]==0.0)
	{
		a[(n-1)*an+n]=tiny;
	}
}

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 rk4(double y[], double dydx[], int n, double x, double h, double yout[])
{
	double yt[11], dyt[11], dym[11];
	double hh = h * 0.5;
	double h6 = h / 6.0;
	double xh = x + hh;
	int i;
	for (i = 1; i<=n; i++)
	{
		yt[i] = y[i] + hh * dydx[i];
	}
	derivs(xh, yt, dyt);
	for (i = 1; i<=n; i++)
	{
		yt[i] = y[i] + hh * dyt[i];
	}
	derivs(xh, yt, dym);
	for (i = 1; i<=n; i++)
	{
		yt[i] = y[i] + h * dym[i];
		dym[i] = dyt[i] + dym[i];
	}
	derivs(x + h, yt, dyt);
	for (i = 1; i<=n; i++)
	{
		yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + 2.0 * dym[i]);
	}
	//erase dym, dyt, yt
}

void rkqc(double y[], double dydx[], int n, double& x, double htry, double eps,
		  double yscal[], double& hdid, double& hnext)
{
	double one = 1.0;
	double safety = 0.9;
	double errcon = 0.0006;
	double fcor = 0.066667;
	double ytemp[11], ysav[11], dysav[11];
	double pgrow = -0.2;
	double pshrnk = -0.25;
	double xsav = x;
	int i,flag;
	for (i = 1; i<=n; i++)
	{
		ysav[i] = y[i];
		dysav[i] = dydx[i];
	}
	double hh,h = htry;
	do
	{
		hh = 0.5 * h;
		rk4(ysav, dysav, n, xsav, hh, ytemp);
		x = xsav + hh;
		derivs(x, ytemp, dydx);
		rk4(ytemp, dydx, n, x, hh, y);
		x = xsav + h;
		if (x == xsav)
		{
			cout<<" stepsize not significant in rkqc."<<endl;
			return;
		}
		rk4(ysav, dysav, n, xsav, h, ytemp);
		double errmax = 0.0;
		for (i = 1; i<=n; i++)
		{
			ytemp[i] = y[i] - ytemp[i];
			if (errmax < fabs(ytemp[i] / yscal[i]))
			{
				errmax = fabs(ytemp[i] / yscal[i]);
			}
		}
		errmax = errmax / eps;
		if (errmax > one)
		{
			h = safety * h * pow(errmax , pshrnk);
			flag = 1;
		}
		else
		{
			hdid = h;
			if (errmax > errcon)
			{
				hnext = safety * h * pow(errmax , pgrow);
			}
			else
			{
				hnext = 4.0 * h;
			}
			flag = 0;
		}
	}while(flag == 1);
	for (i = 1; i<=n; i++)
	{
		y[i] = y[i] + ytemp[i] * fcor;
	}
	//erase dysav, ysav, ytemp
}

void odeint(double ystart[], int nvar, double x1, double x2, double eps,
			double h1, double hmin, int& nok, int& nbad, double xp[],
			double yp[], int ypn)
{
	int kmax=0;
    int maxstp = 10000;
    double two = 2.0;
    double zero = 0.0;
    double tiny = 1e-30;
	double xsav=0;
	double dxsav=0;
    double yscal[11], y[11], dydx[11];
    double x = x1;
    double h = fabs(h1) * sgn(x2 - x1);
    nok = 0;
    nbad = 0;
    int i,kount = 0;
    for (i = 1; i<=nvar; i++)
	{
        y[i] = ystart[i];
    }
    if (kmax > 0)
	{
		xsav = x - dxsav * two;
	}
	int nstp;
    for (nstp = 1; nstp<=maxstp; nstp++)
	{
        derivs(x, y, dydx);
        for (i = 1; i<=nvar; i++)
		{
            yscal[i] = fabs(y[i]) + fabs(h * dydx[i]) + tiny;
        }
        if (kmax > 0)
		{
            if (fabs(x - xsav) > fabs(dxsav))
			{
                if (kount < kmax - 1)
				{
                    kount = kount + 1;
                    xp[kount] = x;
                    for (i = 1; i<=nvar; i++)
					{
                        yp[(i-1)*ypn+kount] = y[i];
                    }
                    xsav = x;
                }
            }
        }
        if ((x + h - x2) * (x + h - x1) > zero)
		{
			h = x2 - x;
		}
		double hnext,hdid;
        rkqc(y, dydx, nvar, x, h, eps, yscal, hdid, hnext);
        if (hdid == h)
		{
            nok = nok + 1;
		}
        else
		{
            nbad = nbad + 1;
        }
        if ((x - x2) * (x2 - x1) >= zero)
		{
            for (i = 1; i<=nvar; i++)
			{
                ystart[i] = y[i];
            }
            if (kmax != 0)
			{
                kount = kount + 1;
                xp[kount] = x;
                for (i = 1; i<=nvar; i++)
				{
                    yp[(i-1)*ypn+kount] = y[i];
                }
            }
            //erase dydx, y, yscal
            return;
        }
        if (fabs(hnext) < hmin)
		{
            cout<<" stepsize smaller than minimum."<<endl;
            return;
        }
        h = hnext;
    }
    cout<<" too many steps."<<endl;
}

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 shootf(int nvar, double v1[], double v2[], double delv1[], double delv2[],
			int n1, int n2, double x1, double x2, double xf, double eps,
			double h1, double hmin, double f[], double dv1[], double dv2[])
{
    double y[21], f1[21], f2[21], dfdv[401], xp[201], yp[2001];
	int indx[21];
	int nok,nbad,ypn=200;
    load1(x1, v1, y);
    odeint(y, nvar, x1, xf, eps, h1, hmin, nok, nbad, xp, yp,ypn);
    score(xf, y, f1);
    load2(x2, v2, y);
    odeint(y, nvar, x2, xf, eps, h1, hmin, nok, nbad, xp, yp,ypn);
    score(xf, y, f2);
    int iv,i,j = 0;
	double sav;
    for (iv = 1; iv<=n2; iv++)
	{
        j = j + 1;
        sav = v1[iv];
        v1[iv] = v1[iv] + delv1[iv];
        load1(x1, v1, y);
        odeint(y, nvar, x1, xf, eps, h1, hmin, nok, nbad, xp, yp,ypn);
        score(xf, y, f);
        for (i = 1; i<=nvar; i++)
		{
            dfdv[(i-1)*20+j] = (f[i] - f1[i]) / delv1[iv];
        }
        v1[iv] = sav;
    }
    for (iv = 1; iv<=n1; iv++)
	{
        j = j + 1;
        sav = v2[iv];
        v2[iv] = v2[iv] + delv2[iv];
        load2(x2, v2, y);
        odeint(y, nvar, x2, xf, eps, h1, hmin, nok, nbad, xp, yp,ypn);
        score(xf, y, f);
        for (i = 1; i<=nvar; i++)
		{
            dfdv[(i-1)*20+j] = (f2[i] - f[i]) / delv2[iv];
        }
        v2[iv] = sav;
    }
    for (i = 1; i<=nvar; i++)
	{
        f[i] = f1[i] - f2[i];
        f1[i] = -f[i];
    }
	int det;
    ludcmp(dfdv, nvar, 20, indx, det);
    lubksb(dfdv, nvar, 20, indx, f1);
    j = 0;
    for (iv = 1; iv<=n2; iv++)
	{
        j = j + 1;
        v1[iv] = v1[iv] + f1[j];
        dv1[iv] = f1[j];
    }
    for (iv = 1; iv<=n1; iv++)
	{
        j = j + 1;
        v2[iv] = v2[iv] + f1[j];
        dv2[iv] = f1[j];
    }
}

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]));
}