poly_roots Subroutine

  1. pf_quadrature.f90
  2. pf_mod_quadrature
  3. poly_roots

public subroutine poly_roots(roots, p0, n)

Subroutine to compute polynomial roots using the Durand-Kerner algorithm. The roots are assumed to be real.

Arguments

Type IntentOptional AttributesName
real(kind=pfqp), intent(out) :: roots(n)
real(kind=pfqp), intent(in) :: p0(0:n)
integer, intent(in), value:: n

Calls

proc~~poly_roots~~CallsGraph proc~poly_roots poly_roots proc~qsort qsort proc~poly_roots->proc~qsort proc~poly_eval_complex poly_eval_complex proc~poly_roots->proc~poly_eval_complex

Called by

proc~~poly_roots~~CalledByGraph proc~poly_roots poly_roots proc~sdc_qnodes sdc_qnodes proc~sdc_qnodes->proc~poly_roots proc~pf_quadrature pf_quadrature proc~pf_quadrature->proc~sdc_qnodes proc~pf_init_sdcmats pf_init_sdcmats proc~pf_init_sdcmats->proc~pf_quadrature proc~pf_level_setup pf_level_setup proc~pf_level_setup->proc~pf_init_sdcmats proc~pf_pfasst_setup pf_pfasst_setup proc~pf_pfasst_setup->proc~pf_level_setup

Contents

Source Code

poly_roots

Source Code

  subroutine poly_roots(roots, p0, n) 
    integer,  intent(in), value   :: n
    real(pfqp),        intent(out)  :: roots(n)
    real(pfqp),        intent(in)   :: p0(0:n)

    integer     :: i, j, k
    complex(pfqp) :: num, den, z0(n), z1(n)
    real(pfqp)    :: p(0:n)
    
    real(pfqp) ::  eps 

    eps = epsilon(1.0_pfqp)*100.0_pfqp
    p = p0 / p0(n)

    ! initial guess
    do i = 1, n
       z0(i) = (0.4_pfqp, 0.9_pfqp)**i
    end do

    ! durand-kerner-weierstrass iterations
    z1 = z0
    do k = 1, 100
       do i = 1, n

          ! evaluate poly at z0(i)
          num = poly_eval_complex(p, n, z0(i))

          ! evaluate denominator
          den = 1.0_pfqp
          do j = 1, n
             if (j == i) cycle
             den = den * (z0(i) - z0(j))
          end do

          ! update
          z0(i) = z0(i) - num / den
       end do

       ! converged?
       if (sum(abs(z0 - z1)) < eps) exit

       z1 = z0
    end do

    roots = real(z0)
    where (abs(roots) < eps) roots = 0.0_pfqp
    call qsort(roots)

  end subroutine poly_roots