d2r4.cpp

Visual 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.0;
          u[1] = 0.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] + double(2);
          y2[i] = (sig - double(1)) / 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] = (double(6) * (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 = (double(3)/ (x[n] - x[n - 1])) * aaa;
      }
      y2[n] = (un - qn * u[n - 1]) / (qn * y2[n - 1] + 1);
      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 main()
{
    //program d2r4
    //driver for routine splint,which calls spline
    int nfunc,i,np;
    double xa[11], ya[11], y2[11],pi,yp1,ypn,x,f,y;
	np = 10;
    pi = 3.141593;
    for( nfunc = 1; nfunc<=2; nfunc++)
	{
        if (nfunc == 1 )
		{
            cout<<endl;
            cout<<"Sine function from 0 to pi"<<endl;
            for (i = 1;i<=np;i++)
			{
                xa[i] = double(i * pi) /double( np);
                ya[i] = sin(xa[i]);
            }
            yp1 = cos(xa[1]);
            ypn = cos(xa[np]);
		}
        else 
		{
			if (nfunc == 2)
			{
				cout<<endl;
				cout<<"Exponential function from 0 to 1"<<endl;
				for (i = 1;i<=np;i++)
				{
					xa[i] = double(1* i) /double(np);
					ya[i] = exp(xa[i]);
				}
				yp1 = exp(xa[1]);
				ypn = exp(xa[np]);
			}
           
		}
		//call spline to get second derivatives
		spline(xa, ya, np, yp1, ypn, y2);
		//call splint for interpolations
		cout<< "       x            f(x)     interpolation"<<endl;
		for (i = 1; i<=10; i++)
		{
			if (nfunc == 1)
			{
				x = (-0.05 + double(i) / 10.0) * pi;
				f = sin(x);
			}
			else
			{
				if( nfunc == 2 )
				{
					x = (-0.05 + double(i) / 10.0);
					f = exp(x);
				}
			}
			splint(xa, ya, y2, np, x, y);
			cout<<setprecision(6)<<setiosflags(ios::fixed);
			cout<<setw(12)<<x;
			cout<<setw(14)<<f;
			cout<<setw(14)<<y<<endl;
		}
        cout<<"**********************************"<<endl;
    }
}