Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

Computation of the Legendre polynomials and derivatives:

 

(define (LegendreP n x)

(let compute ([n n] [Pn-1 x] [Pn-2 1])

Computation of the Legendre polynomial roots:

 

(define (LegendreP-root n i)

; newton-raphson step

(define (newton-step x)

; initial guess

(define x0 (cos (* pi (/ (- i 1/4) (+ n 1/2)))))

(if (< (abs (/ (- x′ x) (+ x′ x))) 1e-15)

x′

 

Computation of Gauss-Legendre nodes and weights

 

(define (Gauss-Legendre-quadrature n)

(define roots

(for/list ([i (in-range (floor (/ n 2)))])

(LegendreP-root n (+ i 1))))

(define weights

(for/list ([x (in-list roots)])

(/ 2 (- 1 (sqr x)) (sqr (LegendreP′ n x)))))

</lang>

 

Integration using Gauss-Legendre quadratures:

 

(define (integrate f a b #:nodes (n 5))

(define m (/ (+ a b) 2))

(define d (/ (- b a) 2))

(define (g x) (f (+ m (* d x))))

(* d (+ (apply + (map * w (map g x))))))

 

Usage:

 

> (Gauss-Legendre-quadrature 5)

'(-0.906179845938664 -0.5384693101056831 0 0.5384693101056831 0.906179845938664)

Accuracy of the method:

 

> (require plot)

> (parameterize ([plot-x-label "Number of Gaussian nodes"]

[plot-y-label "Integration error"]

[plot-y-ticks (log-ticks #:base 10)])

</lang>

[[File:gauss.png]]