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