spline.c revision 7c478bd95313f5f23a4c958a745db2134aa03244
/*
* CDDL HEADER START
*
* The contents of this file are subject to the terms of the
* Common Development and Distribution License, Version 1.0 only
* (the "License"). You may not use this file except in compliance
* with the License.
*
* You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
* See the License for the specific language governing permissions
* and limitations under the License.
*
* When distributing Covered Code, include this CDDL HEADER in each
* file and include the License file at usr/src/OPENSOLARIS.LICENSE.
* If applicable, add the following below this CDDL HEADER, with the
* fields enclosed by brackets "[]" replaced with your own identifying
* information: Portions Copyright [yyyy] [name of copyright owner]
*
* CDDL HEADER END
*/
/* Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T */
/* All Rights Reserved */
/* Portions Copyright (c) 1988, Sun Microsystems, Inc. */
/* All Rights Reserved. */
#ident "%Z%%M% %I% %E% SMI" /* SVr4.0 1.5 */
#include <stdio.h>
#include <math.h>
#define NP 1000
#define REAL double
int n;
int auta;
int periodic;
/* Spline fit technique
let x,y be vectors of abscissas and ordinates
h be vector of differences h9i8=x9i8-x9i-1988
y" be vector of 2nd derivs of approx function
If the points are numbered 0,1,2,...,n+1 then y" satisfies
(R W Hamming, Numerical Methods for Engineers and Scientists,
2nd Ed, p349ff)
h9i8y"9i-1988+2(h9i8+h9i+18)y"9i8+h9i+18y"9i+18
= 6[(y9i+18-y9i8)/h9i+18-(y9i8-y9i-18)/h9i8] i=1,2,...,n
where y"908 = y"9n+18 = 0
This is a symmetric tridiagonal system of the form
| a918 h928 | |y"918| |b918|
| h928 a928 h938 | |y"928| |b928|
| h938 a938 h948 | |y"938| = |b938|
| . | | .| | .|
| . | | .| | .|
It can be triangularized into
| d918 h928 | |y"918| |r918|
| d928 h938 | |y"928| |r928|
| d938 h948 | |y"938| = |r938|
| . | | .| | .|
| . | | .| | .|
where
d918 = a918
r908 = 0
d9i8 = a9i8 - h9i8829/d9i-18 1<i<_n
r9i8 = b9i8 - h9i8r9i-18/d9i-1i8 1<_i<_n
the back solution is
y"9n8 = r9n8/d9n8
y"9i8 = (r9i8-h9i+18y"9i+18)/d9i8 1<_i<n
superficially, d9i8 and r9i8 don't have to be stored for they can be
recalculated backward by the formulas
d9i-18 = h9i8829/(a9i8-d9i8) 1<i<_n
r9i-18 = (b9i8-r9i8)d9i-18/h9i8 1<i<_n
unhappily it turns out that the recursion forward for d
is quite strongly geometrically convergent--and is wildly
unstable going backward.
There's similar trouble with r, so the intermediate
results must be kept.
Note that n-1 in the program below plays the role of n+1 in the theory
Other boundary conditions_________________________
The boundary conditions are easily generalized to handle
y908" = ky918", y9n+18" = ky9n8"
for some constant k. The above analysis was for k = 0;
k = 1 fits parabolas perfectly as well as stright lines;
k = 1/2 has been recommended as somehow pleasant.
All that is necessary is to add h918 to a918 and h9n+18 to a9n8.
Periodic case_____________
To do this, add 1 more row and column thus
| a918 h928 h918 | |y918"| |b918|
| h928 a928 h938 | |y928"| |b928|
| h938 a948 h948 | |y938"| |b938|
| | | .| = | .|
| . | | .| | .|
| h918 h908 a908 | | .| | .|
where h908=_ h9n+18
The same diagonalization procedure works, except for
the effect of the 2 corner elements. Let s9i8 be the part
of the last element in the i8th9 "diagonalized" row that
arises from the extra top corner element.
s918 = h918
s9i8 = -s9i-18h9i8/d9i-18 2<_i<_n+1
After "diagonalizing", the lower corner element remains.
Call t9i8 the bottom element that appears in the i8th9 colomn
as the bottom element to its left is eliminated
t918 = h918
t9i8 = -t9i-18h9i8/d9i-18
Evidently t9i8 = s9i8.
Elimination along the bottom row
introduces further corrections to the bottom right element
and to the last element of the right hand side.
Call these corrections u and v.
u918 = v918 = 0
u9i8 = u9i-18-s9i-18*t9i-18/d9i-18
v9i8 = v9i-18-r9i-18*t9i-18/d9i-18 2<_i<_n+1
The back solution is now obtained as follows
y"9n+18 = (r9n+18+v9n+18)/(d9n+18+s9n+18+t9n+18+u9n+18)
y"9i8 = (r9i8-h9i+18*y9i+18-s9i8*y9n+18)/d9i8 1<_i<_n
Interpolation in the interval x9i8<_x<_x9i+18 is by the formula
y = y9i8x9+8 + y9i+18x9-8 -(h8299i+18/6)[y"9i8(x9+8-x9+8839)+y"9i+18(x9-8-x9-8839)]
where
x9+8 = x9i+18-x
x9-8 = x-x9i8
*/
rhs(i){
int i_;
double zz;
i_ = i==n-1?0:i;
}
spline(){
REAL h;
int end;
int i,j,m;
if(n<3) return(0);
d = 1;
r[0] = 0;
s = periodic?-1:0;
for(i=0;++i<n-!periodic;){ /* triangularize */
u = i==1?zero:u-s*s/d;
s = -hi*s/d;
diag[i] = d = i==1? a:
}
for(i=n-!periodic;--i>=0;){ /* back substitute */
end = i==n-1;
if(i>0){
if(i>1){
d = diag[i-1];
s = -s*d/hi;
}}
if(!periodic) {
}
if(end) continue;
if(m<=0) m = 1;
h = hi1/m;
for(j=m;j>0||i==0&&j==0;j--){ /* interpolate */
}
}
return(1);
}
readin() {
for(n=0;n<NP;n++){
getfloat(p)
REAL *p;{
char buf[30];
register c;
int i;
extern double atof();
for(;;){
c = getchar();
if (c==EOF) {
*buf = '\0';
return(0);
}
*buf = c;
switch(*buf){
case ' ':
case '\t':
case '\n':
continue;}
break;}
for(i=1;i<30;i++){
c = getchar();
if (c==EOF) {
buf[i] = '\0';
break;
}
buf[i] = c;
if('0'<=c && c<='9') continue;
switch(c) {
case '.':
case '+':
case '-':
case 'E':
case 'e':
continue;}
break; }
buf[i] = ' ';
return(1); }
getlim(p)
struct proj *p; {
int i;
for(i=0;i<n;i++) {
}
char *argv[];{
extern char *malloc();
int i;
while(--argc > 0) {
argv++;
case '-':
argv[0]++;
goto again;
case 'a':
auta = 1;
break;
case 'k':
break;
case 'n':
break;
case 'p':
periodic = 1;
break;
case 'x':
x.lbf = 1;
x.ubf = 1;
break;
default:
exit(1);
}
}
readin();
getlim(&x);
getlim(&y);
i = (n+1)*sizeof(dx);
exit(0);
}
int *argcp;
char ***argvp;{
double atof();
char c;
if(*argcp<=1) return(0);
c = (*argvp)[1][0];
if(!('0'<=c&&c<='9' || c=='-' || c== '.' )) return(0);
(*argcp)--;
(*argvp)++;
return(1); }