Gauss-Legendre Quadrature is a numerical integration technique designed to yield maximum algebraic precision with a minimal number of function evaluations. Unlike standard Newton-Cotes methods (such as the Trapezoidal or Simpson's rules) that use fixed, equally spaced evaluation points, Gauss-Legendre quadrature treats both the evaluation locations (nodes) and their corresponding multipliers (weights) as free variables.
Mathematical Foundation To evaluate a definite integral over a standardized interval [−1,1]:
∫−11f(x),dx≈i=1∑nwif(xi)
The nodes xi are selected as the roots of the n-th degree Legendre polynomial Pn(x). The weights wi are calculated by integrating the Lagrange interpolating polynomials over the interval:
wi=(1−xi2)[Pn′(xi)]22
An n-point Gauss-Legendre rule integrates polynomials of degree up to 2n - 1 exactly.
Domain Transformation (Mapping) For integrals over an arbitrary real interval [a, b], a linear change of variables maps the domain to [-1, 1]:
x(t)=2b−at+2a+b
dx=2b−adt
The transformation yields:
∫abf(x),dx=2b−ai=1∑nwif(2b−axi+2a+b)
Applied Examples (Solved via Simple Solver)
.. Admonition:: Example 1 : Exact Polynomial Integration (2-Point Rule)
Evaluate
∫02(x3−3x+2),dx
Example 1C#
1
2
var I = Integral(x => Polyval([1.0, 0.0, -3.0, 2.0], x), 0, 2);
Console.WriteLine($"I = {I:F4}");
OutputFrom the book
I = 2.0000
Exact Analytical Result:
[4x4−23x2+2x]02=(4−6+4)−0=2
.. Admonition:: Example 2 : Transcendental Function
∫0πsin(x),dx
Example 2C#
1
2
var I = Integral(x => Sin(x), 0, pi);
Console.WriteLine($"I = {I:F4}");
OutputFrom the book
I = 2.0000
Exact Analytical Result:
[−cos(x)]0π=1−(−1)=2
.. Admonition:: Example 3 : Infinite bound
Evaluate
∫0∞e−x2(ln(x))2,dx
Example 3C#
1
2
var I = Integral(x => Exp(-x * x) * Pow(Log(x), 2), 0, inf);
Console.WriteLine($"I = {I:F5}");
OutputFrom the book
I = 1.94752
Exact Analytical Result:
∫0∞e−x2(ln(x))2,dx=16π(π2+2(γ+ln(4))2)=1.94752
Multiple Integral
SepalSolver can handle multiple integral via IntegralN. IntegralN is built on Integral running over multiply direction recursively. It is abole to compute upto 4 dimensional integrals efficiently.
.. Admonition:: Example 4 : Hypercube domain in 3 dimensions
Integrate the function f(x,y,z)=xyz over the region where x ranges from 0 to 1, y ranges from 1 to 2, and z ranges from 2 to 3, which can be expressed as:
∫01∫12∫23xyzdzdydx
Example 4C#
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// Define the function to integrate
Func<double, double, double, double> f = (x, y, z) => x * y * z;
// Set the lower bound of x
double x_1 = 0;
// Set the upper bound of x
double x_2 = 1;
// Set the lower bound of y
double y_1 = 1;
// Set the upper bound of y
double y_2 = 2;
// Set the lower bound of z
double z1 = 2;
// Set the upper bound of z
double z2 = 3;
// Calculate the integral
double integral = Integral3(f, x_1, x_2, y_1, y_2, z1, z2);
// Print the result
Console.WriteLine($"The triple integral of x*y*z is approximately: {integral:F4}");
OutputFrom the book
The triple integral of x*y*z is approximately: 1.8750
.. Admonition:: Example 5 :
Integrate the function f(x,y)=xy over the region where x ranges from 0 to 1, y ranges from x2 to 2, which can be expressed as:
∫01∫x22xydydx
Example 5C#
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
// Define the function to integrate
Func<double, double, double> f = (x, y) => x * y;
// Set the lower bound of x
double x_1 = 0;
// Set the upper bound of x
double x_2 = 1;
// Define the lower bound of y as a function of x
Func<double, double> y_1 = x => x * x;
// Set the upper bound of y
double y_2 = 2;
// Calculate the integral
double integral = Integral2(f, x_1, x_2, y_1, y_2);
// Print the result
Console.WriteLine($"The triple integral of x*y*z is approximately: {integral:F4}");
OutputFrom the book
The triple integral of x*y*z is approximately: 0.9167
.. Admonition:: Example 6 : Functional boundary
Integrate the function f(x,y,z)=xyz over the region where x ranges from 0 to 1, y ranges from x2 to x, and z ranges from 2 to 3, which can be expressed as:
∫01∫x2x∫23xyzdzdydx
Example 6C#
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// Define the function to integrate
Func<double, double, double, double> f = (x, y, z) => x * y * z;
// Set the lower bound of x
double x_1 = 0;
// Set the upper bound of x
double x_2 = 1;
// Set the lower bound of y as a function of x
Func<double, double> y_1 = x => x * x;
// Define the upper bound of y as a function of x
Func<double, double> y_2 = x => Sqrt(x);
// Set the lower bound of z
double z_1 = 2;
// Set the upper bound of z
double z_2 = 3;
// Calculate the integral
double integral = Integral3(f, x_1, x_2, y_1, y_2, z_1, z_2);
// Print the result
Console.WriteLine($"The triple integral of x*y*z is approximately: {integral:F4}");
OutputFrom the book
The triple integral of x*y*z is approximately: 0.2083
.. Admonition:: Example 7 : Functional boundary
Integrate the function f(x,y,z)=xyz over the region where x ranges from 0 to 1, y ranges from x2 to x, and z ranges from xy to x+y, which can be expressed as:
∫x1x2∫x2x∫xyx+yxyzdzdydx
Example 7C#
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// Define the function to integrate
Func<double, double, double, double> f = (x, y, z) => x * y * z;
// Set the lower bound of x
double x_1 = 0;
// Set the upper bound of x
double x_2 = 1;
// Define the lower bound of y as a function of x
Func<double, double> y_1 = x => x * x;
// Define the upper bound of y as a function of x
Func<double, double> y_2 = x => Sqrt(x);
// Define the lower bound of z as a function of x and y
double integral = Integral3(f, x_1, x_2, y_1, y_2, z_1, z_2);
// Print the result
Console.WriteLine($"The triple integral of x*y*z is approximately: {integral:F4}");
OutputFrom the book
The triple integral of x*y*z is approximately: 0.0641
The hard way, then the Reineur Solver way
How many points do you need?
Integrate e^(−x²) from 0 to 3. It has no closed form, so every answer is an approximation — the only question is how good, and at what cost.
By handThe trapezoidal rule with a fixed width. Add more points, get more digits.
1Func<double, double> f = x => Exp(-x * x);2double a = 0, b = 3;34foreach (int n in new[] { 10, 100, 1000, 10000 })5{6double h = (b - a) / n, sum = 0;7for (int i = 1; i < n; i++) sum += f(a + i * h);8double area = (h / 2) * (f(a) + 2 * sum + f(b));9 Console.WriteLine($"n = {n,-6} area = {area:F12}");10}
Output
n = 10 area = 0.886202033637
n = 100 area = 0.886207292750
n = 1000 area = 0.886207347704
n = 10000 area = 0.886207348254
With Reineur SolverGauss-Legendre picks its own points, and picks them well.
1double area = Integral(x => Exp(-x * x), 0, 3);2Console.WriteLine($"Integral area = {area:F12}");