! Last change: DOS 29 Jul 2000 12:22 pm ! *** copyright 2000 *** ! *** filename splint.f95 *** John F. Monahan ** ! ********************** program splint ! splint -- demonstration of cubic spline interpolation ! implicit none integer, parameter :: n = 7 integer, parameter :: npts = 100 real, dimension(n) :: xk,z,mn,md,h real xi,s,fn,fd real splev integer i ! 21 format(4f12.6) 24 format(f9.6,2x,4f10.6) ! open( unit=6, file='splint.out') open( unit=9, file='splint.dat') ! get n points of curve do i = 1,n ! x's equally spaced on (-4,4) xi = 4.*real(2*i-1-n)/real(n-1) xk(i) = xi ! function interpolated is simple-looking z(i) = 1./(1.+xi*xi) end do ! loop on i ! find natural spline first call splstn(xk,z,mn,h,n) WRITE(6,21) mn ! now find spline with derivative conditions md(1) = -8./289. md(n) = 8./289. ! call splstd(xk,z,md,h,n) WRITE(6,21) md ! ! write out results for 2 interpolants do i = 1,npts xi = 4.*real(2*i-1-npts)/real(npts-1) s = 1./(1.+xi*xi) fn = splev(xi,xk,z,mn,h,n) fd = splev(xi,xk,z,md,h,n) write(9,24) xi,s,fn,fd end do ! loop on i stop end program splint SUBROUTINE SPLSTN(XK,Y,M,H,N) ! NATURAL CUBIC SPLINES WITH UNEQUALLY SPACED KNOTS ! M(J) GIVES SECOND DERIV AT XKNOT(J) J=1,...,N ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( N ), INTENT(IN) :: XK,Y REAL, DIMENSION( N ), INTENT(OUT) :: M,H REAL, DIMENSION(3*N) :: W REAL S INTEGER J,NM1 ! NM1 = N-1 ! COMPUTE DIFFERENCES H(J) AND FIRST DIFS OF Y'S DO J = 2,N H(J) = XK(J)-XK(J-1) M(J) = (Y(J)-Y(J-1))/H(J) ! S(J) END DO ! LOOP ON J ! FILL MATRIX TRIDIAGONAL MATRIX A IN VECTOR W ! DIAGONAL W(J) = A(J,J) = 2 ! SUBDIAGONAL W(N+J) = A(J,J-1) = LAMBDA(J) ! SUPERDIAGONAL W(2*N+J) = A(J,J+1) = MU(J) ! ! M(J) NOW HOLDS FIRST DIFFERENCES ! ! FIX TWO ENDS DO J = 2,NM1 S = H(J)+H(J+1) W(2*N+J) = H(J+1)/S ! SUPERDIAGONAL W(N+J) = H(J)/S ! SUBDIAGONAL W(J) = 2. ! DIAGONAL M(J) = 6.*(M(J+1)-M(J))/S ! SECOND DIFFERENCES END DO ! LOOP ON J ! FINISH ENDS W(1) = 2. W(N) = 2. W(2*N+1) = 0. W(2*N) = 0. M(1) = 0. M(N) = 0. ! SOLVE TRIDIAGONAL SYSTEM CALL STRIDS(W,M,N) RETURN END SUBROUTINE SPLSTN SUBROUTINE SPLSTD(XK,Y,M,H,N) ! CUBIC SPLINES WITH UNEQUALLY SPACED KNOTS WITH DERIVATIVE CONDITION ! DERIVATIVES AT XK(1) & XK(N) ARE GIVEN IN M(1) & M(N) AS INPUT ! M(J) GIVES SECOND DERIV AT XKNOT(J) J=1,...,N ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( N ), INTENT(IN) :: XK,Y REAL, DIMENSION( N ), INTENT(OUT) :: M,H REAL, DIMENSION(3*N) :: W REAL S INTEGER J,NM1 ! NM1 = N-1 ! COMPUTE DIFFERENCES H(J) AND FIRST DIFS OF Y'S W(N) = M(N) ! HOLD THIS DERIV FOR LATER DO J = 2,N H(J) = XK(J)-XK(J-1) M(J) = (Y(J)-Y(J-1))/H(J) ! S(J) END DO ! LOOP ON J ! FILL MATRIX TRIDIAGONAL MATRIX A IN VECTOR W ! DIAGONAL W(J) = A(J,J) = 2 ! SUBDIAGONAL W(N+J) = A(J,J-1) = LAMBDA(J) ! SUPERDIAGONAL W(2*N+J) = A(J,J+1) = MU(J) ! ! M(J) NOW HOLDS FIRST DIFFERENCES ! ! FIX TWO ENDS M(1) = 6.*(M(2)-M(1))/H(2) ! SECOND DIFFERENCE DO J = 2,NM1 S = H(J)+H(J+1) W(2*N+J) = H(J+1)/S ! SUPERDIAGONAL W(N+J) = H(J)/S ! SUBDIAGONAL W(J) = 2. ! DIAGONAL M(J) = 6.*(M(J+1)-M(J))/S ! SECOND DIFFERENCES END DO ! LOOP ON J M(N) = 6.*(W(N)-M(N))/H(N) ! SECOND DIFFERENCE ** YES W(N) ! FINISH ENDS W(1) = 2. W(N) = 2. W(2*N+1) = 1. W(2*N) = 1. ! SOLVE TRIDIAGONAL SYSTEM CALL STRIDS(W,M,N) RETURN END SUBROUTINE SPLSTD SUBROUTINE STRIDS(W,X,N) ! SOLVES TRIDIAGONAL SYSTEM: A*(X-OUT)=(X-IN) ! WHERE DIAGONAL ELEMENTS IN W(I) = A(I,I) ! SUBDIAGONAL ELEMENTS IN W(N+I) = A(I,I-1) ! SUPERDIAGONAL ELEMENTS IN W(2*N+I) = A(I,I+1) ! ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( 3*N ), INTENT(IN OUT) :: W REAL, DIMENSION( N ), INTENT(IN OUT) :: X INTEGER I,K ! COMPUTE LU FACTORIZATION (L IS UNIT LOWER) DO K = 2,N W(N+K) = W(N+K)/W(K-1) W(K) = W(K) - W(N+K)*W(2*N+K-1) END DO ! LOOP ON K ! SOLVE LOWER TRIANGULAR (BIDIAGONAL) SYSTEM DO I = 2,N X(I) = X(I) - W(N+I)*X(I-1) END DO ! LOOP ON I ! SOLVE UPPER TRIANGULAR SYSTEM (START AT BOTTOM) X(N) = X(N) / W(N) DO K = 2,N I = N+1-K X(I) = ( X(I) - W(2*N+I)*X(I+1) ) / W(I) END DO ! LOOP ON K RETURN END SUBROUTINE STRIDS REAL FUNCTION SPLEV(X,XK,Y,M,H,N) ! EVALUATES SPLINE AT X XK(JC-1) LT X LE XK(JC) ! CUBIC SPLINES WITH UNEQUALLY SPACED KNOTS ! M(J) GIVES SECOND DERIV AT XKNOT(J) J=1,...,N ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(IN) :: X REAL, DIMENSION( N ), INTENT(IN) :: XK,Y,M,H REAL DA,DB INTEGER, SAVE :: J = 0 INTEGER JFIND ! INITIAL VALUE -- RETURN THIS IF OUT OF RANGE SPLEV = 0. ! IF X OUT OF RANGE, SPLEV=0 AND QUIT IF( (X .LT. XK(1)) .OR. (X .GT. XK(N)) ) RETURN ! IF LAST CALL IN SAME INTERVAL, NO NEED TO SEARCH IF( ( J .LE. 1 ) .OR. ( J .GT. N ) ) THEN J = JFIND(X,XK,N) ! FIND THE INTERVAL ELSE IF( (X .LE. XK(J-1)) .OR. (X .GT. XK(J)) ) THEN J = JFIND(X,XK,N) ! FIND THE INTERVAL ENDIF ENDIF ! COMPUTE DIFFERENCES FROM NEARBY KNOTS DB = X-XK(J-1) DA = XK(J)-X ! EVALUATE CUBIC SPLINE SPLEV = ( (M(J-1)*DA**3 + M(J)*DB**3)/6. & & + (Y(J-1)-M(J-1)*H(J)*H(J)/6.)*DA & & + (Y(J)-M(J)*H(J)*H(J)/6.)*DB ) / H(J) RETURN END FUNCTION SPLEV INTEGER FUNCTION JFIND(X,XK,N) ! FIND J SUCH THAT XK(J-1) LT X LE XK(J) ! JFIND=2 IF X EQ XK(1) JFIND=0 IF X IS OUT OF RANGE ! J F MONAHAN AUGUST 1981 DEPT OF STAT, NCSU, RALEIGH, NC 27650 ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, INTENT(IN) :: X REAL, DIMENSION( N ), INTENT(IN) :: XK INTEGER JU,JL,I JFIND = 0 ! IF X OUT OF RANGE, JFIND=0 AND QUIT IF( (X .LT. XK(1)) .OR. (X .GT. XK(N)) ) RETURN ! CHECK SPECIAL CASE IF( X .EQ. XK(1) ) THEN JFIND = 2 RETURN END IF ! INITIALIZE POINTERS JU = N JL = 1 DO WHILE( JU-JL .GT. 1 ) ! TEST THIS ONE I = ( JU + JL ) / 2 IF( X .LE. XK(I) ) THEN JU = I ELSE JL = I END IF END DO ! WHILE( JU-JL .GT. 1 ) JFIND = JU RETURN END FUNCTION JFIND ! *** end of filename splint.f95 *********************