Berechnung von Lagrange-Polynomen

Kein Verweis

Berechnung von Lagrange-Polynomen

Das Programm berechnet den Funktionswert an einer Stelle x0, wenn die Funktion an (n + 1) Stützstellen bekannt ist.
Zur Interpolation wird ein Lagrange-Polynom n-ten Grades durch die (n + 1) Stützstellen verwendet.

module mysubs

contains

  subroutine init(n, x, fx, x0)
    use kinds
    use fileutil
    implicit none
    real(kind=REAL8), dimension(:), pointer :: x, fx
    integer, intent(out)          :: n
    real(kind=REAL8), intent(out) :: x0

    integer, parameter  :: fileid = 1
    character(len=256)  :: fileName
    integer             :: status, i

    write (*,*) 'Berechnet ein Lagrange-Polynom n-ten Grades'
    write (*,*)
    write (*,*) 'In welcher Datei liegen die Stuetzstellen?&
          & (In lagrange.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(x(0:n))
    allocate(fx(0:n))

    read(fileid, *, iostat=status) (x(i), fx(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 (*,*) 'An welcher Stelle soll das Polynom berechnet werden? '
    read (*,*) x0

    return
  end subroutine init

  subroutine result(x, y)
    use kinds
    implicit none
    real(kind=REAL8), intent(in) :: x, y

    write (*,'(A, F10.2, A, F10.2, A)') 'Das Lagrange-Polynom an der Stelle ',&
         &   x, ' hat den Wert ', y, '.'

    return
  end subroutine result
  
end module mysubs


program Lagrange
  ! Berechnet ein Lagrange-Polynom n-ten Grades anhand von n + 1 
  ! vorgegebenen Stuetzstellen.

  use kinds
  use mysubs
  implicit none

  real(kind=REAL8), dimension(:), pointer :: x, fx
  real(kind=REAL8)  :: x0       ! an dieser Stelle wird das Polynom berechnet
  real(kind=REAL8)  :: y        ! enthaelt das Ergebnis
  real(kind=REAL8)  :: produkt  ! Zwischenergebnis: der i-te Produktterm
  integer           :: n, i, j

  call init(n, x, fx, x0)
  y = 0                         ! Ergebnis initialisieren
  do i = 0, n 
     produkt = fx(i)
     do j = 0, n
        if (i /= j) then
           produkt = produkt * (x0 - x(j)) / (x(i) - x(j))
        end if
     end do
     y = y + produkt         
  end do
  call result(x0, y)
  stop
end program Lagrange