Legendre polynomial generation in Fortran

by Yi Zhang

In spectral-based numerical method, orthogonal polynomial basis is commonly used to generate “nice” linear systems. A popular candidate of such polynomial family is Legendre polynomial.  As other Jacobi polynomials,  recursive definition makes it easy to generate a high order member.  The following is  a Fortran snippet to generate Legendre polynomial, following Bonnet’s recursion formula.

! Copyright (C) 2012 Yi Zhang

! Author: Yi Zhang <zhangy2@onid.orst.edu>
! Keywords: numerical, tools, polynomial

! This file is free software; you can redistribute it and/or modify
! it under the terms of the FreeBSD License. So essentially do
! whatever you want.

module legendre_poly

  type polycoef
     real, pointer :: p(:)
  end type polycoef

contains

  subroutine gen_legendre_poly()
 ! q_n = p_(n-1) = 1/(n-1)*[(2n-3)xp_(n-2)-(n-2)p_(n-3)]
 !     = 1/(n-1)*[(2n-3)xq_(n-1)-(n-2)q_(n-2)]
    implicit none
    integer, parameter :: n = 11
    type(polycoef) :: a(n)
    integer :: i, j

    DO i = 1, n
       allocate(a(i) % p(i))
       a(i) % p = 0
    ENDDO
    a(1) % p(1) = 1
    a(2) % p    = (/0,1/)
    do i = 3, n
       do j = 1, i-2
          a(i) % p(j) = -(i-2) * a(i-2) % p(j)
       enddo
       do j = 1, i-1
          a(i) % p(j+1) = a(i) % p(j+1) + (2*i-3) * a(i-1) % p(j)
       enddo
       a(i) % p = a(i) % p / (i-1)
       write(*,*) a(i) % p
    enddo
  end subroutine gen_legendre_poly

end module legendre_poly

DG methods usually can use up to order 10 for high order scheme. So this snippet can be used to generate polynomial for such applications’ basis function, usually projected to a canonical element.

Advertisements