d2r11.cpp

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

# include<iostream.h>
# include<math.h>
# include<iomanip.h>

void spline(double x[], double y[], int n, double yp1,double  ypn, double y2[])
{
    double u[100],aaa,sig,p,bbb,ccc,qn,un;
	int i,k;
    if( yp1 > 9.9e+29 )
	{
        y2[1] = 0;
        u[1] = 0;
	}
    else
	{
        y2[1] = -0.5;
        aaa = (y[2] - y[1]) / (x[2] - x[1]);
        u[1] = (3.0 / (x[2] - x[1])) * (aaa - yp1);
    }
    for (i = 2;i<=n-1;i++)
	{
        sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
        p = sig * y2[i - 1] + 2.0;
        y2[i] = (sig - 1.0) / p;
        aaa = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
        bbb = (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
        ccc = x[i + 1] - x[i - 1];
        u[i] = (6.0 * (aaa - bbb) / ccc - sig * u[i - 1]) / p;
    }
    if (ypn > 9.9e+29 )
	{
        qn = 0.0;
        un = 0.0;
	}
    else
	{
        qn = 0.5;
        aaa = ypn - (y[n] - y[n - 1]) / (x[n] - x[n - 1]);
        un = (3.0/ (x[n] - x[n - 1])) * aaa;
    }
    y2[n] = (un - qn * u[n - 1]) / (qn * y2[n - 1] + 1.0);
    for (k = n - 1;k>=1;k--)
        y2[k] = y2[k] * y2[k + 1] + u[k];
}

void splint(double xa[], double ya[], double y2a[], int n, double& x, double& y)
{
	int klo,khi,k;
	double h,a,b,aaa,bbb;
    klo = 1;
    khi = n;
loop:   if (khi - klo > 1 )
		{
			 k = (khi + klo) / 2;
			 if (xa[k] > x)
				khi = k;
			 else
			{
				klo = k;
			}
        goto loop;
		}
    h = xa[khi] - xa[klo];
    if (h == 0 )
	{
        cout<<"  pause  'bad  xa  input'"<<endl;
        return;
    }
    a = (xa[khi] - x) / h;
    b = (x - xa[klo]) / h;
    aaa = a * ya[klo] + b * ya[khi];
    bbb = (a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi];
    y = aaa + bbb * h*h /6.0;
}

void splin2(double x1a[], double x2a[], double ya[][11], double y2a[][11], int m, int n, double x1, double x2, double& y)
{ 
	double ytmp[100], y2tmp[100], yytmp[100];
	int j,k;
    for (j = 1; j<=m; j++)
	{
        for (k = 1; k<=n; k++)
		{
            ytmp[k] = ya[j][k];
            y2tmp[k] = y2a[j][k];
		}
        splint(x2a, ytmp, y2tmp, n, x2, yytmp[j]);
    }
    spline(x1a, yytmp, m, 1e+30, 1e+30, y2tmp);
    splint(x1a, yytmp, y2tmp, m, x1, y);
}

void splie2(double x1a[], double x2a[], double ya[][11], int m, int n, double y2a[][11])
{
    double ytmp[101], y2tmp[101];
	int j,k;
    for (j = 1; j<=m; j++)
	{
        for (k = 1; k<=n; k++)
            ytmp[k] = ya[j][k];
        spline(x2a, ytmp, n, 1e+30, 1e+30, y2tmp);
        for (k = 1; k<=n; k++)
            y2a[j][k] = y2tmp[k];
    }
}
 
void main()
{
    //program d2r11
    //driver for routine splin2
	int m,n,i,j;
    double x1[11],x2[11],y[11][11],y2[11][11],x1x2,xx1,xx2,ff,f;
	m = 10;
    n = 10;
	f = 0.0;
    for (i = 1; i<=m; i++)
        x1[i] = 0.2 * i;
    for (i = 1; i<=n; i++)
        x2[i] = 0.2 * i;
    
    for (i = 1; i<=m; i++)
	{
        for (j = 1; j<=n; j++)
		{
            x1x2 = x1[i] * x2[j];
            y[i][j] = x1x2 * exp(-x1x2);
        }
    }
    splie2(x1, x2, y, m, n, y2);
    cout<<endl;
    cout<<"        x1           x2         splin2       actual"<<endl;
    for (i = 1; i<=m; i++) 
	{
        xx1 = 0.1 * i;
        xx2 = xx1*xx1;
        splin2(x1, x2, y, y2, m, n, xx1, xx2, f);
        x1x2 = xx1 * xx2;
        ff = x1x2 * exp(-x1x2);
		cout<<setprecision(6)<<setiosflags(ios::fixed);
        cout<<setw(13)<<xx1;
        cout<<setw(13)<<xx2;
        cout<<setw(13)<<f;
        cout<<setw(13)<<ff<<endl;
    }
}