sfroid.cpp

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

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

double x[42], h, c2, anorm;
int mm, n;

void pinvs(int ie1, int ie2, int je1, int jsf, int jc1, int k, double c[], int nci, int ncj, int nck, double s[], int nsi, int nsj)
{
    double zero = 0.0; int one = 1;
    double pscl[11];
	int indxr[11];
    int je2 = je1 + ie2 - ie1;
    int i,j,js1 = je2 + 1;
    for (i = ie1; i<=ie2; i++)
	{
        double big = zero;
        for (j = je1; j<=je2; j++)
		{
            if (fabs(s[(i-1)*nsj+j]) > big)
			{
				big = fabs(s[(i-1)*nsj+j]);
			}
        }
        if (big == zero)
		{
			cout<<"singular matrix, row all 0"<<endl;
		}
        pscl[i] = one / big;
        indxr[i] = 0;
    }
	int id,jp,ipiv,jpiv;
	double pivinv,big;
    for (id = ie1; id<=ie2; id++)
	{
        double piv = zero;
        for (i = ie1; i<=ie2; i++)
		{
            if (indxr[i] == 0)
			{
                big = zero;
                for (j = je1; j<=je2; j++)
				{
                    if (fabs(s[(i-1)*nsj+j]) > big)
					{
                        jp = j;
                        big = fabs(s[(i-1)*nsj+j]);
                    }
                }
                if (big * pscl[i] > piv)
				{
                    ipiv = i;
                    jpiv = jp;
                    piv = big * pscl[i];
                }
            }
        }
        if (s[(ipiv-1)*nsj+jpiv] == zero)
		{
			cout<<"singular matrix";
		}
        indxr[ipiv] = jpiv;
        pivinv = (double)one / s[(ipiv-1)*nsj+jpiv];
        for (j = je1; j<=jsf; j++)
		{
            s[(ipiv-1)*nsj+j] = s[(ipiv-1)*nsj+j] * pivinv;
        }
        s[(ipiv-1)*nsj+jpiv] = (double)one;
        for (i = ie1; i<=ie2; i++)
		{
            if (indxr[i] != jpiv)
			{
                if (s[(i-1)*nsj+jpiv] != zero)
				{
                    double dum = s[(i-1)*nsj+jpiv];
                    for (j = je1; j<=jsf; j++)
					{
                        s[(i-1)*nsj+j] = s[(i-1)*nsj+j] - dum * s[(ipiv-1)*nsj+j];
                    }
                    s[(i-1)*nsj+jpiv] = zero;
                }
            }
        }
    }
    int jcoff = jc1 - js1;
    int icoff = ie1 - je1;
    for (i = ie1; i<=ie2; i++)
	{
        int irow = indxr[i] + icoff;
        for (j = js1; j<=jsf; j++)
		{
            c[((irow-1)*ncj+ j + jcoff-1)*nck+ k] = s[(i-1)*nsj+j];
        }
    }
}

void red(int iz1, int iz2, int jz1, int jz2, int jm1, int jm2, int jmf,
		 int ic1, int jc1, int jcf, int kc, double c[], int nci, int ncj,
		 int nck, double s[], int nsi, int nsj)
{
    int loff = jc1 - jm1;
    int j,l,i,ic = ic1;
	double vx;
    for (j = jz1; j<=jz2; j++)
	{
        for (l = jm1; l<=jm2; l++)
		{
            vx = c[((ic - 1) * ncj + l + loff - 1) * nck + kc];
            for (i = iz1; i<=iz2; i++)
			{
                s[(i-1)*nsj+l] = s[(i-1)*nsj+l] - s[(i-1)*nsj+j] * vx;
            }
        }
        vx = c[((ic - 1) * ncj + jcf - 1) * nck + kc];
        for (i = iz1; i<=iz2; i++)
		{
            s[(i-1)*nsj+jmf] = s[(i-1)*nsj+jmf] - s[(i-1)*nsj+j] * vx;
        }
        ic = ic + 1;
    }
}

void difeq(int k, int k1, int k2, int jsf, int is1, int isf, int indexv[],
		   int ne, double s[], int nsi, int nsj, double y[], int nyj, int nyk)
{
    int m = 41;
    if (k == k1)
	{
        if (n + mm % 2 == 1)
		{
            s[2*nsj+3 + indexv[1]] = 1.0;  //方程(15-37)
            s[2*nsj+3 + indexv[2]] = 0.0;
            s[2*nsj+3 + indexv[3]] = 0.0;
            s[2*nsj+jsf] = y[1];  //方程(15-31)
		}
        else
		{
            s[2*nsj+3 + indexv[1]] = 0.0;  //方程(15-37)
            s[2*nsj+3 + indexv[2]] = 1.0;
            s[2*nsj+3 + indexv[3]] = 0.0;
            s[2*nsj+jsf] = y[nyk+1];  //方程(15-31)
        }
	}
    else
	{
		if (k > k2)
		{
			s[3 + indexv[1]] = -(y[2*nyk+m] - c2) / (2.0 * double(mm + 1.0));  //方程(15-38)
			s[3 + indexv[2]] = 1.0;
			s[3 + indexv[3]] = -y[m] / (2.0 * double(mm + 1.0));
			s[jsf] = y[nyk+m] - (y[2*nyk+m] - c2) * y[m] / (2.0 * double(mm + 1.0)); //方程(15-32)
			s[nsj + 3 + indexv[1]] = 1.0;  //方程(15-39)
			s[nsj + 3 + indexv[2]] = 0.0;
			s[nsj + 3 + indexv[3]] = 0.0;
			s[nsj + jsf] = y[m] - anorm;  //方程(15-33)
		}
		else
		{
			s[indexv[1]] = -1.0;  //方程(15-34)
			s[indexv[2]] = -0.5 * h;
			s[indexv[3]] = 0.0;
			s[3 + indexv[1]] = 1.0;
			s[3 + indexv[2]] = -0.5 * h;
			s[3 + indexv[3]] = 0.0;
			double temp = h / (1.0 - pow((x[k] + x[k - 1]) , 2) * 0.25);
			double temp2 = 0.5 * (y[2*nyk+k] + y[2*nyk+k - 1]) - c2 * 0.25 * pow((x[k] + x[k - 1]) , 2);
			s[nsj+indexv[1]] = temp * temp2 * 0.5;  //方程(15-35)
			s[nsj+indexv[2]] = -1.0 - 0.5 * temp * double(mm + 1.0) * (x[k] + x[k - 1]);
			s[nsj+indexv[3]] = 0.25 * temp * (y[k] + y[k - 1]);
			s[nsj+3 + indexv[1]] = s[nsj+indexv[1]];
			s[nsj+3 + indexv[2]] = 2.0 + s[nsj+indexv[2]];
			s[nsj+3 + indexv[3]] = s[nsj+indexv[3]];
			s[2*nsj+indexv[1]] = 0.0;  //方程(15-36)
			s[2*nsj+indexv[2]] = 0.0;
			s[2*nsj+indexv[3]] = -1.0;
			s[2*nsj+3 + indexv[1]] = 0.0;
			s[2*nsj+3 + indexv[2]] = 0.0;
			s[2*nsj+3 + indexv[3]] = 1.0;
			s[jsf] = y[k] - y[k - 1] - 0.5 * h * (y[nyk+k] + y[nyk+k - 1]);  //方程(15-26)
			double dum = (x[k] + x[k - 1]) * 0.5 * double(mm + 1.0) * (y[nyk+k] + y[nyk+k - 1]);    //方程(15-27)
			dum = dum - temp2 * 0.5 * (y[k] + y[k - 1]);
			s[nsj+jsf] = y[nyk+k] - y[nyk+k - 1] - temp * dum;
			s[2*nsj+jsf] = y[2*nyk+k] - y[2*nyk+k - 1];  //方程(15-30)
		}
	}
}

void bksub(int ne, int nb, int jf, int k1, int k2, double c[], int nci, int ncj, int nck)
{
    int k,j,i,kp,nbf = ne - nb;
	double xx;
    for (k = k2; k>=k1; k--)
	{
        kp = k + 1;
        for (j = 1; j<=nbf; j++)
		{
            xx = c[((j-1)*ncj+jf-1)*nck+kp];
            for (i = 1; i<=ne; i++)
			{
                c[((i-1)*ncj+jf-1)*nck+k] = c[((i-1)*ncj+jf-1)*nck+k] - c[((i-1)*ncj+j-1)*nck+k] * xx;
            }
        }
    }
    for (k = k1; k<=k2; k++)
	{
        kp = k + 1;
        for (i = 1; i<=nb; i++)
		{
            c[(i-1)*ncj*nck+k] = c[((i + nbf-1)*ncj+jf-1)*nck+k];
		}
        for (i = 1; i<=nbf; i++)
		{
            c[(i + nb-1)*ncj*nck+k] = c[((i-1)*ncj+jf-1)*nck+kp];
        }
    }
}

void solvde(int itmax, double conv, double slowc, double scalv[], int indexv[],
			int ne, int nb, int m, double y[], int nyj, int nyk, double c[],
			int nci, int ncj, int nck, double s[], int nsi, int nsj)
{
    double ermax[11];
	int kmax[11];
    int k1 = 1;
    int k2 = m;
    int nvars = ne * m;
    int j1 = 1;
    int j2 = nb;
    int j3 = nb + 1;
    int j4 = ne;
    int j5 = j4 + j1;
    int j6 = j4 + j2;
    int j7 = j4 + j3;
    int j8 = j4 + j4;
    int j9 = j8 + j1;
    int ic1 = 1;
    int ic2 = ne - nb;
    int ic3 = ic2 + 1;
    int ic4 = ne;
    int jc1 = 1;
    int k,j,kp,it,jcf = ic3;
    for (it = 1; it<=itmax; it++)
	{
        k = k1;
        difeq(k, k1, k2, j9, ic3, ic4, indexv, ne, s, nsi, nsj, y, nyj, nyk);
        pinvs(ic3, ic4, j5, j9, jc1, k1, c, nci, ncj, nck, s, nsi, nsj);
        for (k = k1 + 1; k<=k2; k++)
		{
            kp = k - 1;
            difeq(k, k1, k2, j9, ic1, ic4, indexv, ne, s, nsi, nsj, y, nyj, nyk);
            red(ic1, ic4, j1, j2, j3, j4, j9, ic3, jc1, jcf, kp, c, nci, ncj, nck, s, nsi, nsj);
			pinvs(ic1, ic4, j3, j9, jc1, k, c, nci, ncj, nck, s, nsi, nsj);
        }
        k = k2 + 1;
        difeq(k, k1, k2, j9, ic1, ic2, indexv, ne, s, nsi, nsj, y, nyj, nyk);
        red(ic1, ic2, j5, j6, j7, j8, j9, ic3, jc1, jcf, k2, c, nci, ncj, nck, s, nsi, nsj);
        pinvs(ic1, ic2, j7, j9, jcf, k2 + 1, c, nci, ncj, nck, s, nsi, nsj);
        bksub(ne, nb, jcf, k1, k2, c, nci, ncj, nck);
        double erq = 0.0;
		int jv,km;
		double errj,vmax;
        for (j = 1; j<=ne; j++)
		{
            jv = indexv[j];
            ermax[j] = 0.0;
            errj = 0.0;
            kmax[j] = 0;
            vmax = 0.0;
            for (k = k1; k<=k2; k++)
			{
				double temp;
                double vz = fabs(c[((j-1)*ncj)*nck+k]);
				temp=vz-vmax;
                if (vz > vmax)
				{
                    vmax = vz;
                    km = k;
                }
                errj = errj + vz;
            }
            erq = erq + errj / scalv[jv];
            ermax[j] = c[(j-1)*ncj*nck+km] / scalv[jv];
            kmax[j] = km;
        }
        erq = erq / nvars;
		double dum;
        if (erq > slowc)
		{
            dum = erq;
		}
        else
		{
            dum = slowc;
        }
		double fac;
        fac = slowc / dum;
        for (jv = 1; jv<=ne; jv++)
		{
            j = indexv[jv];
            for (k = k1; k<=k2; k++)
			{
                y[(j-1)*nyk+k] = y[(j-1)*nyk+k] - fac * c[((jv-1)*ncj)*nck+k];
            }
        }
		cout<<"  ";
		cout.width(1);
        cout<<it;
		cout.width(12);
		cout<<erq;
		cout.width(14);
        cout<<fac;
        cout<<endl;
        for (j = 1; j<=ne; j++)
		{
            cout<<"   ";
			cout.width(12);
			cout<<kmax[j];
			cout.width(14);
			cout<<ermax[j]<<endl;
        }
        cout<<endl;
        if (erq < conv)
		{
			break;
		}
    }
}

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 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];
	int indexv[4];
	double y[124], c[379], s[22];
    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);
}