static char *sccsid
= "@(#)spline.c 4.3 (Berkeley) %G%";
struct proj
{ int lbf
,ubf
; float a
,b
,lb
,ub
,quant
,mult
,val
[NP
]; } x
,y
;
let x,y be vectors of abscissas and ordinates
h be vector of differences hi=xi-xi-1
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,
hiy"i-1+2(hi+hi+1)y"i+hi+1y"i+1
= 6[(yi+1-yi)/hi+1-(yi-yi-1)/hi] i=1,2,...,n
This is a symmetric tridiagonal system of the form
| h3 a3 h4 | |y"3| = |b3|
It can be triangularized into
di = ai - hi2/di-1 1<i<_n
ri = bi - hiri-1/di-1\ei 1<_i<_n
y"i = (ri-hi+1y"i+1)/di 1<_i<n
superficially, di and ri don't have to be stored for they can be
recalculated backward by the formulas
di-1 = hi2/(ai-di) 1<i<_n
ri-1 = (bi-ri)di-1/hi 1<i<_n
unhappily it turns out that the recursion forward for d
is quite strongly geometrically convergent--and is wildly
There's similar trouble with r, so the intermediate
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
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 h1 to a1 and hn+1 to an.
Periodic case_____________
To do this, add 1 more row and column thus
The same diagonalization procedure works, except for
the effect of the 2 corner elements. Let si be the part
of the last element in the ith "diagonalized" row that
arises from the extra top corner element.
si = -si-1hi/di-1 2<_i<_n+1
After "diagonalizing", the lower corner element remains.
Call ti the bottom element that appears in the ith colomn
as the bottom element to its left is eliminated
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.
vi = vi-1-ri-1*ti-1/di-1 2<_i<_n+1
The back solution is now obtained as follows
y"n+1 = (rn+1+vn+1)/(dn+1+sn+1+tn+1+un+1)
y"i = (ri-hi+1*yi+1-si*yn+1)/di 1<_i<_n
Interpolation in the interval xi<_x<_xi+1 is by the formula
y = yix+ + yi+1x- -(h2i+1/6)[y"i(x+-x+3)+y"i+1(x--x-3)]
zz
= (y
.val
[i
]-y
.val
[i
-1])/(x
.val
[i
]-x
.val
[i
-1]);
return(6*((y
.val
[i_
+1]-y
.val
[i_
])/(x
.val
[i
+1]-x
.val
[i
]) - zz
));
float D2yi
,D2yi1
,D2yn1
,x0
,x1
,yy
,a
;
for(i
=0;++i
<n
-!periodic
;){ /* triangularize */
hi
= x
.val
[i
]-x
.val
[i
-1];
hi1
= i
==n
-1?x
.val
[1]-x
.val
[0]:
v
= i
==1?zero
:v
-s
*r
[i
-1]/d
;
r
[i
] = rhs(i
)-hi
*r
[i
-1]/d
;
if(i
==n
-2) a
+= konst
*hi1
;
for(i
=n
-!periodic
;--i
>=0;){ /* back substitute */
hi1
= end
?x
.val
[1]-x
.val
[0]:
hi
= x
.val
[i
]-x
.val
[i
-1];
D2yi
= (end
*v
+r
[i
]-hi1
*D2yi1
-s
*D2yn1
)/
if(i
==n
-2) a
+= konst
*hi1
;
if(i
==0) D2yi
= konst
*D2yi1
;
if(i
==n
-2) D2yi1
= konst
*D2yi
;
m
= 1.001*m
*hi1
/(x
.ub
-x
.lb
);
for(j
=m
;j
>0||i
==0&&j
==0;j
--){ /* interpolate */
yy
= D2yi
*(x0
-x0
*x0
*x0
)+D2yi1
*(x1
-x1
*x1
*x1
);
yy
= y
.val
[i
]*x0
+y
.val
[i
+1]*x1
-hi1
*hi1
*yy
/6;
printf("%f ",x
.val
[i
]+j
*h
);
if(auta
) x
.val
[n
] = n
*dx
+x
.lb
;
else if(!getfloat(&x
.val
[n
])) break;
if(!getfloat(&y
.val
[n
])) break; } }
if('0'<=c
&& c
<='9') continue;
if(!p
->lbf
&& p
->lb
>(p
->val
[i
])) p
->lb
= p
->val
[i
];
if(!p
->ubf
&& p
->ub
<(p
->val
[i
])) p
->ub
= p
->val
[i
]; }
x
.lbf
= x
.ubf
= y
.lbf
= y
.ubf
= 0;
again
: switch(argv
[0][0]) {
numb(&konst
,&argc
,&argv
);
if(!numb(&x
.lb
,&argc
,&argv
)) break;
if(!numb(&x
.ub
,&argc
,&argv
)) break;
fprintf(stderr
, "Bad agrument\n");
if(auta
&&!x
.lbf
) x
.lb
= 0;
diag
= (float *)malloc((unsigned)i
);
r
= (float *)malloc((unsigned)i
);
if(r
==NULL
||!spline()) for(i
=0;i
<n
;i
++){
printf("%f\n",y
.val
[i
]); }
if(!('0'<=c
&&c
<='9' || c
=='-' || c
== '.' )) return(0);