Skip to content

Instantly share code, notes, and snippets.

@Thraetaona
Last active November 26, 2024 17:33
Show Gist options
  • Save Thraetaona/7efb4f21adf2644fe833bc4184484324 to your computer and use it in GitHub Desktop.
Save Thraetaona/7efb4f21adf2644fe833bc4184484324 to your computer and use it in GitHub Desktop.
FORTRAN 90 code to compute the Spline for a given set of x and f(x) arrays.
! ----------------------------------------------------------------------
! SPDX-License-Identifier: CC0-1.0
!
! spline.f90: Spline Interpolation
! Authored by Fereydoun Memarzanjany
!
! To the maximum extent possible under law, the Author waives all
! copyright and related or neighboring rights to this code.
!
! You should have received a copy of the CC0 legalcode along with this
! work; if not, see <http://creativecommons.org/publicdomain/zero/1.0/>
! ----------------------------------------------------------------------
program spline_coefficients
implicit none ! I want types to be explicit, not Hungarian notation.
integer, parameter :: n = 6 ! Number of data points
real :: xj(n), fj(n) ! Points and function values
real :: hj(n-1), aj(n), bj(n-1), cj(n), dj(n-1)
real :: alpha(n-1), l(n), mu(n), z(n)
integer :: i ! Counter for loop indices (1-indexed)
! Input data points
xj = (/ 0.0, 1.0, 2.0, 10.0, 50.0, 100.0 /)
fj = (/ 100.0, 45.0, 39.0, 22.0, 5.0, 0.05 /)
! Compute hj (i.e., heights or distances between intervals)
do i = 1, n-1
hj(i) = xj(i+1) - xj(i)
end do
! ------------------------------------------------------------------
! Solve the Sparse Matrix
! ------------------------------------------------------------------
! Compute alpha (temporary variable)
alpha(1) = 0.0
do i = 2, n-1
alpha(i) = (3.0/hj(i) ) * (fj(i+1) - fj(i) ) &
- (3.0/hj(i-1)) * (fj(i) - fj(i-1))
end do
! Tridiagonal system initialization
l(1) = 1.0 ! Diagonal Element
mu(1) = 0.0 ! Multipliers (Greek \Mu letter)
z(1) = 0.0 ! Temporary holders
! Forward elimination
do i = 2, n-1
l(i) = 2.0 * (xj(i+1) - xj(i-1)) - hj(i-1) * mu(i-1)
mu(i) = hj(i) / l(i)
z(i) = (alpha(i) - hj(i-1) * z(i-1)) / l(i)
end do
! Boundary conditions
l(n) = 1.0
z(n) = 0.0
cj(n) = 0.0
! Back substitution
do i = n-1, 1, -1
cj(i) = z(i) - mu(i) * cj(i+1)
bj(i) = (fj(i+1) - fj(i)) / hj(i) &
- (hj(i) * (cj(i+1) + 2.0 * cj(i))) / 3.0
dj(i) = (cj(i+1) - cj(i)) / (3.0 * hj(i))
end do
cj(1) = 0.0
! ------------------------------------------------------------------
! Output the coefficients
print '(A)', 'j cj bj dj'
do i = 1, n-1
print '(I2,2X,F10.6,2X,F10.6,2X,F12.6)', &
i-1, cj(i), bj(i), dj(i)
end do
! Print the last cj alone.
print '(I2,2X,F10.6)', &
n-1, cj(n)
end program spline_coefficients
! ----------------------------------------------------------------------
! END OF FILE: spline.f90
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
! Sample Output:
! ----------------------------------------------------------------------
! j cj bj dj
! 0 0.000000 -67.375244 12.375247
! 1 37.125740 -30.249504 -12.876235
! 2 -1.502965 5.373271 0.070710
! 3 0.194078 -5.097822 -0.001931
! 4 -0.037695 1.157507 0.000251
! 5 0.000000
! ----------------------------------------------------------------------
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment