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