Last active
November 26, 2024 17:33
-
-
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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
! ---------------------------------------------------------------------- | |
! 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