gjp_rec Module Reference

This is derived from the recurrence relations in chebfun. More...

Functions/Subroutines
subroutine	gauss_jacobi_rec (npts, alpha, beta, x, wts)
subroutine	recurrence (npts, n2, alpha, beta, x, PP)
subroutine	eval_jacobi_poly (x, npts, alpha, beta, P, Pp)

Detailed Description

This is derived from the recurrence relations in chebfun.

Function/Subroutine Documentation

eval_jacobi_poly()

subroutine gjp_rec::eval_jacobi_poly	(	real(dp), dimension(:), intent(in)	x,
		integer, intent(in)	npts,
		real(dp), intent(in)	alpha,
		real(dp), intent(in)	beta,
		real(dp), dimension(size(x)), intent(out)	P,
		real(dp), dimension(size(x)), intent(out)	Pp
	)

Definition at line 91 of file gjp_rec.f90.

    integer, intent(in) :: npts
    real(dp), intent(in) :: alpha, beta
    real(dp), intent(in) :: x(:)
    real(dp), dimension(size(x)), intent(out) :: P, Pp
    real(dp), dimension(size(x)) :: Pm1, Ppm1, Pa1, Ppa1
    integer :: k, i
    real(dp) :: A_val, B_val, C_val, D_val
 
    p = 0.5_dp * (alpha - beta + (alpha + beta + 2) * x)
    pm1 = 1.0_dp
    pp = 0.5_dp * (alpha + beta + 2)
    ppm1 = 0.0_dp
 
    if (npts == 0) then
        p = pm1
        pp = ppm1
    end if
 
    do k = 1, npts - 1
        a_val = 2 * (k + 1) * (k + alpha + beta + 1) * (2 * k + alpha + beta)
        b_val = (2 * k + alpha + beta + 1) * (alpha**2 - beta**2)
        c_val = product([(2 * k + alpha + beta + i, i=0, 2)])
        d_val = 2 * (k + alpha) * (k + beta) * (2 * k + alpha + beta + 2)
 
        pa1 = ((b_val + c_val * x) * p - d_val * pm1) / a_val
        ppa1 = ((b_val + c_val * x) * pp + c_val * p - d_val * ppm1) / a_val
 
        pm1 = p
        p = pa1
        ppm1 = pp
        pp = ppa1
    end do

gauss_jacobi_rec()

subroutine gjp_rec::gauss_jacobi_rec	(	integer, intent(in)	npts,
		real(dp), intent(in)	alpha,
		real(dp), intent(in)	beta,
		real(dp), dimension(npts), intent(out)	x,
		real(dp), dimension(npts), intent(out)	wts
	)

Definition at line 35 of file gjp_rec.f90.

    integer, intent(in) :: npts
    real(dp), intent(in) :: alpha, beta
    real(dp), intent(out) :: x(npts), wts(npts)
    real(dp), dimension(ceiling(npts/2._dp)) :: x1, ders1
    real(dp), dimension(npts/2) :: x2, ders2
    real(dp) :: ders(npts), C
    integer :: idx
    call recurrence(npts, ceiling(npts / 2._dp), alpha, beta, x1, ders1)
    call recurrence(npts, npts / 2, beta, alpha, x2, ders2)
    do idx = 1, npts / 2
        x(idx) = -x2(npts / 2 - idx + 1)
        ders(idx) = ders2(npts / 2 - idx + 1)
    end do
    do idx = 1, ceiling(npts / 2._dp)
        x(npts / 2 + idx) = x1(idx)
        ders(npts / 2 + idx) = ders1(idx)
    end do
    wts = 1.0_dp / ((1.0_dp - x**2) * ders**2)
    c = 2**(alpha + beta + 1) * exp(log_gamma(npts + alpha + 1) - &
                                    log_gamma(npts + alpha + beta + 1) + &
                                    log_gamma(npts + beta + 1) - log_gamma(npts + 1._dp))
    wts = wts * c

recurrence()

subroutine gjp_rec::recurrence	(	integer, intent(in)	npts,
		integer, intent(in)	n2,
		real(dp), intent(in)	alpha,
		real(dp), intent(in)	beta,
		real(dp), dimension(n2), intent(out)	x,
		real(dp), dimension(n2), intent(out)	PP
	)

Definition at line 60 of file gjp_rec.f90.

    integer, intent(in) :: npts, n2
    real(dp), intent(in) :: alpha, beta
    real(dp), intent(out) :: x(n2), PP(n2)
    real(dp) :: dx(n2), P(n2)
    integer :: r(n2), l, i
    real(dp) :: C, T
 
    do i = 1, n2
        r(i) = n2 - i + 1
    end do
 
    do i = 1, n2
        c = (2 * r(i) + alpha - 0.5_dp) * pi / (2 * npts + alpha + beta + 1)
     t = c + 1 / (2 * npts + alpha + beta + 1)**2 * ((0.25_dp - alpha**2) / tan(0.5_dp * c) - (0.25_dp - beta**2) * tan(0.5_dp * c))
        x(i) = cos(t)
    end do
 
    dx = 1.0_dp
    l = 0
 
    do while (maxval(abs(dx)) > sqrt(epsilon(1.0_dp)) / 1000 .and. l < 10)
        l = l + 1
        call eval_jacobi_poly(x, npts, alpha, beta, p, pp)
        dx = -p / pp
        x = x + dx
    end do
 
    call eval_jacobi_poly(x, npts, alpha, beta, p, pp)