Gaussian Quadrature is a numerical integration method that provides an optimal way to choose
abscissas (evaluation points) and weights for approximating definite integrals. Unlike simpler
methods like the trapezoidal rule or Simpson's rule, Gaussian Quadrature can achieve higher
accuracy with fewer function evaluations by carefully selecting both the points and their weights.
How It Works
Transformation:
Transform the integral from [a,b] to [-1,1]
Use change of variables: x = ((b-a)t + (b+a))/2
Adjust function and differential accordingly
Point Selection:
Use roots of Legendre polynomials as abscissas
These points are optimal for polynomial integration
Points are symmetric around zero
Weight Calculation:
Weights are pre-computed for efficiency
Based on Legendre polynomial properties
Weights sum to 2 on [-1,1] interval
Integration:
Evaluate function at selected points
Multiply by corresponding weights
Sum the products and scale result
History
Gaussian Quadrature was developed by Carl Friedrich Gauss in the early 19th century. The method
emerged from his work on numerical integration and orthogonal polynomials. Gauss discovered that
by choosing specific points and weights, he could achieve exact integration for polynomials of
much higher degree than was possible with equally spaced points. This breakthrough significantly
improved the efficiency of numerical integration and has since become a fundamental tool in
numerical analysis and computational physics.
import numpy as np
def gaussian_quadrature(f, a, b, n):
"""
Compute definite integral using Gaussian quadrature.
Args:
f: Function to integrate
a, b: Integration interval [a,b]
n: Number of points (2-5)
Returns:
Approximation of the definite integral
"""
# Pre-computed Gauss-Legendre points and weights
points = {
2: [-0.5773502692, 0.5773502692],
3: [-0.7745966692, 0.0000000000, 0.7745966692],
4: [-0.8611363116, -0.3399810436, 0.3399810436, 0.8611363116],
5: [-0.9061798459, -0.5384693101, 0.0000000000, 0.5384693101, 0.9061798459]
}
weights = {
2: [1.0000000000, 1.0000000000],
3: [0.5555555556, 0.8888888889, 0.5555555556],
4: [0.3478548451, 0.6521451549, 0.6521451549, 0.3478548451],
5: [0.2369268850, 0.4786286705, 0.5688888889, 0.4786286705, 0.2369268850]
}
if n not in points:
raise ValueError("Number of points must be between 2 and 5")
# Transform points from [-1,1] to [a,b]
x = [(b-a)*t/2 + (b+a)/2 for t in points[n]]
w = [(b-a)*w/2 for w in weights[n]]
# Compute integral
result = sum(w[i] * f(x[i]) for i in range(n))
return result
# Example usage
if __name__ == "__main__":
# Example function: x²
def f(x): return x * x
# Integrate from -1 to 1
result = gaussian_quadrature(f, -1, 1, 3)
print(f"∫x² dx from -1 to 1 ≈ {result:.6f}")
print(f"Exact value: {2/3:.6f}")
#include <vector>
#include <map>
#include <functional>
#include <stdexcept>
#include <iostream>
class GaussianQuadrature {
private:
std::map<int, std::vector<double>> points;
std::map<int, std::vector<double>> weights;
public:
GaussianQuadrature() {
// Initialize Gauss-Legendre points and weights
points[2] = {-0.5773502692, 0.5773502692};
points[3] = {-0.7745966692, 0.0000000000, 0.7745966692};
points[4] = {-0.8611363116, -0.3399810436,
0.3399810436, 0.8611363116};
points[5] = {-0.9061798459, -0.5384693101, 0.0000000000,
0.5384693101, 0.9061798459};
weights[2] = {1.0000000000, 1.0000000000};
weights[3] = {0.5555555556, 0.8888888889, 0.5555555556};
weights[4] = {0.3478548451, 0.6521451549,
0.6521451549, 0.3478548451};
weights[5] = {0.2369268850, 0.4786286705, 0.5688888889,
0.4786286705, 0.2369268850};
}
double integrate(std::function<double(double)> f,
double a, double b, int n) {
if (points.find(n) == points.end()) {
throw std::invalid_argument(
"Number of points must be between 2 and 5"
);
}
double result = 0.0;
// Transform points and weights
for (int i = 0; i < n; i++) {
double x = (b-a)/2 * points[n][i] + (b+a)/2;
double w = (b-a)/2 * weights[n][i];
result += w * f(x);
}
return result;
}
};
int main() {
// Example function: x²
auto f = [](double x) { return x * x; };
GaussianQuadrature quad;
try {
double result = quad.integrate(f, -1.0, 1.0, 3);
std::cout << "∫x² dx from -1 to 1 ≈ "
<< result << std::endl;
std::cout << "Exact value: "
<< 2.0/3.0 << std::endl;
}
catch (const std::exception& e) {
std::cerr << "Error: " << e.what() << std::endl;
}
return 0;
}
public class GaussianQuadrature {
private readonly Dictionary<int, double[]> points;
private readonly Dictionary<int, double[]> weights;
public GaussianQuadrature() {
points = new Dictionary<int, double[]> {
{2, new[] {-0.5773502692, 0.5773502692}},
{3, new[] {-0.7745966692, 0.0000000000, 0.7745966692}},
{4, new[] {-0.8611363116, -0.3399810436,
0.3399810436, 0.8611363116}},
{5, new[] {-0.9061798459, -0.5384693101, 0.0000000000,
0.5384693101, 0.9061798459}}
};
weights = new Dictionary<int, double[]> {
{2, new[] {1.0000000000, 1.0000000000}},
{3, new[] {0.5555555556, 0.8888888889, 0.5555555556}},
{4, new[] {0.3478548451, 0.6521451549,
0.6521451549, 0.3478548451}},
{5, new[] {0.2369268850, 0.4786286705, 0.5688888889,
0.4786286705, 0.2369268850}}
};
}
public double Integrate(Func<double, double> f,
double a, double b, int n) {
if (!points.ContainsKey(n)) {
throw new ArgumentException(
"Number of points must be between 2 and 5"
);
}
double result = 0.0;
// Transform points and weights
for (int i = 0; i < n; i++) {
double x = (b-a)/2 * points[n][i] + (b+a)/2;
double w = (b-a)/2 * weights[n][i];
result += w * f(x);
}
return result;
}
public static void Main() {
// Example function: x²
double f(double x) => x * x;
var quad = new GaussianQuadrature();
try {
double result = quad.Integrate(f, -1.0, 1.0, 3);
Console.WriteLine($"∫x² dx from -1 to 1 ≈ {result:F6}");
Console.WriteLine($"Exact value: {2.0/3.0:F6}");
}
catch (Exception e) {
Console.WriteLine($"Error: {e.Message}");
}
}
}
