static char *sccsid
= "@(#)spline.c 4.1 (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 h\e9i\e8=x\e9i\e8-x\e9i-1\e\e9\e8\e8
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,
h\e9i\e8y"\b\e9i-1\e9\e8\e8+2(h\e9i\e8+h\e9i+1\e8)y"\b\e9i\e8+h\e9i+1\e8y"\b\e9i+1\e8
= 6[(y\e9i+1\e8-y\e9i\e8)/h\e9i+1\e8-(y\e9i\e8-y\e9i-1\e8)/h\e9i\e8] i=1,2,...,n
where y"\b\e90\e8 = y"\b\e9n+1\e8 = 0
This is a symmetric tridiagonal system of the form
| a\e91\e8 h\e92\e8 | |y"\b\e91\e8| |b\e91\e8|
| h\e92\e8 a\e92\e8 h\e93\e8 | |y"\b\e92\e8| |b\e92\e8|
| h\e93\e8 a\e93\e8 h\e94\e8 | |y"\b\e93\e8| = |b\e93\e8|
It can be triangularized into
| d\e91\e8 h\e92\e8 | |y"\b\e91\e8| |r\e91\e8|
| d\e92\e8 h\e93\e8 | |y"\b\e92\e8| |r\e92\e8|
| d\e93\e8 h\e94\e8 | |y"\b\e93\e8| = |r\e93\e8|
d\e9i\e8 = a\e9i\e8 - h\e9i\e8\b\e82\e9/d\e9i-1\e8 1<i<\b_n
r\e9i\e8 = b\e9i\e8 - h\e9i\e8r\e9i-1\e8/d\e9i-1\ei\e8 1<\b_i<\b_n
y"\b\e9n\e8 = r\e9n\e8/d\e9n\e8
y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8y"\b\e9i+1\e8)/d\e9i\e8 1<\b_i<n
superficially, d\e9i\e8 and r\e9i\e8 don't have to be stored for they can be
recalculated backward by the formulas
d\e9i-1\e8 = h\e9i\e8\b\e82\e9/(a\e9i\e8-d\e9i\e8) 1<i<\b_n
r\e9i-1\e8 = (b\e9i\e8-r\e9i\e8)d\e9i-1\e8/h\e9i\e8 1<i<\b_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\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b_________________________
The boundary conditions are easily generalized to handle
y\e90\e8\b" = ky\e91\e8\b", y\e9n+1\e8\b\b\b" = ky\e9n\e8\b"
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 h\e91\e8 to a\e91\e8 and h\e9n+1\e8 to a\e9n\e8.
Periodic case\b\b\b\b\b\b\b\b\b\b\b\b\b_____________
To do this, add 1 more row and column thus
| a\e91\e8 h\e92\e8 h\e91\e8 | |y\e91\e8\b"| |b\e91\e8|
| h\e92\e8 a\e92\e8 h\e93\e8 | |y\e92\e8\b"| |b\e92\e8|
| h\e93\e8 a\e94\e8 h\e94\e8 | |y\e93\e8\b"| |b\e93\e8|
| h\e91\e8 h\e90\e8 a\e90\e8 | | .| | .|
where h\e90\e8=\b_ h\e9n+1\e8
The same diagonalization procedure works, except for
the effect of the 2 corner elements. Let s\e9i\e8 be the part
of the last element in the i\e8th\e9 "diagonalized" row that
arises from the extra top corner element.
s\e9i\e8 = -s\e9i-1\e8h\e9i\e8/d\e9i-1\e8 2<\b_i<\b_n+1
After "diagonalizing", the lower corner element remains.
Call t\e9i\e8 the bottom element that appears in the i\e8th\e9 colomn
as the bottom element to its left is eliminated
t\e9i\e8 = -t\e9i-1\e8h\e9i\e8/d\e9i-1\e8
Evidently t\e9i\e8 = s\e9i\e8.
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.
u\e9i\e8 = u\e9i-1\e8-s\e9i-1\e8*t\e9i-1\e8/d\e9i-1\e8
v\e9i\e8 = v\e9i-1\e8-r\e9i-1\e8*t\e9i-1\e8/d\e9i-1\e8 2<\b_i<\b_n+1
The back solution is now obtained as follows
y"\b\e9n+1\e8 = (r\e9n+1\e8+v\e9n+1\e8)/(d\e9n+1\e8+s\e9n+1\e8+t\e9n+1\e8+u\e9n+1\e8)
y"\b\e9i\e8 = (r\e9i\e8-h\e9i+1\e8*y\e9i+1\e8-s\e9i\e8*y\e9n+1\e8)/d\e9i\e8 1<\b_i<\b_n
Interpolation in the interval x\e9i\e8<\b_x<\b_x\e9i+1\e8 is by the formula
y = y\e9i\e8x\e9+\e8 + y\e9i+1\e8x\e9-\e8 -(h\e82\e9\b\e9i+1\e8/6)[y"\b\e9i\e8(x\e9+\e8-x\e9+\e8\e8\b3\e9)+y"\b\e9i+1\e8(x\e9-\e8-x\e9-\e8\b\e83\e9)]
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);