Created
December 17, 2024 23:53
-
-
Save Thraetaona/804578aa54534804475b41a82fcfadf3 to your computer and use it in GitHub Desktop.
Spline interpolation matrices in FORTRAN
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 | |
! 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 | |
! ---------------------------------------------------------------------- |
Author
Thraetaona
commented
Dec 17, 2024
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment