[插值拟合] 关于差值方法

Program www_fcode_cn
  use SplineMod
  Implicit None
  Integer , parameter :: N = 10
  Integer , parameter :: M = 37
  Real :: x(N) , y(N) , xx(M) , yy(M)
  integer :: i
  Open( 12 , File = 'a.txt' )
  Open( 13 , File = 'b.txt' )
  Do i = 1 , N
    Read( 12 , * ) x(i) , y(i)
  End Do
  call SplineMod_Init( N, x , y )
  Do i = 1 , M
    xx(i) = (i-1) * 0.5
    yy(i) = SplineMod_Interp( N , x , y , xx(i) )
  End Do
  Do i = 1 , M
    Write( 13 , * ) xx(i) , yy(i)
  End Do
End Program www_fcode_cn

Module SplineMod
  Implicit None  
  Integer , private ::  i , N
  Real , private , Allocatable :: S2(:) , Dely(:) , S3(:)

Contains

  Subroutine SplineMod_Init( iN , x , y )
    Integer , Intent( IN )  :: iN
    Real , Intent( IN )      :: x( iN ) , y( iN )
    Real :: B( iN ) , C( iN ) , H( iN ) , H1( iN ) , Delsqy
    Integer jj , N1
    N  = iN
    N1 = N - 1
    Allocate( S2( N ) )
    Allocate( Dely( N ) )
    Allocate( S3( N ) )
    Do i = 1 , N1
      H( i ) = x( i+1 ) - x( i )
      if ( abs(H(i)) < 1./3600 ) H(i) = 1./3600
      Dely( i ) = ( y(i+1) - y(i) ) / H( i ) 
    End Do
    Do i = 2 , N1 
      H1( i ) = H( i-1 ) + H( i ) 
      B( i )  = 0.5 * H(i-1) / H1( i ) 
      Delsqy = ( Dely(i) - Dely(i-1) ) / H1( i ) 
      S2( i ) = 2.0 * Delsqy
      C( i ) = 3.0 * Delsqy
    End Do
    S2( 1 ) = 0.0 
    S2( N ) = 0.0 
    Do jj = 1 , 26 
      Do i = 2 , N1 
        S2(i) = (C(i)-B(i)*S2(i-1)-(0.5-B(i))*S2(i+1)-S2(i))*1.0717968+S2(i) 
      End Do
    End Do
    Do i = 1 , N1 
      S3( i ) = ( S2(i+1) - S2(i) ) / H( i ) 
    End Do
  End Subroutine SplineMod_Init

  Real Function SplineMod_Interp( iN , x , y , T )
    Integer , Intent( IN ) :: iN
    Real , Intent( IN ) :: T
    Real , Intent( IN )      :: x( iN ) , y( iN )
    Integer i
    Real                 :: ht1 , ht2 , Delsqs    
    i = 1 
    if( ( T - x(i) ) <= 0.0 ) goto 17 
    if( ( T - x(N) ) <  0.0 )  goto 57 
    goto 59 
 56 if( ( T - x(i) ) <  0.0 ) goto 60 
    if( ( T - x(i) ) == 0.0 ) goto 17 
 57 i = i + 1 
    GOTO 56 
 59 i = N 
 60 i = i - 1 
 17 HT1 = T - x(i)
    HT2 = T - x(i+1)
    Delsqs = ( 2.0 * S2(i) + S2(i+1) + HT1 * S3(i) ) / 6.0 
    SplineMod_Interp = y(i) + HT1 * Dely( i ) + HT1 * HT2 * Delsqs 
  End Function SplineMod_Interp

  Subroutine SplineMod_UnInit()
    DeAllocate( S2 )
    DeAllocate( Dely )
    DeAllocate( S3 )
  End Subroutine SplineMod_UnInit
  
End Module SplineMod