Integration Fundamentals

import numpy as np
from scipy import integrate

# Define the function to integrate
def f(x):
    return x**2 * np.exp(-x)

# Compute the integral from 0 to infinity (improper integral)
result, error = integrate.quad(f, 0, np.inf)
print(f"Integral of x^2 * exp(-x) from 0 to inf:")
print(f"  Result: {result:.6f}")
print(f"  Error estimate: {error:.2e}")
print(f"  Exact (2! = 2): 2.000000")

# Definite integral from 0 to 1
result, error = integrate.quad(lambda x: np.sin(x), 0, np.pi)
print(f"\nIntegral of sin(x) from 0 to pi: {result:.6f} (exact: 2)")

import numpy as np
from scipy import integrate

# Define test function: f(x) = x^2
f = lambda x: x**2
a, b = 0, 1
exact = 1/3

# Trapezoidal rule
for n in [10, 100, 1000]:
    x = np.linspace(a, b, n + 1)
    trap_result = np.trapz(f(x), x)
    print(f"Trapezoidal (n={n:4d}): {trap_result:.8f}  error: {abs(trap_result - exact):.2e}")

print()

# Simpson's rule
for n in [10, 100, 1000]:
    x = np.linspace(a, b, n + 1)
    simp_result = integrate.simpson(f(x), x=x)
    print(f"Simpson's   (n={n:4d}): {simp_result:.8f}  error: {abs(simp_result - exact):.2e}")

print()

# Gaussian quadrature (scipy)
for n in [5, 10, 20]:
    result, error = integrate.fixed_quad(lambda x: x**2, a, b, n=n)
    print(f"Gauss quad  (n={n:4d}): {result:.8f}  error: {abs(result - exact):.2e}")

# scipy.integrate.quad (adaptive)
result, error = integrate.quad(f, a, b)
print(f"\nAdaptive quad:      {result:.8f}  error: {error:.2e}")

import sympy as sp

x = sp.Symbol('x')

# Symbolic indefinite integral
f = x**2 * sp.exp(x)
F = sp.integrate(f, x)
print(f"Indefinite integral of x^2 * e^x: {F}")

# Definite integral
result = sp.integrate(x**2, (x, 0, 1))
print(f"Definite integral of x^2 from 0 to 1: {result}")

# Improper integral
result = sp.integrate(sp.exp(-x**2), (x, -sp.oo, sp.oo))
print(f"Integral of e^(-x^2) from -inf to inf: {result}")
print(f"  = sqrt(pi) = {sp.sqrt(sp.pi)}")

# Verify Fundamental Theorem
F = sp.integrate(sp.sin(x), x)
print(f"\nAntiderivative of sin(x): {F}")
print(f"FTC: F(pi) - F(0) = {F.subs(x, sp.pi) - F.subs(x, 0)}")

import numpy as np

def monte_carlo_integrate(f, bounds, n_samples=100000):
    """Monte Carlo integration for arbitrary dimensions."""
    dim = len(bounds)
    samples = np.random.uniform(
        low=[b[0] for b in bounds],
        high=[b[1] for b in bounds],
        size=(n_samples, dim)
    )
    values = np.array([f(s) for s in samples])
    volume = np.prod([b[1] - b[0] for b in bounds])
    mean_val = np.mean(values)
    std_err = np.std(values) / np.sqrt(n_samples)
    return mean_val * volume, std_err

# Example: integral of x^2 + y^2 over [0,1] x [0,1]
f = lambda p: p[0]**2 + p[1]**2
result, error = monte_carlo_integrate(f, [(0, 1), (0, 1)])
exact = 2/3  # integral of x^2 + y^2 over [0,1]^2
print(f"Monte Carlo estimate: {result:.6f} +/- {error:.6f}")
print(f"Exact value:          {exact:.6f}")