"""
Numerical Integration with SciPy
Author: RSK World
Website: https://rskworld.in
Email: help@rskworld.in
Phone: +91 93305 39277
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

# Example 1: Basic Definite Integration
def example_basic_integration():
    """
    Basic definite integration examples
    """
    print("=" * 60)
    print("Example 1: Basic Definite Integration")
    print("=" * 60)
    
    # Function: f(x) = x²
    def func(x):
        return x**2
    
    # Integrate from 0 to 2: ∫x²dx from 0 to 2 = x³/3|₀² = 8/3
    result, error = integrate.quad(func, 0, 2)
    exact = 8/3
    
    print(f"Function: f(x) = x²")
    print(f"Integration from 0 to 2")
    print(f"Numerical result: {result:.8f}")
    print(f"Exact result: {exact:.8f}")
    print(f"Absolute error: {abs(result - exact):.2e}")
    print(f"Estimated error: {error:.2e}")
    
    # Visualize
    x = np.linspace(0, 2, 100)
    y = func(x)
    
    plt.figure(figsize=(10, 6))
    plt.plot(x, y, 'b-', linewidth=2, label='f(x) = x²')
    plt.fill_between(x, 0, y, alpha=0.3, color='blue', label=f'Area = {result:.4f}')
    plt.xlabel('x', fontsize=12)
    plt.ylabel('f(x)', fontsize=12)
    plt.title('Definite Integration: ∫x²dx from 0 to 2', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig('integration_basic.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'integration_basic.png'")
    plt.close()


# Example 2: Integration with Infinite Limits
def example_infinite_integration():
    """
    Integration with infinite limits
    """
    print("\n" + "=" * 60)
    print("Example 2: Integration with Infinite Limits")
    print("=" * 60)
    
    # Gaussian function: f(x) = e^(-x²)
    def gaussian(x):
        return np.exp(-x**2)
    
    # Integrate from -∞ to +∞: ∫e^(-x²)dx = √π
    result, error = integrate.quad(gaussian, -np.inf, np.inf)
    exact = np.sqrt(np.pi)
    
    print(f"Function: f(x) = e^(-x²)")
    print(f"Integration from -∞ to +∞")
    print(f"Numerical result: {result:.8f}")
    print(f"Exact result: {exact:.8f}")
    print(f"Absolute error: {abs(result - exact):.2e}")
    
    # Visualize
    x = np.linspace(-4, 4, 200)
    y = gaussian(x)
    
    plt.figure(figsize=(10, 6))
    plt.plot(x, y, 'b-', linewidth=2, label='f(x) = e^(-x²)')
    plt.fill_between(x, 0, y, alpha=0.3, color='blue', label=f'Area = {result:.4f}')
    plt.xlabel('x', fontsize=12)
    plt.ylabel('f(x)', fontsize=12)
    plt.title('Gaussian Integration: ∫e^(-x²)dx from -∞ to +∞', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig('integration_infinite.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'integration_infinite.png'")
    plt.close()


# Example 3: Double Integration
def example_double_integration():
    """
    Double integration (multiple integrals)
    """
    print("\n" + "=" * 60)
    print("Example 3: Double Integration")
    print("=" * 60)
    
    # Function: f(x,y) = x² + y²
    def func(x, y):
        return x**2 + y**2
    
    # Integrate over rectangle: x from 0 to 1, y from 0 to 1
    # ∫∫(x² + y²)dxdy = ∫[x²y + y³/3]₀¹dx = ∫[x² + 1/3]dx = [x³/3 + x/3]₀¹ = 2/3
    result, error = integrate.dblquad(func, 0, 1, lambda x: 0, lambda x: 1)
    exact = 2/3
    
    print(f"Function: f(x,y) = x² + y²")
    print(f"Integration over [0,1] × [0,1]")
    print(f"Numerical result: {result:.8f}")
    print(f"Exact result: {exact:.8f}")
    print(f"Absolute error: {abs(result - exact):.2e}")
    
    # Visualize
    x = np.linspace(0, 1, 50)
    y = np.linspace(0, 1, 50)
    X, Y = np.meshgrid(x, y)
    Z = func(X, Y)
    
    fig = plt.figure(figsize=(12, 5))
    
    # Contour plot
    ax1 = plt.subplot(1, 2, 1)
    contour = plt.contour(X, Y, Z, levels=20, cmap='viridis')
    plt.colorbar(contour, label='f(x,y)')
    plt.xlabel('x', fontsize=12)
    plt.ylabel('y', fontsize=12)
    plt.title('Function: f(x,y) = x² + y²', fontsize=12, fontweight='bold')
    plt.grid(True, alpha=0.3)
    
    # 3D surface plot
    ax2 = plt.subplot(1, 2, 2, projection='3d')
    surf = ax2.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)
    ax2.set_xlabel('x', fontsize=10)
    ax2.set_ylabel('y', fontsize=10)
    ax2.set_zlabel('f(x,y)', fontsize=10)
    ax2.set_title('3D Surface Plot', fontsize=12, fontweight='bold')
    
    plt.tight_layout()
    plt.savefig('integration_double.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'integration_double.png'")
    plt.close()


# Example 4: Integration with Singularities
def example_singular_integration():
    """
    Integration with singularities
    """
    print("\n" + "=" * 60)
    print("Example 4: Integration with Singularities")
    print("=" * 60)
    
    # Function with singularity: f(x) = 1/√x
    def func(x):
        return 1 / np.sqrt(x)
    
    # Integrate from 0 to 1: ∫1/√x dx = 2√x|₀¹ = 2
    # Note: singularity at x=0
    result, error = integrate.quad(func, 0, 1, points=[0])
    exact = 2.0
    
    print(f"Function: f(x) = 1/√x")
    print(f"Integration from 0 to 1 (singularity at x=0)")
    print(f"Numerical result: {result:.8f}")
    print(f"Exact result: {exact:.8f}")
    print(f"Absolute error: {abs(result - exact):.2e}")
    
    # Visualize
    x = np.linspace(0.01, 1, 100)
    y = func(x)
    
    plt.figure(figsize=(10, 6))
    plt.plot(x, y, 'b-', linewidth=2, label='f(x) = 1/√x')
    plt.fill_between(x, 0, y, alpha=0.3, color='blue', label=f'Area = {result:.4f}')
    plt.xlabel('x', fontsize=12)
    plt.ylabel('f(x)', fontsize=12)
    plt.title('Integration with Singularity', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig('integration_singular.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'integration_singular.png'")
    plt.close()


# Example 5: Numerical Integration Methods Comparison
def example_integration_methods():
    """
    Compare different numerical integration methods
    """
    print("\n" + "=" * 60)
    print("Example 5: Integration Methods Comparison")
    print("=" * 60)
    
    # Function: f(x) = sin(x)
    def func(x):
        return np.sin(x)
    
    # Integration from 0 to π: ∫sin(x)dx = -cos(x)|₀^π = 2
    a, b = 0, np.pi
    exact = 2.0
    
    # Different methods
    methods = {
        'quad': integrate.quad(func, a, b),
        'fixed_quad': (integrate.fixed_quad(func, a, b, n=10)[0], 0),
        'trapezoid': (integrate.trapezoid(func(np.linspace(a, b, 100)), 
                                          np.linspace(a, b, 100)), 0),
        'simpson': (integrate.simpson(func(np.linspace(a, b, 100)), 
                                     np.linspace(a, b, 100)), 0),
    }
    
    print(f"Function: f(x) = sin(x)")
    print(f"Integration from 0 to π (exact = 2.0)")
    print(f"\n{'Method':<15} {'Result':<15} {'Error':<15}")
    print("-" * 45)
    
    for method, (result, error) in methods.items():
        err = abs(result - exact)
        print(f"{method:<15} {result:<15.8f} {err:<15.2e}")
    
    # Visualize
    x = np.linspace(0, np.pi, 100)
    y = func(x)
    
    plt.figure(figsize=(10, 6))
    plt.plot(x, y, 'b-', linewidth=2, label='f(x) = sin(x)')
    plt.fill_between(x, 0, y, alpha=0.3, color='blue', label=f'Area = {exact:.4f}')
    plt.xlabel('x', fontsize=12)
    plt.ylabel('f(x)', fontsize=12)
    plt.title('Integration Methods Comparison: ∫sin(x)dx from 0 to π', 
              fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig('integration_methods.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'integration_methods.png'")
    plt.close()


# Example 6: Cumulative Integration
def example_cumulative_integration():
    """
    Cumulative integration (antiderivative)
    """
    print("\n" + "=" * 60)
    print("Example 6: Cumulative Integration")
    print("=" * 60)
    
    # Function: f(x) = x²
    def func(x):
        return x**2
    
    # Generate data points
    x = np.linspace(0, 3, 100)
    y = func(x)
    
    # Cumulative integration
    y_cumulative = integrate.cumulative_trapezoid(y, x, initial=0)
    
    # Exact antiderivative: F(x) = x³/3
    y_exact = x**3 / 3
    
    print(f"Function: f(x) = x²")
    print(f"Cumulative integration (antiderivative): F(x) = x³/3")
    print(f"At x = 3: F(3) = {y_cumulative[-1]:.6f} (exact: {y_exact[-1]:.6f})")
    
    # Visualize
    plt.figure(figsize=(12, 5))
    
    plt.subplot(1, 2, 1)
    plt.plot(x, y, 'b-', linewidth=2, label='f(x) = x²')
    plt.xlabel('x', fontsize=12)
    plt.ylabel('f(x)', fontsize=12)
    plt.title('Original Function', fontsize=12, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.subplot(1, 2, 2)
    plt.plot(x, y_cumulative, 'r-', linewidth=2, label='Numerical F(x)')
    plt.plot(x, y_exact, 'g--', linewidth=2, label='Exact F(x) = x³/3')
    plt.xlabel('x', fontsize=12)
    plt.ylabel('F(x)', fontsize=12)
    plt.title('Cumulative Integration (Antiderivative)', fontsize=12, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('integration_cumulative.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'integration_cumulative.png'")
    plt.close()


def main():
    """
    Main function to run all integration examples
    """
    print("\n" + "=" * 60)
    print("SciPy Numerical Integration Examples")
    print("Author: RSK World - https://rskworld.in")
    print("=" * 60 + "\n")
    
    example_basic_integration()
    example_infinite_integration()
    example_double_integration()
    example_singular_integration()
    example_integration_methods()
    example_cumulative_integration()
    
    print("\n" + "=" * 60)
    print("All integration examples completed!")
    print("=" * 60)


if __name__ == "__main__":
    main()