'******************************************************** '* CUBIC SPLINE INTERPOLATION OF A DISCRETE FUNCTION * '* F(X), GIVEN BY N POINTS X(I), Y(I) * '* ---------------------------------------------------- * '* SAMPLE RUN: * '* * '* (Interpolate F(X), defined by 11 points: * '* 0.0 0.01 * '* 1.0 1.01 * '* 2.0 3.99 * '* 3.0 8.85 * '* 4.0 15.0 * '* 5.0 25.1 * '* 6.0 37.0 * '* 7.0 50.0 * '* 8.0 63.9 * '* 9.0 81.2 * '* 10.0 100.5 for X = 8.7 ) * '* * '* CUBIC SPLINE INTERPOLATION: * '* * '* For X= 8.70000000 Y= 75.67009329 * '* * '* ---------------------------------------------------- * '* Ref.: From Numath Library By Tuan Dang Trong in * '* Fortran 77. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '******************************************************** DEFDBL A-H, O-Z DEFINT I-N N = 11 'Number of given points DIM X(N), Y(N) 'Coordinates of N points DIM B(N), C(N), D(N) 'Coefficients of cubic spline ' U,V as Double Coordinates of interpolated point X(1) = 0#: X(2) = 1#: X(3) = 2#: X(4) = 3#: X(5) = 4#: X(6) = 5# X(7) = 6#: X(8) = 7#: X(9) = 8#: X(10) = 9#: X(11) = 10# Y(1) = .01#: Y(2) = 1.01#: Y(3) = 3.99: Y(4) = 8.85: Y(5) = 15#: Y(6) = 25.1: Y(7) = 37#: Y(8) = 50#: Y(9) = 63.9: Y(10) = 81.2: Y(11) = 100.5 GOSUB 2000 'call SPLINE (N,X,Y,B,C,D) U = 8.7 GOSUB 1000 'V = SEVAL(N,U,X,Y,B,C,D) F$ = " For X=####.######## Y=####.########" CLS PRINT PRINT " CUBIC SPLINE INTERPOLATION:" PRINT USING F$; U; V PRINT END 'of main 1000 'SEVAL(N, U,X, Y, B, C, D) '----------------------------------------------------------------------- ' EVALUATE A CUBIC SPLINE INTERPOLATION OF A DISCRETE FUNCTION F(X), ' GIVEN IN N POINTS X(I), Y(I). THE B, C AND D COEFFICIENTS DEFINING ' THE BEST CUBIC SPLINE FOR THE GIVEN POINTS, ARE CALCULATED BEFORE ' BY THE SPLINE SUBROUTINE. ' ' INPUTS: ' N NUMBER OF POINTS OF CURVE Y = F(X) ' U ABSCISSA OF POINT TO BE INTERPOLATED ' X,Y TABLES OF DIMENSION N, STORING THE COORDINATES OF CURVE F(X) ' B,C,D TABLES STORING THE COEFFICIENTS DEFINING THE CUBIC SPLINE ' ' OUTPUT: ' V INTERPOLATED VALUE ' = Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I))) ' WITH DX = U-X(I), U BETWEEN X(I) AND X(I+1) ' ' REFERENCE : ' FORSYTHE,G.E. (1977) COMPUTER METHODS FOR MATHEMATICAL ' COMPUTATIONS. PRENTICE-HALL,INC. '------------------------------------------------------------------------ ' Labels: 10, 20, 30 I = 1 ' BINARY SEARCH IF (I >= N) THEN I = 1 IF (U < X(I)) THEN GOTO 10 IF (U <= X(I + 1)) THEN GOTO 30 10 I = 1 J = N + 1 20 K = (I + J) / 2 IF (U < X(K)) THEN J = K IF (U >= X(K)) THEN I = K IF (J > I + 1) THEN GOTO 20 ' SPLINE EVALUATION 30 DX = U - X(I) V = Y(I) + DX * (B(I) + DX * (C(I) + DX * D(I))) RETURN 2000 'SPLINE (N, X, Y, B, C, D) '--------------------------------------------------------------------- ' THIS SUBROUTINE CALCULATES THE COEFFICIENTS B,C,D OF A CUBIC ' SPLINE TO BEST APPROXIMATE A discrete FONCTION GIVEN BY N POINTS ' ' INPUTS: ' N NUMBER OF GIVEN POINTS ' X,Y VECTORS OF DIMENSION N, STORING THE COORDINATES ' OF FUNCTION F(X) ' ' OUTPUTS: ' A,B,C VECTORS OF DIMENSION N, STORING THE COEFFICIENTS ' OF THE CUBIC SPLINE ' ' REFERENCE: ' FORSYTHE,G.E. (1977) COMPUTER METHODS FOR MATHEMATICAL ' COMPUTATIONS. PRENTICE-HALL,INC. '--------------------------------------------------------------------- ' Labels: 15, 50 NM1 = N - 1 IF (N < 2) THEN RETURN IF (N < 3) THEN GOTO 50 ' BUILD THE TRIDIAGONAL SYSTEM ' B (DIAGONAL), D (UPPERDIAGONAL) , C (SECOND MEMBER) D(1) = X(2) - X(1) C(2) = (Y(2) - Y(1)) / D(1) FOR I = 2 TO N - 1 D(I) = X(I + 1) - X(I) B(I) = 2# * (D(I - 1) + D(I)) C(I + 1) = (Y(I + 1) - Y(I)) / D(I) C(I) = C(I + 1) - C(I) NEXT I ' CONDITIONS AT LIMITS ' THIRD DERIVATIVES OBTAINED BY DIVIDED DIFFERENCES B(1) = -D(1) B(N) = -D(N - 1) C(1) = 0# C(N) = 0# IF (N = 3) THEN GOTO 15 C(1) = C(3) / (X(4) - X(2)) - C(2) / (X(3) - X(1)) C(N) = C(N - 1) / (X(N) - X(N - 2)) - C(N - 2) / (X(N - 1) - X(N - 3)) C(1) = C(1) * D(1) * D(1) / (X(4) - X(1)) C(N) = -C(N) * D(N - 1) * D(N - 1) / (X(N) - X(N - 3)) ' FORWARD ELIMINATION 15 FOR I = 2 TO N - 1 T = D(I - 1) / B(I - 1) B(I) = B(I) - T * D(I - 1) C(I) = C(I) - T * C(I - 1) NEXT I ' BACK SUBSTITUTION C(N) = C(N) / B(N) FOR L = 1 TO N - 1 I = N - L C(I) = (C(I) - D(I) * C(I + 1)) / B(I) NEXT L ' COEFFICIENTS OF 3RD DEGREE POLYNOMIAL B(N) = (Y(N) - Y(NM1)) / D(NM1) + D(NM1) * (C(NM1) + 2# * C(N)) FOR I = 1 TO N - 1 B(I) = (Y(I + 1) - Y(I)) / D(I) - D(I) * (C(I + 1) + 2# * C(I)) D(I) = (C(I + 1) - C(I)) / D(I) C(I) = 3# * C(I) NEXT I C(N) = 3# * C(N) D(N) = D(NM1) RETURN ' CASE N = 2 50 B(1) = (Y(2) - Y(1)) / (X(2) - X(1)) C(1) = 0# D(1) = 0# B(2) = B(1) C(2) = 0# D(2) = 0# RETURN ' end of file tseval.bas