cubicspline.c

su 的源代码库

	else if (del1*del2<0.0) {
		dmax = 3.0*del1;
		if (ABS(yd[0][1])>ABS(dmax)) yd[0][1] = dmax;
	}

	/* loop over interior points */
	for (i=1; i<n-1; i++) {

		/* compute intervals and slopes */
		h1 = x[i]-x[i-1];
		h2 = x[i+1]-x[i];
		hsum = h1+h2;
		del1 = (y[i]-y[i-1])/h1;
		del2 = (y[i+1]-y[i])/h2;

		/* if not strictly monotonic, zero derivative */
		if (del1*del2<=0.0) {
			yd[i][1] = 0.0;
		
		/*
		 * else, if strictly monotonic, use Butland's formula:
		 *      3*(h1+h2)*del1*del2 
		 * -------------------------------
		 * ((2*h1+h2)*del1+(h1+2*h2)*del2)
		 * computed as follows to avoid roundoff error
		 */
		} else {
			hsum3 = hsum+hsum+hsum;
			w1 = (hsum+h1)/hsum3;
			w2 = (hsum+h2)/hsum3;
			dmin = MIN(ABS(del1),ABS(del2));
			dmax = MAX(ABS(del1),ABS(del2));
			drat1 = del1/dmax;
			drat2 = del2/dmax;
			yd[i][1] = dmin/(w1*drat1+w2*drat2);
		}
	}

	/* set right end derivative via shape-preserving 3-point formula */
	w1 = -h2/hsum;
	w2 = (h2+hsum)/hsum;
	yd[n-1][1] = w1*del1+w2*del2;
	if (yd[n-1][1]*del2<=0.0)
		yd[n-1][1] = 0.0;
	else if (del1*del2<0.0) {
		dmax = 3.0*del2;
		if (ABS(yd[n-1][1])>ABS(dmax)) yd[n-1][1] = dmax;
	}

	/* compute 2nd and 3rd derivatives of cubic polynomials */
	for (i=0; i<n-1; i++) {
		h2 = x[i+1]-x[i];
		del2 = (y[i+1]-y[i])/h2;
		divdf3 = yd[i][1]+yd[i+1][1]-2.0*del2;
		yd[i][2] = 2.0*(del2-yd[i][1]-divdf3)/h2;
		yd[i][3] = (divdf3/h2)*(6.0/h2);
	}
	yd[n-1][2] = yd[n-2][2]+(x[n-1]-x[n-2])*yd[n-2][3];
	yd[n-1][3] = yd[n-2][3];
}

void csplin (int n, float x[], float y[], float yd[][4])
/*****************************************************************************
compute cubic spline interpolation coefficients for interpolation with
continuous second derivatives
******************************************************************************
Input:
n		number of samples
x  		array[n] of monotonically increasing or decreasing abscissae
y		array[n] of ordinates

Output:
yd		array[n][4] of cubic interpolation coefficients (see notes)
******************************************************************************
Notes:
The computed cubic spline coefficients are as follows:
yd[i][0] = y(x[i])    (the value of y at x = x[i])
yd[i][1] = y'(x[i])   (the 1st derivative of y at x = x[i])
yd[i][2] = y''(x[i])  (the 2nd derivative of y at x = x[i])
yd[i][3] = y'''(x[i]) (the 3rd derivative of y at x = x[i])

To evaluate y(x) for x between x[i] and x[i+1] and h = x-x[i],
use the computed coefficients as follows:
y(x) = yd[i][0]+h*(yd[i][1]+h*(yd[i][2]/2.0+h*yd[i][3]/6.0))
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 10/03/89
Modified:  Dave Hale, Colorado School of Mines, 02/28/91
	changed to work for n=1.
Modified:  Dave Hale, Colorado School of Mines, 08/04/91
	fixed bug in computation of left end derivative
*****************************************************************************/
{
	int i;
	float h1,h2,del1,del2,dmax,hsum,w1,w2,divdf3,sleft,sright,alpha,t;

	/* if n=1, then use constant interpolation */
	if (n==1) {
		yd[0][0] = y[0];
		yd[0][1] = 0.0;
		yd[0][2] = 0.0;
		yd[0][3] = 0.0;
		return;

	/* else, if n=2, then use linear interpolation */
	} else if (n==2) {
		yd[0][0] = y[0];  yd[1][0] = y[1];
		yd[0][1] = yd[1][1] = (y[1]-y[0])/(x[1]-x[0]);
		yd[0][2] = yd[1][2] = 0.0;
		yd[0][3] = yd[1][3] = 0.0;
		return;
	}
	
	/* set left end derivative via shape-preserving 3-point formula */
	h1 = x[1]-x[0];
	h2 = x[2]-x[1];
	hsum = h1+h2;
	del1 = (y[1]-y[0])/h1;
	del2 = (y[2]-y[1])/h2;
	w1 = (h1+hsum)/hsum;
	w2 = -h1/hsum;
	sleft = w1*del1+w2*del2;
	if (sleft*del1<=0.0)
		sleft = 0.0;
	else if (del1*del2<0.0) {
		dmax = 3.0*del1;
		if (ABS(sleft)>ABS(dmax)) sleft = dmax;
	}

	/* set right end derivative via shape-preserving 3-point formula */
	h1 = x[n-2]-x[n-3];
	h2 = x[n-1]-x[n-2];
	hsum = h1+h2;
	del1 = (y[n-2]-y[n-3])/h1;
	del2 = (y[n-1]-y[n-2])/h2;
	w1 = -h2/hsum;
	w2 = (h2+hsum)/hsum;
	sright = w1*del1+w2*del2;
	if (sright*del2<=0.0)
		sright = 0.0;
	else if (del1*del2<0.0) {
		dmax = 3.0*del2;
		if (ABS(sright)>ABS(dmax)) sright = dmax;
	}
	
	/* compute tridiagonal system coefficients and right-hand-side */	
	yd[0][0] = 1.0;
	yd[0][2] = 2.0*sleft;
	for (i=1; i<n-1; i++) {
		h1 = x[i]-x[i-1];
		h2 = x[i+1]-x[i];
		del1 = (y[i]-y[i-1])/h1;
		del2 = (y[i+1]-y[i])/h2;
		alpha = h2/(h1+h2);
		yd[i][0] = alpha;
		yd[i][2] = 3.0*(alpha*del1+(1.0-alpha)*del2);
	}
	yd[n-1][0] = 0.0;
	yd[n-1][2] = 2.0*sright;
	
	/* solve tridiagonal system for slopes */
	t = 2.0;
	yd[0][1] = yd[0][2]/t;
	for (i=1; i<n; i++) {
		yd[i][3] = (1.0-yd[i-1][0])/t;
		t = 2.0-yd[i][0]*yd[i][3];
		yd[i][1] = (yd[i][2]-yd[i][0]*yd[i-1][1])/t;
	}
	for (i=n-2; i>=0; i--)
		yd[i][1] -= yd[i+1][3]*yd[i+1][1];

	/* copy ordinates into output array */
	for (i=0; i<n; i++)
		yd[i][0] = y[i];

	/* compute 2nd and 3rd derivatives of cubic polynomials */
	for (i=0; i<n-1; i++) {
		h2 = x[i+1]-x[i];
		del2 = (y[i+1]-y[i])/h2;
		divdf3 = yd[i][1]+yd[i+1][1]-2.0*del2;
		yd[i][2] = 2.0*(del2-yd[i][1]-divdf3)/h2;
		yd[i][3] = (divdf3/h2)*(6.0/h2);
	}
	yd[n-1][2] = yd[n-2][2]+(x[n-1]-x[n-2])*yd[n-2][3];
	yd[n-1][3] = yd[n-2][3];
}
#include <cwp.h>

void chermite (int n, float x[], float y[], float yd[][4])
/*****************************************************************************
Compute cubic interpolation coefficients via Hermite polynomial
******************************************************************************
Input:
n		number of samples
x  		array[n] of monotonically increasing or decreasing abscissae
y		array[n] of ordinates

Output:
yd		array[n][4] of cubic interpolation coefficients (see notes)
******************************************************************************
Notes:
The computed cubic spline coefficients are as follows:
yd[i][0] = y(x[i])    (the value of y at x = x[i])
yd[i][1] = y'(x[i])   (the 1st derivative of y at x = x[i])
yd[i][2] = y''(x[i])  (the 2nd derivative of y at x = x[i])
yd[i][3] = y'''(x[i]) (the 3rd derivative of y at x = x[i])

To evaluate y(x) for x between x[i] and x[i+1] and h = x-x[i],
use the computed coefficients as follows:
y(x) = yd[i][0]+h*(yd[i][1]+h*(yd[i][2]/2.0+h*yd[i][3]/6.0))

*****************************************************************************/
{
	int i;
	float dx ,f1,f3;

	/* copy ordinates into output array */
	for (i=0; i<n; i++){
		yd[i][0] = y[i];
        }

	/* if n=1, then use constant interpolation */
	if (n==1) {
		yd[0][1] = 0.0;
		yd[0][2] = 0.0;
		yd[0][3] = 0.0;
		return;

	/* else, if n=2, then use linear interpolation */
	} else if (n==2) {
		yd[0][1] = yd[1][1] = (y[1]-y[0])/(x[1]-x[0]);
		yd[0][2] = yd[1][2] = 0.0;
		yd[0][3] = yd[1][3] = 0.0;
		return;
	}



        /* compute forward & backwards differences at the ends */

        yd[0][1]   = (y[2] - y[0])     / (x[2] - x[0]);
        yd[n-1][1] = (y[n-1] - y[n-2]) / (x[n-1] - x[n-2]);

        /* compute central differences in between */

        for( i=1; i<n-1; i++ ){

           yd[i][1] = 0.5 * (y[i] - y[i-1]) / (x[i] - x[i-1])
	            + 0.5 * (y[i+1] - y[i]) / (x[i+1] - x[i]);
        }

        /* calculate 2nd and 3rd derivatives */

        for( i=0; i<n-1; i++ ){

           dx       = x[i+1] - x[i];
           f1       = (y[i+1] - y[i]) / dx;
           f3       = yd[i][1] + yd[i+1][1] - 2.0 * f1;
           yd[i][2] = 2.0*(f1 - yd[i][1] - f3) / dx;
           yd[i][3] = 6.0*f3 / (dx*dx);

        }

}