module mysubs contains subroutine init(n, t, y, x, u, v, w, nd, hd, c) use kinds use fileutil implicit none real(kind=REAL8), dimension(:), pointer :: t, y real(kind=REAL8), dimension(:), pointer :: u, v, w, nd, hd, c integer, intent(out) :: n real(kind=REAL8), intent(out) :: x integer, parameter :: fileid = 1 character(len=256) :: fileName integer :: status, i write (*,*) 'Berechnet ein interpoliertes kubisches natuerliches Spline' write (*,*) write (*,*) 'Welche Datei enthaelt die Punktepaare?& & (In splines.dat liegt ein Beispiel.)' read (*,*) fileName open(unit=fileId, file=fileName, action='read', iostat=status) if (status /= 0) then write (*,*) 'Kann Datei ', fileName, 'nicht öffnen.' stop end if n = anzahlReals2(fileId) - 1 allocate(t(0:n), y(0:n)) allocate(u(0:n), v(0:n), w(0:n), nd(0:n), hd(0:n), c(0:n)) read(fileid, *, iostat=status) (t(i), y(i), i=0, n) if (status /= 0) then write (*,*) 'Konnte Array nicht lesen.' stop end if close(unit=fileid) write (*,*) 'Insgesamt ', n + 1, ' Stuetzstellen gelesen.' write (*,*) write (*, '(A)', advance='no') & & 'An welcher Stelle soll das Spline berechnet werden? ' read (*,*) x write (*,*) return end subroutine init subroutine Randbedingungen(c, u, v, w, y, t) use kinds implicit none real(kind=REAL8), dimension(:), pointer :: c, u, v, w, y, t real(kind=REAL8) :: h1, hn, d1, dn, b1, b2 real(kind=REAL8) :: alpha, beta, det real(kind=REAL8) :: a11, a12, a21, a22 integer :: i, n n = size(t) - 1 h1 = t(1) - t(0) hn = t(n) - t(n - 1) d1 = (y(1) - y(0)) / h1 dn = (y(n) - y(n - 1)) / hn b1 = 3 * d1 - u(1) b2 = 3 * dn - u(n - 1) a11 = 2 * v(1) a12 = w(1) a21 = v(n - 1) a22 = 2 + w(n - 1) det = a11 * a22 - a12 * a21 alpha = (b1 * a22 - a12 * b2) / det beta = (a11 * b2 - b1 * a21) / det c(0) = alpha do i = 1, n - 1 c(i) = u(i) + alpha * v(i) + beta * w(i) end do c(n) = beta return end subroutine Randbedingungen function Spline(x, t, y, c) use kinds implicit none real(kind=REAL8), dimension(:), pointer :: t, y, c real(kind=REAL8), intent(in) :: x real(kind=REAL8) :: Spline real(kind=REAL8) :: h, delta, a0, a1, a2, a3 integer :: k, n n = size(t) - 1 k = n + 1 do k = k - 1 if (x > t(k)) exit end do k = k + 1 h = t(k) - t(k - 1) delta = (y(k) - y(k - 1)) / h a0 = y(k) a1 = c(k) a2 = (c(k) - delta) / h a3 = (c(k) - 2 * delta + c(k - 1)) / (h*h) Spline = ((a3 * (x - t(k - 1)) + a2) * (x - t(k)) + a1) * (x- t(k)) + a0 return end function Spline subroutine result(x, t, y, c) use kinds implicit none real(kind=REAL8), dimension(:), pointer :: t, y, c real(kind=REAL8), intent(in) :: x write (*,'(A, F12.5)') 'Das Ergenis hat den Wert ', Spline(x, t, y, c) return end subroutine result end module mysubs program Splines ! Berechnung eines interpolierten kubischen natuerlichen Splines} use kinds use mysubs implicit none real(kind=REAL8), dimension(:), pointer :: t, y real(kind=REAL8), dimension(:), pointer :: u, v, w, nd, hd, c real(kind=REAL8) :: x, h0, h1, a, b, q integer :: n, i call init(n, t, y, x, u, v, w, nd, hd, c) h0 = t(1) - t(0) a = (y(1) - y(0)) / (h0*h0) v(1) = -1 / h0 do i = 1, n - 1 h1 = t(i + 1) - t(i) b = (y(i + 1) - y(i)) / (h1*h1) u(i) = 3 * (a + b) hd(i) = 2 * (1 / h0 + 1 / h1) nd(i) = 1 / h1 a = b h0 = h1 end do w(n - 1) = -1 / h1 ! Vorwaertselimination do i = 1, n - 2 q = -nd(i) / hd(i) hd(i + 1) = hd(i + 1) + nd(i) * q u(i + 1) = u(i + 1) + q * u(i) v(i + 1) = q * v(i) end do ! Rueckwaertselimination do i = n - 1, 2, -1 q = hd(i) u(i) = u(i) / q v(i) = v(i) / q w(i) = w(i) / q q = nd(i - 1) u(i - 1) = u(i - 1) - q * u(i) v(i - 1) = v(i - 1) - q * v(i) w(i - 1) = -q * w(i) end do q = hd(1) u(1) = u(1) / q v(1) = v(1) / q w(1) = w(1) / q call Randbedingungen(c, u, v, w, y, t) call result(x, t, y, c) stop end program Splines