How to resolve the algorithm Numerical integration/Gauss-Legendre Quadrature step by step in the Racket programming language
How to resolve the algorithm Numerical integration/Gauss-Legendre Quadrature step by step in the Racket programming language
Table of Contents
Problem Statement
For this, we first need to calculate the nodes and the weights, but after we have them, we can reuse them for numerious integral evaluations, which greatly speeds up the calculation compared to more simple numerical integration methods.
Task description Similar to the task Numerical Integration, the task here is to calculate the definite integral of a function
f ( x )
{\displaystyle f(x)}
, but by applying an n-point Gauss-Legendre quadrature rule, as described here, for example. The input values should be an function f to integrate, the bounds of the integration interval a and b, and the number of gaussian evaluation points n. An reference implementation in Common Lisp is provided for comparison. To demonstrate the calculation, compute the weights and nodes for an 5-point quadrature rule and then use them to compute:
Let's start with the solution:
Step by Step solution about How to resolve the algorithm Numerical integration/Gauss-Legendre Quadrature step by step in the Racket programming language
Source code in the racket programming language
(define (LegendreP n x)
(let compute ([n n] [Pn-1 x] [Pn-2 1])
(case n
[(0) Pn-2]
[(1) Pn-1]
[else (compute (- n 1)
(/ (- (* (- (* 2 n) 1) x Pn-1)
(* (- n 1) Pn-2)) n)
Pn-1)])))
(define (LegendreP′ n x)
(* (/ n (- (* x x) 1))
(- (* x (LegendreP n x))
(LegendreP (- n 1) x))))
(define (LegendreP-root n i)
; newton-raphson step
(define (newton-step x)
(- x (/ (LegendreP n x) (LegendreP′ n x))))
; initial guess
(define x0 (cos (* pi (/ (- i 1/4) (+ n 1/2)))))
; computation of a root with relative accuracy 1e-15
(if (< (abs x0) 1e-15)
0
(let next ([x′ (newton-step x0)] [x x0])
(if (< (abs (/ (- x′ x) (+ x′ x))) 1e-15)
x′
(next (newton-step x′) x′)))))
(define (Gauss-Legendre-quadrature n)
;; positive roots
(define roots
(for/list ([i (in-range (floor (/ n 2)))])
(LegendreP-root n (+ i 1))))
;; weights for positive roots
(define weights
(for/list ([x (in-list roots)])
(/ 2 (- 1 (sqr x)) (sqr (LegendreP′ n x)))))
;; all roots and weights
(values (append (map - roots)
(if (odd? n) (list 0) '())
(reverse roots))
(append weights
(if (odd? n) (list (/ 2 (sqr (LegendreP′ n 0)))) '())
(reverse weights))))
(define (integrate f a b #:nodes (n 5))
(define m (/ (+ a b) 2))
(define d (/ (- b a) 2))
(define-values [x w] (Gauss-Legendre-quadrature n))
(define (g x) (f (+ m (* d x))))
(* d (+ (apply + (map * w (map g x))))))
> (Gauss-Legendre-quadrature 5)
'(-0.906179845938664 -0.5384693101056831 0 0.5384693101056831 0.906179845938664)
'(0.23692688505618875 0.47862867049936625 128/225 0.47862867049936625 0.23692688505618875)
> (integrate exp -3 3)
20.035577718385547
> (- (exp 3) (exp -3)
20.035749854819805
> (require plot)
> (parameterize ([plot-x-label "Number of Gaussian nodes"]
[plot-y-label "Integration error"]
[plot-y-transform log-transform]
[plot-y-ticks (log-ticks #:base 10)])
(plot (points (for/list ([n (in-range 2 11)])
(list n (abs (- (integrate exp -3 3 #:nodes n)
(- (exp 3) (exp -3)))))))))
You may also check:How to resolve the algorithm Run-length encoding step by step in the Elena programming language
You may also check:How to resolve the algorithm Factors of a Mersenne number step by step in the J programming language
You may also check:How to resolve the algorithm Factors of an integer step by step in the ALGOL-M programming language
You may also check:How to resolve the algorithm Count in octal step by step in the Befunge programming language
You may also check:How to resolve the algorithm Brace expansion step by step in the Tcl programming language