# Gauss-Legendre Algorithm

## Introduction

The Gauss-Legendre algorithm is an iterative numerical integration technique used for approximating definite integrals. The algorithm is based on the idea of interpolating a function over a given interval using Legendre polynomials. The resulting approximation is then used to compute the integral of the function over the same interval. The Gauss-Legendre algorithm is widely used in numerical analysis and scientific computing.

## Implementation

```
let gauss_legendre ?(eps = 3e-11) ?(a = -1.) ?(b = 1.) n =
let m = (n + 1) / 2 in
let n' = float_of_int n in
let x = Array.create_float n in
let w = Array.create_float n in
let xm = 0.5 *. (b +. a) in
let xl = 0.5 *. (b -. a) in
let p1 = ref infinity in
let p2 = ref infinity in
let p3 = ref infinity in
let pp = ref infinity in
let z = ref infinity in
for i = 1 to m do
let i' = float_of_int i in
z := cos (Owl_const.pi *. (i' -. 0.25) /. (n' +. 0.5));
(try
while true do
p1 := 1.;
p2 := 0.;
for j = 1 to n do
p3 := !p2;
p2 := !p1;
let j' = float_of_int j in
p1 := ((((2. *. j') -. 1.) *. !z *. !p2) -. ((j' -. 1.) *. !p3)) /. j'
done;
pp := n' *. ((!z *. !p1) -. !p2) /. ((!z *. !z) -. 1.);
let z1 = !z in
z := z1 -. (!p1 /. !pp);
assert (abs_float (!z -. z1) > eps)
done
with
| _ -> ());
x.(i - 1) <- xm -. (xl *. !z);
x.(n - i) <- xm +. (xl *. !z);
w.(i - 1) <- 2. *. xl /. ((1. -. (!z *. !z)) *. !pp *. !pp);
w.(n - i) <- w.(i - 1)
done;
x, w
let gauss_legendre_cache = Array.init 50 gauss_legendre
let _gauss_laguerre ?(_eps = 3e-11) _a _b _n = ()
let gaussian_fixed ?(n = 10) f a b =
let x, w =
match n < Array.length gauss_legendre_cache with
| true -> gauss_legendre_cache.(n)
| false -> gauss_legendre n
in
let xr = 0.5 *. (b -. a) in
let s = ref 0. in
for i = 0 to n - 1 do
let c = (xr *. (x.(i) +. 1.)) +. a in
s := !s +. (w.(i) *. f c)
done;
!s *. xr
let gaussian ?(n = 50) ?(eps = 1e-6) f a b =
let s_new = ref infinity in
let s_old = ref infinity in
(try
for i = 1 to n do
s_new := gaussian_fixed ~n:i f a b;
assert (abs_float (!s_new -. !s_old) > eps);
s_old := !s_new
done
with
| _ -> ());
!s_new
```

The Gauss-Legendre algorithm is implemented in OCaml as a function `gaussian`

which takes in the following arguments:

`n`

: the number of iterations to perform (default 50)`eps`

: the desired accuracy of the approximation (default 1e-6)`f`

: the function to integrate`a`

: the lower bound of the integration interval`b`

: the upper bound of the integration interval

The implementation is split into two functions: `gaussian_fixed`

and `gaussian`

. The former performs a single iteration of the algorithm, while the latter performs a specified number of iterations until the desired accuracy is achieved. The `gaussian`

function first checks if the desired number of iterations is available in a precomputed cache, otherwise it computes them using the `gauss_legendre`

function.

## Step-by-step Explanation

The Gauss-Legendre algorithm works by approximating the integral of a function `f`

over an interval `[a,b]`

using Legendre polynomials. The algorithm proceeds as follows:

- Compute the Legendre-Gauss nodes and weights for the given interval
`[a,b]`

and number of nodes`n`

. - Transform the integral over
`[a,b]`

into an integral over`[-1,1]`

using the change of variables`x = ((b-a)t + b + a)/2`

. - Approximate the integral over
`[-1,1]`

using the Legendre-Gauss nodes and weights. - Transform the approximation back to the integral over
`[a,b]`

using the same change of variables as in step 2.

The Legendre-Gauss nodes and weights are computed using the `gauss_legendre`

function. This function uses the iterative method described by the algorithm to approximate the roots of the Legendre polynomial of degree `n`

. The roots are then used as the nodes for the Gauss-Legendre quadrature, and the corresponding weights are computed using the formula `w_i = 2 / [(1 - z_i^2) * (P_n'(z_i))^2]`

, where `P_n`

is the Legendre polynomial of degree `n`

and `z_i`

is the `i`

-th root.

The `gaussian_fixed`

function approximates the integral over `[-1,1]`

using the Legendre-Gauss nodes and weights. It first computes the midpoint `xr`

of the interval `[-1,1]`

and then evaluates the function `f`

at each of the nodes `x_i`

transformed to the interval `[a,b]`

using the change of variables described above. The approximation is then computed as the sum of the products of the weights `w_i`

and the function values `f(x_i)`

.

The `gaussian`

function iteratively calls `gaussian_fixed`

with increasing numbers of nodes until the desired accuracy is achieved. It uses the precomputed cache of Legendre-Gauss nodes and weights if available, otherwise it computes them using the `gauss_legendre`

function.

## Complexity Analysis

The time complexity of the Gauss-Legendre algorithm is dominated by the number of iterations `n`

performed by the `gaussian`

function. Each iteration of `gaussian`

requires evaluating the function `f`

at `n`

points and computing the Legendre-Gauss nodes and weights. The latter is precomputed and cached for up to 50 nodes, and computed on-the-fly using the `gauss_legendre`

function for larger values of `n`

.

The time complexity of computing the Legendre-Gauss nodes and weights using `gauss_legendre`

is `O(n^2)`

, due to the nested loop over `j`

and `i`

. However, this function is only called once for each unique `n`

, so its time complexity is not a bottleneck for the overall algorithm.

The time complexity of evaluating the function `f`

at `n`

points is dependent on the implementation of `f`

. Assuming that `f`

has a time complexity of `O(1)`

for each evaluation, the time complexity of evaluating `f`

at `n`

points is `O(n)`

.

Therefore, the overall time complexity of the Gauss-Legendre algorithm is `O(n^3)`

for large values of `n`

. The space complexity is `O(n)`

for storing the Legendre-Gauss nodes and weights.