procedure QUINAT(N1,N2,X,Y,B,C,D,E,F);
value N1,N2;
integer N1,N2;
real array X,Y,B,C,D,E,F;
comment QUINAT COMPUTES THE COEFFICIENTS OF A QUINTIC NATURAL SPLINE
S(X) INTERPOLATING THE ORDINATES Y[I] AT POINTS X[I], I = N1
THROUGH N2. FOR XX IN (X[I],X[I+1]) THE VALUE OF THE SPLINE
FUNCTION S(XX) IS GIVEN BY THE FIFTH DEGREE POLYNOMIAL:
S(XX) = ((((F[I]*T+E[I])*T+D[I])*T+C[I])*T+B[I])*T+Y[I]
WITH T = XX - X[I].
INPUT:
N1,N2 SUBSCRIPT OF FIRST AND LAST DATA POINT RESPECTIVELY,
IT IS REQUIRED THAT N2 > N1 + 1,
X,Y[N1::N2] ARRAYS WITH X[I] AS ABSCISSA AND Y[I] AS ORDINATE
OF THE I-TH DATA POINT. THE ELEMENTS OF THE ARRAY X
MUST BE STRICTLY MONOTONE INCREASING (BUT SEE BELOW FOR
EXCEPTIONS TO THIS).
OUTPUT:
B,C,D,E,F[N1::N2] ARRAYS COLLECTING THE COEFFICIENTS OF THE
QUINTIC NATURAL SPLINE S(XX) AS DESCRIBED ABOVE.
SPECIFICALLY B[I] = S'(X[I]), C[I] = S"(X[I])/2,
D[I] = S"'(X[I])/6, E[I] = S""(X[I])/24,
F[I] = S""'(X[I]+0)/120. F[N2] IS NEITHER USED OR
ALTERED. THE ARRAYS B,C,D,E,F MUST ALWAYS BE DISTINCT.
OPTIONS:
1. THE REQUIREMENT THAT THE ELEMENTS OF THE ARRAY X BE
STRICTLY MONOTONE INCREASING CAN BE RELAXED TO ALLOW TWO
OR THREE CONSECUTIVE ABSCISSAS TO BE EQUAL AND THEN
SPECIFYING VALUES OF THE FIRST AND SECOND DERIVATIVES OF
THE SPLINE FUNCTION AT SOME OF THE INTERPOLATING POINTS.
SPECIFICALLY
IF X[J] = X[J+1] THEN S(X[J]) = Y[J] AND S'(X[J]) = Y[J+1],
IF X[J] = X[J+1] = X[J+2] THEN IN ADDITION S"(X[J]) =Y[J+2].
NOTE THAT S""(X) IS DISCONTINUOUS AT A DOUBLE KNOT AND IN
ADDITION S"'(X) IS DISCONTINUOUS AT A TRIPLE KNOT. AT A
DOUBLE KNOT, X[J] = X[J+1], THE OUTPUT COEFFICIENTS HAVE THE
FOLLOWING VALUES:
B[J] = S'(X[J]) = B[J+1]
C[J] = S"(X[J])/2 = C[J+1]
D[J] = S"'(X[J])/6 = D[J+1]
E[J] = S""(X[J]-0)/24 E[J+1] = S""(X[J]+0)/24
F[J] = S""'(X[J]-0)/120 F[J+1] = S""'(X[J]+0)/120
THE REPRESENTATION OF S(XX) REMAINS VALID IN ALL INTERVALS
PROVIDED THE REDEFINITION Y[J+1] := Y[J] IS MADE
IMMEDIATELY AFTER THE CALL OF THE PROCEDURE QUINAT. AT A
TRIPLE KNOT, X[J] = X[J+1] = X[J+2], THE OUTPUT COEFFICIENTS
HAVE THE FOLLOWING VALUES:
B[J] = S'(X[J]) = B[J+1] = B[J+2]
C[J] = S"(X[J])/2 = C[J+1] = C[J+2]
D[J] = S"'(X[J]-0)/6 D[J+1] = 0 D[J+2] = S"'(X[J]+0)/6
E[J] = S""(X[J]-0)/24 E[J+1] = 0 E[J+2] = S""(X[J]+0)/24
F[J] = S""'(X[J]-0)/120 F[J+1]=0 F[J+2]=S""'(X[J]+0)/120
THE REPRESENTATION OF S(XX) REMAINS VALID IN ALL INTERVALS
PROVIDED THE REDEFINITION Y[J+2] := Y[J+1] := Y[J] IS MADE
IMMEDIATELY AFTER THE CALL OF THE PROCEDURE QUINAT.
2. THE ARRAY X MAY BE MONOTONE DECREASING INSTEAD OF
INCREASING;
if N2 > N1 + 1 then
begin
integer M;
real B1,P,PQ,PQQR,PR,P2,P3,Q,QR,Q2,Q3,R,R2,S,T,U,V;
comment COEFFICIENTS OF A POSITIVE DEFINITE, PENTADIAGONAL
MATRIX STORED IN D,E,F[N1+1::N2-2];
M ≔ 2-2;
Q ≔ [N1+1]-X[N1];
R ≔ [N1+2]-X[N1+1];
Q2 ≔ ×Q;
R2 ≔ ×R;
QR ≔ +R;
E[N1] ≔ 0·0;
D[N1] ≔ E[N1];
D[N1+1] ≔ ̲f Q=0·0 then 0·0 else 6·0×Q×Q2/(QR×QR);
for I ≔ 1+1 step 1 until M do
begin
P ≔ ; Q ≔ ; R ≔ [I+2]-X[I+1];
P2 ≔ 2; Q2 ≔ 2; R2 ≔ ×R; PQ ≔ R; QR ≔ +R;
if Q=0·0 then D[I+1] ≔ [I] ≔ [I-1] ≔ ·0 else
begin
Q3 ≔ 2×Q; PR ≔ ×R; PQQR ≔ Q×QR;
D[I+1] ≔ ·0×Q3/(QR×QR);
D[I] ≔ [I]+(Q+Q)×(15·0×PR×PR+(P+R)×Q×(20·0×PR+7·0×Q2)
+Q2×(8·0×(P2+R2)+21·0×PR+Q2+Q2))/(PQQR×PQQR);
D[I-1] ≔ [I-1]+6·0×Q3/(PQ×PQ);
E[I] ≔ 2×(P×QR+3·0×PQ×(QR+R+R))/(PQQR×QR);
E[I-1] ≔ [I-1]+Q2×(R×PQ+3·0×QR×(PQ+P+P))/(PQQR×PQ);
F[I-1] ≔ 3/PQQR;
end
end;
if ABS(R) > 0·0 then D[M] ≔ [M]+6·0×R×R2/(QR×QR);
comment FIRST AND SECOND ORDER DIVIDED DIFFERENCES OF THE GIVEN
FUNCTION VALUES,STORED IN B[N1+1::N2] AND C[N1+2::N2]
RESPECTIVELY, TAKE CARE OF DOUBLE AND TRIPLE KNOTS;
S ≔ [N1];
for I ≔ 1+1 step 1 until N2 do
if X[I]=X[I-1] then B[I] ≔ [I] else
begin
B[I] ≔ Y[I]-S)/(X[I]-X[I-1]);
S ≔ [I];
end;
comment f̲o̲r̲ I=N=N1+2 s̲t̲e̲p̲ 1 u̲n̲t̲i̲l̲ N2 d̲o̲;
for I ≔ 1+2 step 1 until N2 do
ifX[I]=X[I-2] then
begin C[I] ≔ [I]×0·5; B[I] ≔ [I-1] end
else C[I] ≔ B[I]-B[I-1])/(X[I]-X[I-2]);
comment SOLVE THE LINEAR SYSTEM WITH C[I+2]-C[I+1]
AS RIGHT-HAND SIDE;
if M > N1 then
begin
P ≔ [N1] ≔ [M] ≔ [N1] ≔ [M-1] ≔ [M] ≔ ·0;
C[N1+1] ≔ [N1+3]-C[N1+2]; D[N1+1] ≔ ·0/D[N1+1];
end;
for I ≔ 1+2 step 1 until M do
begin
Q ≔ [I-1]×E[I-1];
D[I] ≔ ·0/(D[I]-P×F[I-2]-Q×E[I-1]);
E[I] ≔ [I]-Q×F[I-1];
C[I] ≔ [I+2]-C[I+1]-P×C[I-2]-Q×C[I-1];
P ≔ [I-1]×F[I-1];
end;
M ≔ 1+1; C[N2-1] ≔ [N2] ≔ ·0;
for I ≔ 2-2 step -1 until M do
C[I] ≔ C[I]-E[I]×C[I+1]-F[I]×C[I+2])×D[I];
comment INTEGRATE THE THIRD DERIVATIVE OF S(X);
M ≔ 2-1;
Q ≔ [N1+1]-X[N1]; R ≔ [N1+2]-X[N1+1]; B1 ≔ N1+1);
Q3 ≔ ×Q×Q; QR ≔ +R;
V ≔ ≔ ̲f QR=0·0 then 0·0 else C[N1+1]/QR;
F[N1] ≔ ̲f Q=0·0 then 0·0 else V/Q;
for I ≔ 1+1 step 1 until M do
begin
P ≔ ; Q ≔ ;
R ≔ ̲f I=N2-1 then 0·0 else X[I+2]-X[I+1];
P3 ≔ 3; Q3 ≔ ×Q×Q; PQ ≔ R; QR ≔ +R;
S ≔ ; T ≔ ̲f QR=0·0 then 0·0 else (C[I+1]-C[I])/QR;
U ≔ ; V ≔ -S;
if PQ=0·0 then
begin C[I] ≔ ·5×Y[I+1]; D[I] ≔ [I] ≔ [I] ≔ ·0 end
else
begin
F[I] ≔ ̲f Q=0·0 then F[I-1] else V/Q;
E[I] ≔ ·0×S;
D[I] ≔ 0·0×(C[I]-Q×S);
C[I] ≔ [I]×(P-Q)+(B[I+1]-B[I]+(U-E[I])×P3
-(V+E[I])×Q3)/PQ;
B[I] ≔ P×(B[I+1]-V×Q3)+Q×(B[I]-U×P3))/PQ
-P×Q×(D[I]+E[I]×(Q-P));
end
end I;
comment END POINTS X[N1] AND X[N2];
P ≔ [N1+1]-X[N1]; S ≔ [N1]×P×P×P;
E[N1] ≔ [N1] ≔ ·0;
C[N1] ≔ [N1+1]-10·0×S;
C[N1] ≔ 1-(C[N1]+S)×P;
Q ≔ [N2]-X[N2-1]; T ≔ [N2-1]×Q×Q×Q;
E[N2] ≔ [N2] ≔ ·0;
C[N2] ≔ [N2-1]+10·0×T;
B[N2] ≔ [N2]+(C[N2]-T)×Q;
end QUINAT;
procedure QUINEQ(N1,N2, Y,B,C,D,E,F);
value N1,N2;
integer N1,N2;
real array Y,B,C,D,E,F;
comment QUINEQ COMPUTES THE COEFFICIENTS OF A QUINTIC NATURAL SPLINE
S(X) INTERPOLATING THE ORDINATES Y[I] AT EQUIDISTANT POINTS X[I],
I = N1 THROUGH N2. FOR XX IN (X[I],X[I+1]) THE VALUE OF THE
SPLINE FUNCTION S(XX) IS GIVEN BY THE FIFTH DEGREE POLYNOMIAL:
S(XX) = ((((F[I]*T+E[I])*T+D[I])*T+C[I])*T+B[I])*T+Y[I]
WITH T = (XX - X[I])/(X[I+1] - X[I]).
INPUT:
N1, N2 SUBSCRIPT OF FIRST AND LAST DATA POINT RESPECTIVELY,
IT IS REQUIRED THAT N2 > N1 + 1,
Y[N1::N2] THE GIVEN FUNCTION VALUES (ORDINATES).
OUTPUT:
B,C,D,E,F[N1::N2] ARRAYS COLLECTING THE COEFFICIENTS OF THE
QUINTIC NATURAL SPLINE S(XX) AS DESCRIBED ABOVE.
SPECIFICALLY B[I] = S'(X[I]), C[I] = S"(X[I])/2,
D[I] = S"'(X[I])/6, E[I] = S""(X[I])/24,
F[I] = S""'(X[I]+0)/120. F[N2] IS NEITHER USED
NOR ALTERED. THE ARRAYS Y,B,C,D MUST ALWAYS BE
DISTINCT. IF E AND F ARE NOT WANTED, THE CALL
QUINEQ(N1,N2,Y,B,C,D,D,D) MAY BE USED TO SAVE STORAGE
LOCATIONS;
if N2>N1+1 then
begin
integer N;
real P,Q,R,S,T,U,V;
N ≔ 2-3; P ≔ ≔ ≔ ≔ ≔ ·0;
for I ≔ 1 step 1 until N do
begin
U ≔ ×R; B[I] ≔ ·0/(66·0-U×R-Q);
C[I] ≔ ≔ 6·0-U;
D[I] ≔ [I+3]-3·0×(Y[I+2]-Y[I+1])-Y[I]-U×S-Q×T;
Q ≔ ; P ≔ [I]; T ≔ ; S ≔ [I]
end I;
D[N+1] ≔ [N+2] ≔ ·0;
for I ≔ step -1 until N1 do
D[I] ≔ D[I]-C[I]×D[I+1]-D[I+2])×B[I];
N ≔ 2-1; Q ≔ ·0; R ≔ ≔ ≔ [N1];
for I ≔ 1+1 step 1 until N do
begin
P ≔ ; Q ≔ ; R ≔ [I]; S ≔ ;
F[I] ≔ ≔ -Q-Q+R;
E[I] ≔ ≔ ·0×(-P+Q);
D[I] ≔ 0·0×(P+Q);
C[I] ≔ ·5×(Y[I+1]+Y[I-1]+S-T)-Y[I]-U;
B[I] ≔ ·5×(Y[I+1]-Y[I-1]-S-T)-D[I]
end I;
F[N1] ≔ ; E[N1] ≔ [N2] ≔ [N1] ≔ [N2] ≔ ·0;
C[N1] ≔ [N1+1]-10·0×V; C[N2] ≔ [N2-1]+10·0×T;
B[N1] ≔ [N1+1]-Y[N1]-C[N1]-V; B[N2] ≔ [N2]-Y[N2-1]+C[N2]-T
end QUINEQ;
procedure QUINDF(N1,N2, X,Y,B,C,D,E,F);
value N1,N2;
integer N1,N2;
real array X,Y,B,C,D,E,F;
comment QUINDF COMPUTES THE COEFFICIENTS OF A QUINTIC NATURAL SPLINE
S(X) FOR WHICH THE ORDINATES Y[I] AND THE FIRST DERIVATIVES B[I]
ARE SPECIFIED AT POINTS X[I], I = N1 THROUGH N2. FOR XX IN
(X[I],X[I+1]) THE VALUE OF THE SPLINE FUNCTION S(XX) IS GIVEN
BY THE FIFTH DEGREE POLYNOMIAL:
S(XX) = (((F[I]*T+E[I])*T+D[I])*T+C[I])*T+B[I])*T+Y[I]
WITH T = XX - X[I].
INPUT:
N1, N2 SUBSCRIPT OF FIRST AND LAST DATA POINT RESPECTIVELY,
IT IS REQUIRED THAT N2 > N1,
X,Y,B[N1::N2] ARRAYS WITH X[I] AS ABSCISSA, Y[I] AS ORDINATE
AND B[I] AS FIRST DERIVATIVE AT THE I-TH DATA POINT.
THE ELEMENTS OF THE ARRAY X MUST BE STRICTLY MONOTONE
INCREASING OR DECREASING.
OUTPUT:
C,D,E,F[N1::N2] ARRAYS COLLECTING THE COEFFICIENTS OF THE
QUINTIC NATURAL SPLINE S(XX) AS DESCRIBED ABOVE.
E[N2] AND F[N2] ARE NEITHER USED NOR ALTERED. THE
ARRAYS C,D,E,F MUST ALWAYS BE DISTINCT;
if N2 > N1 then
begin
integer M2;
real CC,G,H,HH,H2,HH2,P,PP,Q,QQ,R,RR;
M2 ≔ 2-1; CC ≔ H ≔ P ≔ Q ≔ R ≔ ≔ ·0;
for I ≔ 1 step 1 until M2 do
begin
H ≔ ·0/(X[I+1]-X[I]); H2 ≔ ×H; D[I] ≔ ·0×(HH+H) - G×HH;
P ≔ Y[I+1]-Y[I])×H×H2; Q ≔ B[I+1]+B[I])×H2;
R ≔ B[I+1]-B[I])×H2;
C[I] ≔ C ≔ 0·0×(P-PP) - 5·0×(Q-QQ) + R + RR + G×CC;
G ≔ /D[I]; HH ≔ ; HH2 ≔ 2; PP ≔ ; QQ ≔ ; RR ≔
end I;
C[N2] ≔ -10·0×PP + 5·0×QQ + RR + G×CC)/(3·0×HH - G×HH);
for I ≔ 2 step -1 until N1 do
begin
D[I+1] ≔ ·0/(X[I+1]-X[I]); C[I] ≔ C[I] + C[I+1]×D[I+1])/D[I]
end I;
for I ≔ 1 step 1 until M2 do
begin
H ≔ [I+1]; H2 ≔ ×H;
P ≔ Y[I+1]-Y[I])×H×H2 - B[I]×H2 - C[I]×H;
Q ≔ B[I+1]-B[I])×H2 - C[I]×(H+H);
R ≔ C[I+1]-C[I])×H;
G ≔ - 3·0×P; RR ≔ - 3·0×(P+G); QQ ≔ -RR - RR + G;
F[I] ≔ R×H2; E[I] ≔ Q×H; D[I] ≔ -RR - QQ + P
end I;
D[N2] ≔ ·0
end QUINDF;