Vorige Seite Kein Verweis Übersicht Kein Verweis content index Hilfeseiten  

Berechnung eines interpolierenden kubischen Splines

Berechnung eines interpolierenden kubischen Splines mit natürlichen Randbedingungen.
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