|
|
|
|
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
|
|
|