fortran三次样条插值程序实例

Program testspline
    implicit none
    integer :: i 
    integer, parameter :: n = 50, m = 100, dp = 8
    real(dp)           :: x(n), y(n), sx(m), sy(m), dsy(m)
    real(dp)           :: x1 = -1.d0, x2 = 1.d0, pi = acos(-1.d0)
    
    !..==================插值基节点==================
    !..                    x,  y
    !..=============================================
    
    !..=============待插值节点以及一阶导数============
    !..                  sx, sy, sdy
    !..==============================================
    open( 100, file = 'origin.dat' )
    do i = 1, n
        x(i) = x1 + dble(i-1)*(x2-x1)/dble(n-1)
        y(i) = exp( x(i) ) * sin(pi*x(i))   !.. y = e^x + sin(pi*x)
        write( 100,* ) x(i), y(i), exp( x(i) ) * sin(pi*x(i)) + exp( x(i) ) * pi * cos(pi*x(i))  !.. y' = e^x + sin(pi*x) + e^x * pi * cos(pi*x)
    end do
    close( 100 )
    
    sy = 0.d0; dsy = 0.d0
    open( 100, file = 'spline.dat' )
    do i = 1, m
        sx(i) = x1 + dble(i-1)*(x2-x1)/dble(m-1)
        call spline( x, y, n, sx(i), sy(i), dsy(i), 1, n-1, n-2 )
        write( 100,* ) sx(i), sy(i), dsy(i)
    end do
    close( 100 )
    
End program testspline
    

subroutine spline( x, y, n, sx, f, f1, m, n1, n2 )  !.. 三次样条插值程序
    implicit none
    integer            :: i, j, k, m, n, n1, n2
    integer, parameter :: dp = 8
    Real(dp)           :: x(n), y(n), sx(m), f(m), f1(m)
    Real(dp)           :: s2(n), h(n1), dy(n1), s(n1), e(n2)
    Real(dp)           :: z, h1, h2, h3, h4
  
    do i = 1, n1
        h(i)  = x(i+1) - x(i)
        dy(i) = ( y(i+1) - y(i) ) / h(i)
    end do
    
    s2(1) = 0.d0; s2(n) = 0.d0
    do i = 2, n1
        s2(i) = 6.d0 * ( dy(i) - dy(i-1) )
    end do
    
    z = 0.5d0 / ( h(1) + h(2) )
    s(1) = -h(2) * z
    e(1) = s2(2) * z
    do i = 2, n2
        k    = i - 1
        j    = i + 1
        z    = 1.d0 / ( 2.d0*( h(i)+h(j) ) + h(i)*s(k) )
        s(i) = -h(j) * z
        e(i) = ( s2(j)-h(i)*e(k) ) * z
    end do
    
    s2(n1) = e(n2)
    do i = n2, 2, -1
        k     = i - 1
        s2(i) = s(k)*s2(i+1) + e(k)
    end do
    
    do i = 1, n1
        s(i) = ( s2(i+1) - s2(i) ) / h(i)
    end do
    
    i = 2
    k = 1
    do j = 1, m
        do 
            if ( sx(j) > x(i) ) then
                k = i
                i = i + 1
            else
                exit
            end if
        end do
        h1    = sx(j) - x(k)
        h2    = sx(j) - x(i)
        h3    = h1 * h2
        h4    = s2(k) + h1*s(k)
        z     = ( s2(i) + s2(k) + h4 ) / 6.d0
        f(j)  = y(k) + h1*dy(k) + h3*Z
        f1(j) = dy(k) + z*( h1+h2 ) + h3 * s(k) / 6.d0
    end do
    
end subroutine spline