d2r3.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;
        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 main()
{
	double x[21], y[21], y2[21],pi,yp1,ypn ;
	int n, i;
    //program d2r3
    //driver for routine spline
    n = 20;
    pi = 3.1415926;
    //生成插值数组
    for (i = 1;i<=20;i++)
	{
        x[i] = i * pi / n;
        y[i] = sin(x[i]);
    }
    //用三次样条函数计算二阶导数
    yp1 = cos(x[1]);
    ypn = cos(x[n]);
    cout<<endl;   
    cout<<"Second-derivative for sin(x) from 0 to pi"<<endl;
    spline(x, y, n, yp1, ypn, y2);
    //验证计算结果
    cout<<"                     spline        actual"<<endl;
    cout<<"        angle      2nd deriv     2nd deriv"<<endl;
    for( i = 1;i<=n;i++)
	{
		cout<<setprecision(6)<<setiosflags(ios::fixed);
        cout<<setw(14)<<x[i];
        cout<<setw(14)<<y2[i];
        cout<<setw(14)<<-sin(x[i])<<endl;
    } 
}