Integration

from scipy.integrate import (
    quad, dblquad, tplquad, nquad,
    fixed_quad, tanhsinh,
    solve_ivp, odeint,
    trapezoid, simpson, quad_vec,
    cumulative_trapezoid, cumulative_simpson
)
import numpy as np

from scipy.integrate import quad

# Basic usage: integrate f from a to b
val, err = quad(np.sin, 0, np.pi)
# val ≈ 2.0, err ≈ 2.2e-14 (absolute error estimate)

# Lambda
val, err = quad(lambda x: x**2, 0, 1)

# With arguments
def integrand(x, a, b):
    return a * np.exp(-b * x)

val, err = quad(integrand, 0, np.inf, args=(2, 1))

# Infinite limits
val, err = quad(lambda x: np.exp(-x**2), -np.inf, np.inf)  # = sqrt(pi)

# Key options
val, err = quad(f, a, b,
                limit=100,          # max number of subintervals
                epsabs=1.49e-8,     # absolute tolerance
                epsrel=1.49e-8,     # relative tolerance
                points=[0.5, 1.0],  # known singularities/discontinuities
                weight='sin',       # Fourier weights: 'cos', 'sin', 'alg', etc.
                wvar=omega,         # weight variable (e.g. frequency for sin/cos)
                full_output=True)   # more info dict

# Singularities: tell quad where they are
val, err = quad(lambda x: 1/np.sqrt(x), 0, 1, points=[0])  # sqrt singularity at 0

# Complex-valued integrand
from scipy.integrate import quad_vec
val, err = quad_vec(lambda x: np.exp(1j * x), 0, np.pi)

from scipy.integrate import quad_vec

# Integrand returns an array
def f_vec(x):
    return np.array([np.sin(x), np.cos(x), x**2])

val, err = quad_vec(f_vec, 0, np.pi)
# val.shape == (3,)

# Parametric family of integrals
def f_param(x, k):
    return np.sin(k * x)

# Vectorize over k
k_vals = np.arange(1, 6)
val, err = quad_vec(lambda x: np.sin(x * k_vals), 0, np.pi)  # shape (5,)

from scipy.integrate import dblquad, tplquad

# Double integral: int_{a}^{b} int_{gfun(x)}^{hfun(x)} f(y,x) dy dx
# Note the argument order: f(y, x) — inner variable first!
val, err = dblquad(lambda y, x: x*y, 0, 1,
                    lambda x: 0,          # y lower limit (can depend on x)
                    lambda x: x)          # y upper limit

# Constant inner limits
val, err = dblquad(lambda y, x: np.sin(x + y), 0, np.pi, 0, np.pi)

# Triple integral: f(z, y, x)
val, err = tplquad(lambda z, y, x: x*y*z,
                    0, 1,
                    0, 1,
                    0, 1)

from scipy.integrate import nquad

# f(x0, x1, ..., xn): last arg is outermost integration variable
def f(x, y, z):
    return x*y*z

# Fixed limits
val, err = nquad(f, [[0, 1], [0, 1], [0, 1]])

# Dynamic limits (as callables)
def y_lim(x):
    return [0, x]          # y in [0, x]

def z_lim(x, y):
    return [0, x + y]      # z in [0, x+y]

val, err = nquad(f, [z_lim, y_lim, [0, 1]])  # innermost limit first

from scipy.integrate import trapezoid, simpson, cumulative_trapezoid, cumulative_simpson

x = np.linspace(0, np.pi, 101)
y = np.sin(x)

# Trapezoidal rule
val = trapezoid(y, x)                 # ≈ 2.0
val = trapezoid(y, dx=x[1]-x[0])     # with uniform spacing
val = trapezoid(y, axis=0)            # along specific axis

# Simpson's rule (requires odd number of points, or last interval trapezoidal)
val = simpson(y, x=x)
val = simpson(y, dx=np.pi/100)

# Cumulative integrals (running integral)
F = cumulative_trapezoid(y, x)        # shape (n-1,): F[i] = int from x[0] to x[i+1]
F = cumulative_trapezoid(y, x, initial=0)  # shape (n,): starts at 0

F = cumulative_simpson(y, x=x)
F = cumulative_simpson(y, x=x, initial=0)

from scipy.integrate import fixed_quad, tanhsinh

# Gaussian quadrature (fixed n-point, non-adaptive)
val, _ = fixed_quad(np.sin, 0, np.pi, n=10)

# tanh-sinh (double exponential) quadrature — SciPy 1.15+
# Excellent for integrable endpoint singularities
res = tanhsinh(np.sin, 0, np.pi)
res.integral, res.error, res.success

from scipy.integrate import solve_ivp

# dy/dt = f(t, y), y(t0) = y0
def lotka_volterra(t, y, alpha, beta, delta, gamma):
    x, yd = y
    dxdt = alpha*x - beta*x*yd
    dydt = delta*x*yd - gamma*yd
    return [dxdt, dydt]

# Solve on t in [0, 15]
result = solve_ivp(
    fun=lotka_volterra,
    t_span=(0, 15),
    y0=[10, 5],          # initial conditions
    args=(1.5, 1.0, 0.75, 3.0),
    method='RK45',       # default
    t_eval=np.linspace(0, 15, 300),  # output times
    dense_output=True,   # also return solution callable
    rtol=1e-6,
    atol=1e-9
)

print(result.success, result.message)
print(result.t.shape, result.y.shape)  # (300,), (2, 300)

# Dense output: evaluate at any t
sol = result.sol(np.linspace(0, 15, 1000))   # shape (2, 1000)

# Stiff ODE example (use Radau or BDF)
def van_der_pol(t, y, mu):
    return [y[1], mu*(1 - y[0]**2)*y[1] - y[0]]

result = solve_ivp(van_der_pol, (0, 3000), [2, 0], args=(1000,),
                    method='Radau', rtol=1e-6, atol=1e-9)

def hit_ground(t, y):
    return y[0]      # zero when y[0] = 0

hit_ground.terminal = True    # stop integration
hit_ground.direction = -1     # only detect downward crossings (-1, 1, or 0)

def bounce(t, y):
    return y[0] - 1.0

bounce.terminal = False

result = solve_ivp(fun, t_span, y0, events=[hit_ground, bounce])
print(result.t_events[0])     # times of hit_ground events
print(result.y_events[0])     # state at those events

def jac(t, y, mu):
    return [[0, 1],
            [-2*mu*y[0]*y[1] - 1, mu*(1 - y[0]**2)]]

result = solve_ivp(van_der_pol, (0, 100), [2, 0], args=(1000,),
                    method='Radau', jac=jac, rtol=1e-8)

from scipy.integrate import odeint

# f(y, t, ...) — note t is second arg (opposite of solve_ivp)
def dydx(y, t, alpha):
    return -alpha * y

t = np.linspace(0, 10, 200)
y0 = [1.0]
sol = odeint(dydx, y0, t, args=(0.5,),
              rtol=1.49e-8, atol=1.49e-8,
              full_output=False)
# sol.shape = (200, 1)

from scipy.integrate import solve_bvp

# Example: y'' = -y, y(0) = 0, y(pi) = 0
# As first-order system: [y, y'] -> [y', -y]
def fun(x, y):
    return np.vstack([y[1], -y[0]])

def bc(ya, yb):
    return np.array([ya[0], yb[0]])   # y(0)=0, y(pi)=0

x_init = np.linspace(0, np.pi, 10)
y_init = np.zeros((2, 10))
y_init[0] = np.sin(x_init)   # initial guess

result = solve_bvp(fun, bc, x_init, y_init,
                    tol=1e-3, max_nodes=1000, verbose=0)

if result.status == 0:
    x_plot = np.linspace(0, np.pi, 200)
    y_plot = result.sol(x_plot)[0]   # y component only