nach dem Verfahren von Bairstow bestimmt werden.
Gefundene Nullstellen werden abdividiert und das Verfahren fortgesetzt, bis eine lineare oder quadratische Gleichung übrig bleibt.
Gefundene Nullstellen sind durch Rundungs- und Abbruchfehler u.U. stark verfälscht. Deshalb muss jede Nullstelle nachiteriert werden, z.B. mit dem gewöhnlichem Newtonverfahren.
module mysubs contains subroutine init(a, b, c, n, eps, maxit) use kinds implicit none real(kind=REAL8), dimension(:), pointer :: a, b, c integer, intent(out) :: n, maxit real(kind=REAL8), intent(out) :: eps integer :: i write (*,*) 'Darstellung der Berechnung der Nullstellen eines reellen' write (*,*) 'Polynomes mittels des Verfahren von Bairstow' write (*,*) write (*,*) 'Geben Sie bitte den Grad des Polynoms an' read (*,*) n allocate (a(0:n)) allocate (b(0:n)) allocate (c(0:n)) do i=0, n write (*,*) 'Geben Sie nun den Koeffizienten a(', i, ') an' read (*,*) a(i) end do write (*,*) 'Wie genau sollen die Nullstellen berechnet werden?' read (*,*) eps do if (eps > 0) exit write (*,*) 'Nur positive Genauigkeiten machen Sinn. Versuchen Sie & &es noch einmal.' read (*,*) eps end do write (*,*) 'Wieviele Iterationen sollen maximal durchgefuehrt werden?' read (*,*) maxit do if (maxit > 0) exit write (*,*) 'Nur positive Werte sind hier vernuenftig. Versuchen & &Sie es noch einmal.' read (*,*) maxit end do return end subroutine init subroutine loese(r, s) use kinds implicit none real(kind=REAL8), intent(in) :: r, s real(kind=REAL8) :: x0, x1 if (r*r + 4*s < 0) then write (*,*) 'Das Polynom hat komplexe Wurzeln.' stop else x0 = r/2.0 - sqrt(r*r/4.0 + s) x1 = r/2.0 + sqrt(r*r/4.0 + s) write (*, '(F9.5)') x0 write (*, '(F9.5)') x1 end if return end subroutine loese end module mysubs program Bairstow ! Berechnung der Nullstellen eines reellen Polynoms use kinds use mysubs implicit none real(kind=REAL8), dimension(:), pointer :: a, b, c integer :: n, i, deg, iter, maxit real(kind=REAL8) :: r, s ! x^2-rx-s ist ein Teiler des Polynoms real(kind=REAL8) :: eps, dr, ds, eps1, eps2, det logical :: abbruch call init(a, b, c, n, eps, maxit) deg = n iter = 0 do if ((deg <= 2) .or. (iter >= maxit)) exit iter = 0 write (*,*) 'Geben Sie geschaetzte Koeffizienten fuer einen & &quadratischen Teiler des Polynoms an.' read (*,*) r, s do iter = iter + 1 b(deg) = a(deg) b(deg-1) = a(deg-1) + r*b(deg) c(deg) = b(deg) c(deg-1) = b(deg-1) + r*c(deg) i = deg - 2 do if (i < 0) exit b(i) = a(i) + r*b(i+1) + s*b(i+2) c(i) = b(i) + r*c(i+1) + s*c(i+2) i = i-1 end do det = c(2)*c(2) - c(3)*c(1) abbruch = (det == 0) .or. (iter >maxit) if (.not. abbruch) then dr = (-b(1)*c(2) + b(0)*c(3))/det ds = (-b(0)*c(2) + b(1)*c(1))/det eps1 = dr eps2 = ds r = r + dr s = s + ds end if if ( (eps1eps) .or. abbruch) exit end do if (.not. abbruch) then call loese(r, s) deg = deg - 2 do i=0, deg a(i) = b(i+2) end do end if end do if (iter < maxit) then if (deg == 2) then r = -a(1)/a(2) s = -a(0)/a(2) call loese(r, s) else write (*,*) -a(0)/a(1) end if else write (*,*) 'Es werden zu viele Schritte benoetigt.' end if stop end program Bairstow