A Short Python Tutorial[2]

import numpy as np
from scipy import integrate

# Define the integrand
def f(x):
    return np.sin(x) * np.exp(-x)

# Integration Interval
a, b = 0, np.pi

# 1. Using quad (based on adaptive Simpson and Gaussian quadrature methods)
result_quad, error_quad = integrate.quad(f, a, b)
print(f"quad result: {result_quad:.10f}, Estimated Error: {error_quad:.2e}")

# 2. Using fixed_quad (fixed-order Gaussian quadrature)
result_gauss, _ = integrate.fixed_quad(f, a, b, n=5)
print(f"5-Point Gaussian Quadrature Result: {result_gauss:.10f}")

# 3. Using the trapezoidal rule"
x_trap = np.linspace(a, b, 100) # generate 100 nodes
y_trap = f(x_trap)
result_trap = np.trapz(y_trap, x_trap)
print(f"Composite Trapezoidal Rule Result (100 points): {result_trap:.10f}")

# 4. Using Simpson's rule
result_simps = integrate.simpson(y_trap, x_trap)
print(f"Composite Simpson's Rule Result (100 points): {result_simps:.10f}")

# The exact value is provided for comparison (this example has an analytical solution)
exact = (1 + np.exp(-np.pi)) / 2
print(f"Exact Value: {exact:.10f}")

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
# Compare with Scipy's results
from scipy import integrate

def trapezoidal_rule(f, a, b, n):
    """
    Compute the approximate value of definite integral using composite trapezoidal rule

    Parameters:
    f: integrand function
    a: lower limit of integration
    b: upper limit of integration
    n: number of subintervals

    Returns:
    integral: approximate value of the integral
    """
    # Calculate step size
    h = (b - a) / n

    # Calculate x values of all nodes
    x = np.linspace(a, b, n + 1)

    # Calculate function values at all nodes
    y = f(x)

    # Apply trapezoidal rule formula
    # First term f(a) and last term f(b) have coefficient 1, 
    # middle terms have coefficient 2
    integral = (h / 2) * (y[0] + 2 * np.sum(y[1:n]) + y[n])

    return integral

# Test function
def f(x):
    return np.sin(x) + np.cos(x)

# Integration parameters
a, b = 0, np.pi  # Integration interval
n = 100          # Number of subintervals

# Using self-implemented function
my_result = trapezoidal_rule(f, a, b, n)

# Using Scipy's trapezoidal rule
x_scipy = np.linspace(a, b, n + 1)
y_scipy = f(x_scipy)
scipy_result = np.trapz(y_scipy, x_scipy)

# Using Scipy's quad function to obtain high-precision result
quad_result, error = integrate.quad(f, a, b)

print("Method comparison:")
print(f"Custom trapezoidal rule: {my_result:.10f}")
print(f"NumPy trapezoidal rule: {scipy_result:.10f}")
print(f"Scipy quad (high precision): {quad_result:.10f}")
print(f"Difference between two trapezoidal methods: {abs(my_result - scipy_result):.2e}")

import numpy as np

from scipy import integrate




def simpsons_rule(f, a, b, n):

    """Composite Simpson's rule implementation"""

    if n % 2 != 0:

        n += 1

    

    h = (b - a) / n

    x = np.linspace(a, b, n + 1)

    y = f(x)

    

    coefficients = np.ones(n + 1)

    coefficients[1:n:2] = 4

    coefficients[2:n-1:2] = 2

    

    return (h / 3) * np.dot(coefficients, y)




# Test function

def f(x):

    return np.sin(x) + np.cos(x)




# Integration parameters

a, b = 0, np.pi

n = 100  # Number of subintervals




# Compare various methods

my_result = simpsons_rule(f, a, b, n)




# Scipy's Simpson's rule

x_scipy = np.linspace(a, b, n + 1)

y_scipy = f(x_scipy)

scipy_result = integrate.simpson(y_scipy, x_scipy)




# Scipy high-precision integration

quad_result, error = integrate.quad(f, a, b)




# Trapezoidal rule comparison

trapz_result = np.trapz(y_scipy, x_scipy)




print("Method comparison:")

print(f"Custom Simpson's rule: {my_result:.10f}")

print(f"Scipy Simpson's rule: {scipy_result:.10f}")

print(f"Scipy quad (high precision): {quad_result:.10f}")

print(f"NumPy trapezoidal rule: {trapz_result:.10f}")

print()

print(f"Error of Simpson's rule vs exact value: {abs(my_result - quad_result):.2e}")

print(f"Error of trapezoidal rule vs exact value: {abs(trapz_result - quad_result):.2e}")