1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
|
#include <stdio.h>
#include <math.h>
#define PI 3.14159265
void spline(float *, float *, int, float, float, float *);
void splint(float *, float *, float *, int, float, float *);
void main(void)
{
int i, nfunc, np;
float f, x, y, yp1, ypn, *xa, *ya, *y2;
np = 10; /* number of points */
xa = (float *) malloc((np+1)*sizeof(float));
ya = (float *) malloc((np+1)*sizeof(float));
y2 = (float *) malloc((np+1)*sizeof(float));
printf("\nsine function from 0 to pi\n");
for (i=1; i<=np; i++) {
xa[i] = i*PI/np;
ya[i] = sin(xa[i]);
}
yp1 = cos(xa[1]);
ypn = cos(xa[np]);
/* call spline to get second derivatives */
spline(xa, ya, np, yp1, ypn, y2);
/* call splint for interpolations */
printf("\n%9s %13s %17s\n", "x", "f(x)", "interpolation");
for (i=1; i<=10; i++) {
x = (-0.05+i/10.0)*PI;
f = sin(x);
splint(xa, ya, y2, np, x, &y);
printf("%12.6f %12.6f %12.6f\n", x, f, y);
}
}
void spline(float *x, float *y, int n, float yp1, float ypn, float *y2)
{
int i, k;
float p, qn, sig, un, *u;
u = (float *) malloc(n*sizeof(float));
if (yp1 > 0.99e30)
y2[1] = u[1] = 0.0;
else {
y2[1] = -0.5;
u[1] = (3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-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;
u[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
u[i] = (6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
}
if (ypn > 0.99e30)
qn = un = 0.0;
else {
qn = 0.5;
un = (3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
}
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(float *xa, float *ya, float *y2a, int n, float x, float *y)
{
int klo, khi, k;
float h, b, a;
klo=1; khi=n;
while (khi-klo > 1) {
k=(khi+klo) >> 1;
if (xa[k] > x) khi=k;
else klo=k;
}
h = xa[khi]-xa[klo];
if (h == 0.0)
printf("Error: Bad xa input to routine splint.\n");
a = (xa[khi]-x)/h;
b = (x-xa[klo])/h;
*y = a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
}
|