Help us shape the future of Go! Take the 2026 Go Developer Survey

integrate

package module
v0.2.0 Latest Latest
Warning

This package is not in the latest version of its module.

 Go to latest Published: Oct 17, 2026 License: BSD-3-Clause Imports: 6 Imported by: 0

README

go-numerics/integrate

integrate

Docs pkg.go.dev License

Numerical integration and ODE solvers in pure Go (CGO_ENABLED=0): adaptive quadrature (QUADPACK), fixed and sample rules, multiple integrals, and initial value problems with RK45, RK23, DOP853, Radau, BDF and automatic stiffness switching — the functionality of scipy.integrate for the Go scientific stack (go-ndarray, go-plot, go-notebook), with the same algorithms, defaults and diagnostics.

🌐 Website · 📚 Documentation

go get github.com/go-numerics/integrate

Headline numbers

  • Same algorithms, same paths. QUADPACK is ported from the netlib Fortran; the ODE methods follow SciPy's implementations step for step. Against SciPy 1.18.1, 116 of 124 fixture integrals take exactly SciPy's subdivision path (same evaluation count, value to 1e-12), and 110 of 120 solve_ivp runs of the ported methods take exactly SciPy's steps (same nfev, steps, njev, nlu). The remaining differences are explained in docs/verification.md: roundoff-limited paths that depend on the last bit, and three integrals where SciPy's C translation of QUADPACK departs from the Fortran and this port agrees with the Fortran.
  • Faster. 6–8× SciPy on quad (3× without counting Python callbacks), 90–190× on solve_ivp (RK45, DOP853, Radau, BDF; same step sequences), 8–16× on odeint/LSODA (docs/perf.md).
  • Verified. Fixtures from SciPy and mpmath, the netlib Fortran as a second oracle, a negative control on the coefficients, 100 % statement coverage with -race, CI lanes on amd64, arm64, riscv64, loong64, ppc64le, s390x, 386 and arm.

From SciPy to Go

SciPy Go
quad(f, a, b) integrate.Quad(f, a, b, nil)
quad(f, 0, np.inf) integrate.Quad(f, 0, math.Inf(1), nil)
quad(f, a, b, epsabs=0, epsrel=1e-12, limit=200) Quad(f, a, b, &integrate.QuadOptions{EpsAbs: -1, EpsRel: 1e-12, Limit: 200})
quad(f, a, b, points=[c]) Quad(f, a, b, &integrate.QuadOptions{Points: []float64{c}})
quad(f, a, b, weight="cos", wvar=w) (also "sin", infinite b) Quad(f, a, b, &integrate.QuadOptions{Weight: integrate.WeightCos, Omega: w})
quad(f, a, b, weight="alg", wvar=(al, be)) ("alg-loga", "alg-logb", "alg-log") Quad(f, a, b, &integrate.QuadOptions{Weight: integrate.WeightAlg, Alpha: al, Beta: be})
quad(f, a, b, weight="cauchy", wvar=c) Quad(f, a, b, &integrate.QuadOptions{Weight: integrate.WeightCauchy, C: c})
IntegrationWarning (ier 1–5, 7) *integrate.QuadError returned with the result; errors.Is(err, integrate.ErrRoundoff) …
dblquad(f, a, b, gfun, hfun) integrate.DblQuad(f, a, b, gfun, hfun, nil), constants as integrate.Const(v)
tplquad(f, a, b, g, h, q, r) integrate.TplQuad(f, a, b, g, h, q, r, nil)
nquad(f, ranges, opts=...) integrate.NQuad(f, ranges, &integrate.NQuadOptions{...})
fixed_quad(f, a, b, n=5) integrate.FixedQuad(f, a, b, 5)
romberg(f, a, b) (SciPy < 1.15) integrate.Romberg(f, a, b, nil)
trapezoid(y, x), simpson(y, x=x), cumulative_trapezoid(y, x, initial=0), romb(y, dx) Trapezoid(y, x, 0), Simpson(y, x, 0), CumulativeTrapezoid(y, x, 0, true), Romb(y, dx); TrapezoidND, SimpsonND, CumulativeTrapezoidND take an *ndarray.Array and an axis
solve_ivp(fun, (t0, tf), y0) integrate.SolveIVP(f, [2]float64{t0, tf}, y0, nil)
method="Radau", rtol=1e-6, atol=[...], max_step, first_step, jac &integrate.IVPOptions{Method: integrate.Radau, RTol: 1e-6, ATolVec: []float64{...}, MaxStep: …, FirstStep: …, Jac: …}
t_eval=ts, dense_output=True, sol.sol(t) TEval: ts, DenseOutput: true, sol.Sol.Eval(t)
events=[g] with g.terminal = True, g.direction = -1 Events: []integrate.Event{{F: g, Terminal: 1, Direction: -1}}
sol.t, sol.y, sol.t_events, sol.nfev sol.T, sol.Y (sol.YArray() for SciPy's n × len(t) ndarray), sol.TEvents, sol.NFev
odeint(func, y0, t, tfirst=True) integrate.Odeint(f, y0, t, nil) → len(t) × n *ndarray.Array and an info record

A right-hand side writes its derivative instead of returning it:

// SciPy: def lotka(t, y): return [1.5*y[0] - y[0]*y[1], -3*y[1] + y[0]*y[1]]
lotka := func(t float64, y, dydt []float64) {
	dydt[0] = 1.5*y[0] - y[0]*y[1]
	dydt[1] = -3*y[1] + y[0]*y[1]
}
sol, err := integrate.SolveIVP(lotka, [2]float64{0, 15}, []float64{10, 5}, nil)
// sol.NFev == 470, len(sol.T) == 57 — the same as SciPy's solve_ivp.

A stiff problem with an event:

const mu = 1000
vdp := func(t float64, y, dydt []float64) {
	dydt[0] = y[1]
	dydt[1] = mu*(1-y[0]*y[0])*y[1] - y[0]
}
cross := integrate.Event{F: func(t float64, y []float64) float64 { return y[0] }}
sol, err := integrate.SolveIVP(vdp, [2]float64{0, 3000}, []float64{2, 0},
	&integrate.IVPOptions{Method: integrate.Radau, Events: []integrate.Event{cross}, DenseOutput: true})

Quadrature with diagnostics:

r, err := integrate.Quad(func(x float64) float64 { return math.Log(x) / math.Sqrt(x) }, 0, 1, nil)
// r.Value = -4 (to 1.7e-14), r.AbsErr = 2.3e-13, r.NEval = 315, err == nil
_, err = integrate.Quad(func(x float64) float64 { return 1 / x }, 0, 1, nil)
// errors.Is(err, integrate.ErrLimit)

More in example_test.go.

Documentation

Not provided

Complex-valued integrands and states, vectorized right-hand sides, sparse or banded Jacobians, quad_vec, cubature, tanhsinh, solve_bvp and the legacy ode class. LSODA is an automatic Runge–Kutta/BDF switch, not a port of ODEPACK's Adams/BDF LSODA: same purpose, different evaluation counts (see design.md).

License

BSD 3-Clause; see LICENSE.

Documentation

Overview

Package integrate is numerical integration and ODE solving in pure Go (CGO_ENABLED=0): the functionality of scipy.integrate, with the same algorithms, tolerances and diagnostics, so that a notebook written against SciPy finds the same answers here.

Quadrature:

  • Quad: QUADPACK's adaptive Gauss–Kronrod drivers, ported from the netlib Fortran: QAGS (21-point rule with Wynn's epsilon extrapolation), QAGI (infinite ranges), QAGP (break points), QAWO and QAWF (cos/sin weights, Fourier integrals), QAWC (Cauchy principal values) and QAWS (algebraic-logarithmic end-point singularities). QUADPACK's ier diagnostics are returned as *QuadError values that wrap ErrLimit, ErrRoundoff, ErrBadIntegrand, ErrExtrapolation, ErrDivergent, ErrInvalidInput and ErrCycles.
  • DblQuad, TplQuad and NQuad: iterated adaptive quadrature over regions whose inner limits may depend on the outer variables.
  • FixedQuad (Gauss–Legendre), Romberg, Romb, Simpson, Trapezoid and CumulativeTrapezoid on functions or samples, with ndarray variants integrating along an axis (TrapezoidND, SimpsonND, CumulativeTrapezoidND).

Initial value problems:

  • SolveIVP with RK45 (Dormand–Prince), RK23 (Bogacki–Shampine), DOP853 (Hairer), Radau (implicit Radau IIA, order 5), BDF (variable order 1–5, NDF variant) and LSODA (automatic switching between Runge–Kutta and BDF). Events with direction and terminal counts, output at chosen times (TEval), continuous dense output, scalar or per-component tolerances, MaxStep, FirstStep and optional Jacobians (finite differences otherwise).
  • Odeint: SciPy's classic interface over the switching method.

Right-hand sides are func(t float64, y, dydt []float64), writing the derivative into dydt so that a step allocates nothing; SciPy's fun(t, y) -> dydt maps one to one. Each exported function documents the SciPy call it corresponds to. docs/design.md describes the algorithms and their sources, docs/verification.md how they were checked against SciPy, the netlib Fortran and mpmath, docs/perf.md the timings.

Index

Examples

Constants

View Source 
const (
	DefaultEpsAbs = 1.49e-8
	DefaultEpsRel = 1.49e-8
	DefaultLimit  = 50
	DefaultMaxP1  = 50
	DefaultLimLst = 50
)

Default tolerances and limits of Quad, the same as SciPy's quad.

Variables

View Source 
var (
	// ErrStepTooSmall: the required step size fell below the spacing of
	// floating-point numbers at t.
	ErrStepTooSmall = errors.New("integrate: required step size is less than spacing between numbers")
	// ErrODEInput: the problem or the options are invalid.
	ErrODEInput = errors.New("integrate: invalid ODE input")
	// ErrExcessWork: Odeint took more than MxStep steps between two
	// output times.
	ErrExcessWork = errors.New("integrate: excess work done on this call")
)

Errors returned by the ODE solvers.

View Source 
var (
	// ErrLimit: the maximum number of subdivisions (or QAWF cycles) was
	// reached (ier 1).
	ErrLimit = errors.New("integrate: maximum number of subdivisions reached")
	// ErrRoundoff: roundoff error prevents the requested tolerance from
	// being achieved (ier 2).
	ErrRoundoff = errors.New("integrate: roundoff error detected")
	// ErrBadIntegrand: extremely bad integrand behaviour at some points
	// (ier 3).
	ErrBadIntegrand = errors.New("integrate: extremely bad integrand behaviour")
	// ErrExtrapolation: the extrapolation does not converge because of
	// roundoff (ier 4).
	ErrExtrapolation = errors.New("integrate: extrapolation does not converge")
	// ErrDivergent: the integral is probably divergent or slowly
	// convergent (ier 5).
	ErrDivergent = errors.New("integrate: integral is probably divergent")
	// ErrInvalidInput: the input is invalid (ier 6); no estimate is
	// returned.
	ErrInvalidInput = errors.New("integrate: invalid input")
	// ErrCycles: abnormal termination in one or more QAWF cycles (ier 7).
	ErrCycles = errors.New("integrate: bad integrand behaviour in one or more cycles")
)

Sentinel errors, one per QUADPACK diagnostic. A *QuadError wraps one of them, so errors.Is(err, integrate.ErrRoundoff) works.

View Source 
var ErrDivMax = errors.New("integrate: Romberg: DivMax exceeded")

ErrDivMax reports that Romberg reached DivMax before converging.

View Source 
var ErrSamples = errors.New("integrate: invalid samples")

ErrSamples reports sample arrays that do not fit a fixed rule.

Functions

func Const 

func Const(v float64) func(float64) float64

Const returns the function x ↦ v, for constant limits of DblQuad and TplQuad.

func Const2 

func Const2(v float64) func(float64, float64) float64

Const2 returns the function (x, y) ↦ v, for constant inner limits of TplQuad.

func CumulativeTrapezoid 

func CumulativeTrapezoid(y, x []float64, dx float64, initial bool) ([]float64, error)

CumulativeTrapezoid returns the running trapezoidal integral of y. Without initial the result has len(y)-1 entries; with initial a leading 0 is prepended (SciPy's initial=0). It is SciPy's cumulative_trapezoid(y, x=None, dx=1.0, initial=None).

func CumulativeTrapezoidND 

func CumulativeTrapezoidND(y *ndarray.Array, x []float64, dx float64, axis int, initial bool) (*ndarray.Array, error)

CumulativeTrapezoidND is SciPy's cumulative_trapezoid(y, x, dx, axis, initial=0 if initial else None) for an ndarray.

func FixedQuad 

func FixedQuad(f func(float64) float64, a, b float64, n int) (float64, error)

FixedQuad integrates f over [a, b] with the n-point Gauss–Legendre rule. It is SciPy's fixed_quad(func, a, b, n=5) (whose second return value is always None).

func GaussLegendre 

func GaussLegendre(n int) (x, w []float64)

GaussLegendre returns the n nodes (ascending) and weights of the Gauss–Legendre rule on [-1, 1] (SciPy's roots_legendre). Nodes are found by Newton's method on the three-term recurrence, started from Tricomi's asymptotic approximation; results are cached and must not be modified.

func Romb 

func Romb(y []float64, dx float64) (float64, error)

Romb integrates 2^k + 1 equally spaced samples by Romberg extrapolation. It is SciPy's romb(y, dx).

func Romberg 

func Romberg(f func(float64) float64, a, b float64, opts *RombergOptions) (float64, error)

Romberg integrates f over [a, b] by Richardson extrapolation of the trapezoid rule on 1, 2, 4, ... 2^DivMax intervals. It is the scipy.integrate.romberg of SciPy before 1.15, which removed it. If DivMax is reached first the last estimate is returned with ErrDivMax.

func Simpson 

func Simpson(y, x []float64, dx float64) (float64, error)

Simpson integrates samples y with the composite Simpson rule at abscissae x (or spacing dx when x is nil). For an even number of samples the last interval uses Cartwright's correction, as SciPy ≥ 1.11 does. It is SciPy's simpson(y, x=None, dx=1.0).

Example

simpson(y, x=x) 

package main

import (
	"fmt"

	"github.com/go-numerics/integrate"
)

func main() {
	x := []float64{0, 0.5, 1, 1.5, 2}
	y := make([]float64, len(x))
	for i, v := range x {
		y[i] = v * v * v
	}
	v, _ := integrate.Simpson(y, x, 0)
	fmt.Println(v)
}

Output:
4

func SimpsonND 

func SimpsonND(y *ndarray.Array, x []float64, dx float64, axis int) (*ndarray.Array, error)

SimpsonND integrates y along axis with Simpson's rule (SciPy's simpson(y, x=x, dx=dx, axis=axis)).

func Trapezoid 

func Trapezoid(y, x []float64, dx float64) (float64, error)

Trapezoid integrates samples y with the composite trapezoidal rule, at abscissae x (or spacing dx when x is nil). It is SciPy's trapezoid(y, x=None, dx=1.0).

func TrapezoidND 

func TrapezoidND(y *ndarray.Array, x []float64, dx float64, axis int) (*ndarray.Array, error)

TrapezoidND integrates y along axis with the trapezoidal rule, at abscissae x (nil: spacing dx). It is SciPy's trapezoid(y, x, dx, axis) for an ndarray; the integrated axis is removed from the result.

Types

type Event 

type Event struct {
	// F is the event function g(t, y); events are its zero crossings.
	F func(t float64, y []float64) float64
	// Terminal is 0 for an event that is only recorded; k > 0 stops the
	// integration at the k-th occurrence (SciPy's terminal=True is 1).
	Terminal int
	// Direction > 0 keeps only crossings from negative to positive,
	// < 0 only from positive to negative, 0 both.
	Direction float64
}

Event describes a function whose zeros SolveIVP locates: SciPy's event callables with their terminal and direction attributes.

type IVPOptions 

type IVPOptions struct {
	// Method is the integrator (default RK45).
	Method Method
	// TEval, if non-nil, lists the times (sorted in the direction of
	// integration, inside the span) at which the solution is stored;
	// otherwise every step is stored.
	TEval []float64
	// DenseOutput requests a continuous solution (Solution.Sol).
	DenseOutput bool
	// Events are located by root finding on the dense output.
	Events []Event
	// RTol and ATol are the relative and absolute tolerances (defaults
	// 1e-3 and 1e-6). The local error estimate e of each component is kept
	// below ATol + RTol·|y|.
	RTol, ATol float64
	// RTolVec and ATolVec give per-component tolerances; when non-nil they
	// replace RTol and ATol.
	RTolVec, ATolVec []float64
	// MaxStep bounds the step size (default +Inf).
	MaxStep float64
	// FirstStep is the initial step size; 0 lets the solver choose.
	FirstStep float64
	// Jac is the Jacobian for Radau, BDF and LSODA; nil uses finite
	// differences.
	Jac JacFunc
}

IVPOptions tunes SolveIVP. The zero value gives SciPy's defaults.

type JacFunc 

type JacFunc func(t float64, y []float64, jac []float64)

JacFunc evaluates the Jacobian ∂f_i/∂y_j of an ODEFunc into jac, an n×n row-major matrix: jac[i*n+j] = ∂f_i/∂y_j. It is SciPy's jac argument.

type Method 

type Method int

Method selects the integrator of SolveIVP. 

const (
	// RK45 is the explicit Dormand–Prince 5(4) pair with its 4th-order
	// free interpolant (SciPy's default, method="RK45").
	RK45 Method = iota
	// RK23 is the explicit Bogacki–Shampine 3(2) pair ("RK23").
	RK23
	// DOP853 is Hairer's explicit Dormand–Prince 8(5,3) method with its
	// 7th-order interpolant ("DOP853").
	DOP853
	// Radau is the implicit Radau IIA method of order 5 ("Radau").
	Radau
	// BDF is the implicit variable-order (1–5) backward differentiation
	// method in its numerical differentiation formula (NDF) variant, the
	// quasi-constant step size implementation of SciPy's "BDF".
	BDF
	// LSODA switches automatically between an explicit Runge–Kutta method
	// on non-stiff stretches and BDF on stiff ones ("LSODA" in spirit; see
	// docs/design.md for how it differs from ODEPACK's Adams/BDF LSODA).
	LSODA
)

func (Method) String 

func (m Method) String() string

type NQuadOptions 

type NQuadOptions struct {
	// Opts gives the QuadOptions of each level, innermost first (index i
	// applies to x[i]). A single entry applies to every level; nil uses
	// the defaults. Points of a level are interpreted in that level's
	// variable.
	Opts []*QuadOptions
}

NQuadOptions tunes NQuad.

type ODEFunc 

type ODEFunc func(t float64, y, dydt []float64)

ODEFunc evaluates the right-hand side of dy/dt = f(t, y), writing f(t, y) into dydt. It must not retain y or dydt. SciPy's fun(t, y) returns the derivative; here it is written into a caller-owned slice so that a step allocates nothing.

type ODESolution 

type ODESolution struct {
	// contains filtered or unexported fields
}

ODESolution is a continuous solution assembled from the local interpolants of the steps (SciPy's OdeSolution).

func (*ODESolution) Eval 

func (s *ODESolution) Eval(t float64) []float64

Eval returns the solution at t.

func (*ODESolution) EvalArray 

func (s *ODESolution) EvalArray(ts []float64) *ndarray.Array

EvalArray evaluates the solution at the times ts and returns an n×len(ts) ndarray, as SciPy's sol(ts) does.

func (*ODESolution) EvalInto 

func (s *ODESolution) EvalInto(t float64, y []float64)

EvalInto writes the solution at t into y. Outside [TMin, TMax] the end interpolants extrapolate.

func (*ODESolution) TMax 

func (s *ODESolution) TMax() float64

TMax is the upper end of the interval covered.

func (*ODESolution) TMin 

func (s *ODESolution) TMin() float64

TMin and TMax are the ends of the interval covered.

type OdeintInfo 

type OdeintInfo struct {
	// HU is the step size last used, TCur the time reached.
	HU, TCur []float64
	// NST, NFE and NJE are the cumulative numbers of steps, right-hand side
	// evaluations and Jacobian evaluations.
	NST, NFE, NJE []int
	// MUsed is the method of the last step: 1 non-stiff (Runge–Kutta),
	// 2 stiff (BDF).
	MUsed []int
	// Message describes the outcome.
	Message string
}

OdeintInfo is SciPy's odeint infodict (full_output=True): one entry per output time after the first.

func Odeint 

func Odeint(f ODEFunc, y0 []float64, t []float64, opts *OdeintOptions) (*ndarray.Array, *OdeintInfo, error)

Odeint integrates dy/dt = f(t, y) from t[0] and returns y at every time of t (which must be monotonic) as a len(t)×n ndarray whose first row is y0. It is SciPy's odeint(func, y0, t, Dfun=...) over the LSODA-style switching method; SciPy's func(y, t) takes the state first, here f(t, y, dydt) takes time first, as with odeint(..., tfirst=True). On failure the rows not reached are NaN and the error wraps ErrExcessWork or ErrStepTooSmall.

Example

odeint(robertson, [1, 0, 0], t, tfirst=True) 

package main

import (
	"fmt"

	"github.com/go-numerics/integrate"
)

func main() {
	robertson := func(t float64, y, dydt []float64) {
		dydt[0] = -0.04*y[0] + 1e4*y[1]*y[2]
		dydt[1] = 0.04*y[0] - 1e4*y[1]*y[2] - 3e7*y[1]*y[1]
		dydt[2] = 3e7 * y[1] * y[1]
	}
	y, info, _ := integrate.Odeint(robertson, []float64{1, 0, 0}, []float64{0, 40, 4e4}, nil)
	fmt.Printf("%.4f %.4f %d\n", y.At(1, 0), y.At(2, 2), info.MUsed[1])
}

Output:
0.7158 0.9610 2

type OdeintOptions 

type OdeintOptions struct {
	// RTol and ATol default to 1.49012e-8. A negative ATol means 0.
	RTol, ATol float64
	// RTolVec and ATolVec give per-component tolerances.
	RTolVec, ATolVec []float64
	// Jac is the Jacobian (SciPy's Dfun, with col_deriv=False).
	Jac JacFunc
	// H0 is the first step (0: automatic) and HMax the largest step
	// (0: unbounded).
	H0, HMax float64
	// MxStep bounds the number of steps between two output times
	// (default 500).
	MxStep int
}

OdeintOptions tunes Odeint. Zero values select SciPy's odeint defaults.

type QuadError 

type QuadError struct {
	// Ier is QUADPACK's error code.
	Ier int
	// Detail says which routine or which argument is concerned.
	Detail string
}

QuadError reports a QUADPACK diagnostic. For ier 1–5 and 7 the QuadResult returned alongside still carries the best estimate available, as SciPy returns it with an IntegrationWarning.

Example

quad(f, 0, 1) on 1/x reports a diagnostic, as SciPy's IntegrationWarning. 

package main

import (
	"fmt"

	"github.com/go-numerics/integrate"
)

func main() {
	_, err := integrate.Quad(func(x float64) float64 { return 1 / x }, 0, 1, nil)
	fmt.Println(err)
}

Output:
integrate: maximum number of subdivisions reached

func (*QuadError) Error 

func (e *QuadError) Error() string

func (*QuadError) Unwrap 

func (e *QuadError) Unwrap() error

Unwrap returns the sentinel error of the diagnostic.

type QuadOptions 

type QuadOptions struct {
	// EpsAbs and EpsRel are the absolute and relative error targets: Quad
	// tries to make |Value - I| ≤ max(EpsAbs, EpsRel·|I|). Zero selects
	// 1.49e-8. A negative value switches that criterion off, which is
	// SciPy's epsabs=0 (or epsrel=0).
	EpsAbs, EpsRel float64
	// Limit bounds the number of subintervals (default 50).
	Limit int
	// Points are interior break points where the integrand is known to be
	// difficult (discontinuities, kinks, integrable singularities). Only
	// for finite ranges without a weight (QAGP). Points outside (a, b) are
	// ignored, as SciPy does.
	Points []float64
	// Weight selects a weight function; see Weight.
	Weight Weight
	// Omega is the angular frequency of WeightCos and WeightSin.
	Omega float64
	// Alpha and Beta are the exponents of the algebraic weights (> -1).
	Alpha, Beta float64
	// C is the pole of WeightCauchy, strictly inside (a, b).
	C float64
	// MaxP1 bounds the number of Chebyshev moment sets kept by QAWO
	// (default 50).
	MaxP1 int
	// LimLst bounds the number of cycles of QAWF (default 50, at least 3).
	LimLst int
}

QuadOptions tunes Quad. The zero value (or a nil pointer) gives SciPy's defaults.

type QuadResult 

type QuadResult struct {
	// Value is the estimate of the integral.
	Value float64
	// AbsErr is the estimated absolute error.
	AbsErr float64
	// NEval is the number of integrand evaluations.
	NEval int
	// NIntervals is the number of subintervals the adaptive routine used
	// (QUADPACK's last); for QAWF it is the number of cycles.
	NIntervals int
	// CycleIer holds the QUADPACK ier of each QAWF cycle (SciPy's
	// infodict["ierlst"]); nil for other routines.
	CycleIer []int
}

QuadResult is the outcome of a quadrature.

func DblQuad 

func DblQuad(f func(y, x float64) float64, a, b float64, gfun, hfun func(x float64) float64, opts *QuadOptions) (QuadResult, error)

DblQuad computes the double integral of f(y, x) for x from a to b and y from gfun(x) to hfun(x). It is SciPy's dblquad(func, a, b, gfun, hfun, epsabs=..., epsrel=...), with constant limits written Const(v): 

dblquad(f, 0, 1, 0, 2)                → DblQuad(f, 0, 1, Const(0), Const(2), nil)
dblquad(f, 0, 1, lambda x: 0, lambda x: x) → DblQuad(f, 0, 1, Const(0), func(x float64) float64 { return x }, nil)

opts (tolerances, Limit) apply to both levels, as in SciPy.

Example

dblquad(lambda y, x: x*y*y, 0, 2, 0, 1) 

package main

import (
	"fmt"

	"github.com/go-numerics/integrate"
)

func main() {
	r, _ := integrate.DblQuad(func(y, x float64) float64 { return x * y * y }, 0, 2,
		integrate.Const(0), integrate.Const(1), nil)
	fmt.Printf("%.12f\n", r.Value)
}

Output:
0.666666666667

func NQuad 

func NQuad(f func(x []float64) float64, ranges []Range, opts *NQuadOptions) (QuadResult, error)

NQuad integrates f over a hyper-rectangle or a more general region whose inner limits depend on the outer variables. ranges[i] gives the limits of x[i]; x[0] is integrated innermost. It is SciPy's scipy.integrate.nquad(func, ranges, opts=opts): 

nquad(f, [[0, 1], [0, 2]])              → NQuad(f, []Range{ConstRange(0, 1), ConstRange(0, 2)}, nil)
nquad(f, [lambda x1: (0, x1), [0, 1]])  → NQuad(f, []Range{func(o []float64) (float64, float64) { return 0, o[0] }, ConstRange(0, 1)}, nil)

f receives x with len(ranges) entries (SciPy passes them as separate arguments). AbsErr is the largest error estimate of any one-dimensional quadrature, as in SciPy; NEval counts the evaluations of f. The first QUADPACK diagnostic met, if any, is returned with the result.

func Quad 

func Quad(f func(float64) float64, a, b float64, opts *QuadOptions) (QuadResult, error)

Quad computes the definite integral of f from a to b, where a and b may be infinite. It is SciPy's scipy.integrate.quad: 

quad(f, a, b)                                   → Quad(f, a, b, nil)
quad(f, a, b, epsabs=0, epsrel=1e-12, limit=200) → Quad(f, a, b, &QuadOptions{EpsAbs: -1, EpsRel: 1e-12, Limit: 200})
quad(f, a, b, points=[c])                       → Quad(f, a, b, &QuadOptions{Points: []float64{c}})
quad(f, a, b, weight="cos", wvar=w)             → Quad(f, a, b, &QuadOptions{Weight: WeightCos, Omega: w})
quad(f, a, b, weight="alg", wvar=(al, be))      → Quad(f, a, b, &QuadOptions{Weight: WeightAlg, Alpha: al, Beta: be})
quad(f, a, b, weight="cauchy", wvar=c)          → Quad(f, a, b, &QuadOptions{Weight: WeightCauchy, C: c})

The routine is chosen as SciPy chooses it: QAGS on finite ranges, QAGI on infinite ones, QAGP with break points, QAWO/QAWF/QAWS/QAWC with a weight. When QUADPACK reports a problem the result is returned together with a *QuadError; for invalid input (ier 6) the result is zero.

Example

quad(lambda x: math.log(x) / math.sqrt(x), 0, 1) 

package main

import (
	"fmt"
	"math"

	"github.com/go-numerics/integrate"
)

func main() {
	r, err := integrate.Quad(func(x float64) float64 { return math.Log(x) / math.Sqrt(x) }, 0, 1, nil)
	fmt.Printf("%.12f %v %d\n", r.Value, err, r.NEval)
}

Output:
-4.000000000000 <nil> 315
Example (Fourier)

quad(lambda x: math.exp(-x), 0, np.inf, weight="cos", wvar=2) 

package main

import (
	"fmt"
	"math"

	"github.com/go-numerics/integrate"
)

func main() {
	r, _ := integrate.Quad(func(x float64) float64 { return math.Exp(-x) }, 0, math.Inf(1),
		&integrate.QuadOptions{Weight: integrate.WeightCos, Omega: 2})
	fmt.Printf("%.12f\n", r.Value)
}

Output:
0.200000000000
Example (Infinite)

quad(lambda x: math.exp(-x*x), -np.inf, np.inf) 

package main

import (
	"fmt"
	"math"

	"github.com/go-numerics/integrate"
)

func main() {
	r, _ := integrate.Quad(func(x float64) float64 { return math.Exp(-x * x) }, math.Inf(-1), math.Inf(1), nil)
	fmt.Printf("%.15f\n", r.Value)
}

Output:
1.772453850905516

func TplQuad 

func TplQuad(f func(z, y, x float64) float64, a, b float64, gfun, hfun func(x float64) float64, qfun, rfun func(x, y float64) float64, opts *QuadOptions) (QuadResult, error)

TplQuad computes the triple integral of f(z, y, x) for x in [a, b], y in [gfun(x), hfun(x)] and z in [qfun(x, y), rfun(x, y)]. It is SciPy's tplquad(func, a, b, gfun, hfun, qfun, rfun).

type Range 

type Range func(outer []float64) (lo, hi float64)

Range gives the integration limits of one variable of NQuad as a function of the outer variables: outer[0] is the variable integrated just outside this one, outer[1] the next, and so on (the arguments SciPy's nquad passes to a callable range).

func ConstRange 

func ConstRange(lo, hi float64) Range

ConstRange is a Range with fixed limits.

type RombergOptions 

type RombergOptions struct {
	Tol    float64 // absolute tolerance (default 1.48e-8)
	RTol   float64 // relative tolerance (default 1.48e-8)
	DivMax int     // maximum order of extrapolation (default 10)
}

RombergOptions tunes Romberg; zero values select SciPy's defaults.

type Solution 

type Solution struct {
	// T holds the time points.
	T []float64
	// Y holds the states: Y[k] is y(T[k]). SciPy's sol.y is the transpose
	// (component first); YArray returns that layout.
	Y [][]float64
	// Sol is the continuous solution when DenseOutput was requested.
	Sol *ODESolution
	// TEvents[i] and YEvents[i] hold the times and states at which event i
	// occurred.
	TEvents [][]float64
	YEvents [][][]float64
	// NFev, NJev and NLU count right-hand side evaluations, Jacobian
	// evaluations and LU decompositions, as SciPy counts them: NFev does
	// not include the evaluations spent on finite-difference Jacobians,
	// which NFevJac reports.
	NFev, NJev, NLU, NFevJac int
	// NSteps is the number of accepted steps.
	NSteps int
	// Status is 0 when the end of the span was reached, 1 when a terminal
	// event occurred and -1 when integration failed.
	Status int
	// Message describes the reason for termination.
	Message string
	// Success is Status >= 0.
	Success bool
}

Solution is the result of SolveIVP (SciPy's OdeResult).

func SolveIVP 

func SolveIVP(f ODEFunc, tSpan [2]float64, y0 []float64, opts *IVPOptions) (*Solution, error)

SolveIVP integrates the initial value problem dy/dt = f(t, y), y(t0) = y0 over tSpan = {t0, tf} (tf < t0 integrates backwards). It is SciPy's scipy.integrate.solve_ivp: 

solve_ivp(f, (t0, tf), y0)                        → SolveIVP(f, [2]float64{t0, tf}, y0, nil)
solve_ivp(f, span, y0, method="Radau", rtol=1e-6) → SolveIVP(f, span, y0, &IVPOptions{Method: Radau, RTol: 1e-6})
solve_ivp(f, span, y0, t_eval=ts, dense_output=True, events=[g])
                                                   → SolveIVP(f, span, y0, &IVPOptions{TEval: ts, DenseOutput: true, Events: []Event{{F: g}}})

where SciPy's fun(t, y) returning dydt becomes f(t, y, dydt) writing it. Invalid input returns a nil Solution and an error wrapping ErrODEInput. If the solver fails mid-way the partial Solution (Success false, Status -1) is returned together with the error.

Example

solve_ivp(lotka, (0, 15), [10, 5]) 

package main

import (
	"fmt"

	"github.com/go-numerics/integrate"
)

func main() {
	lotka := func(t float64, y, dydt []float64) {
		dydt[0] = 1.5*y[0] - y[0]*y[1]
		dydt[1] = -3*y[1] + y[0]*y[1]
	}
	sol, _ := integrate.SolveIVP(lotka, [2]float64{0, 15}, []float64{10, 5}, nil)
	fmt.Println(sol.Message, sol.NFev, len(sol.T))
}

Output:
The solver successfully reached the end of the integration interval. 470 57
Example (Events)

A terminal event: solve_ivp(ball, (0, 10), [10, 5], events=hit) with hit.terminal = True, hit.direction = -1. 

package main

import (
	"fmt"

	"github.com/go-numerics/integrate"
)

func main() {
	ball := func(t float64, y, dydt []float64) { dydt[0], dydt[1] = y[1], -9.81 }
	hit := integrate.Event{F: func(t float64, y []float64) float64 { return y[0] }, Terminal: 1, Direction: -1}
	sol, _ := integrate.SolveIVP(ball, [2]float64{0, 10}, []float64{10, 5},
		&integrate.IVPOptions{Events: []integrate.Event{hit}, RTol: 1e-8})
	fmt.Printf("%.6f %d\n", sol.TEvents[0][0], sol.Status)
}

Output:
2.025769 1
Example (Stiff)

solve_ivp(vdp, (0, 3000), [2, 0], method="Radau") with μ = 1000 (stiff) 

package main

import (
	"fmt"

	"github.com/go-numerics/integrate"
)

func main() {
	const mu = 1000
	vdp := func(t float64, y, dydt []float64) {
		dydt[0] = y[1]
		dydt[1] = mu*(1-y[0]*y[0])*y[1] - y[0]
	}
	sol, err := integrate.SolveIVP(vdp, [2]float64{0, 3000}, []float64{2, 0}, &integrate.IVPOptions{Method: integrate.Radau})
	fmt.Println(err, sol.Success, sol.NSteps < 300)
}

Output:
<nil> true true

func (*Solution) Component 

func (s *Solution) Component(i int) []float64

Component returns the history of component i, SciPy's sol.y[i].

func (*Solution) TArray 

func (s *Solution) TArray() *ndarray.Array

TArray returns the time points as a one-dimensional ndarray.

func (*Solution) YArray 

func (s *Solution) YArray() *ndarray.Array

YArray returns the states as an n×len(T) ndarray, the layout of SciPy's sol.y.

type Weight 

type Weight int

Weight selects a weight function w(x) for Quad, which then integrates f(x)·w(x) with a QUADPACK routine specialised for that weight. It is SciPy's quad(..., weight=...) argument; the SciPy wvar parameter is spread over the typed fields of QuadOptions (Omega, Alpha and Beta, C). 

const (
	// WeightNone integrates f itself (QAGS, QAGI, or QAGP with Points).
	WeightNone Weight = iota
	// WeightCos is cos(Omega·x): QAWO on a finite range, QAWF on
	// [a, +inf) or (-inf, b]. SciPy weight="cos", wvar=omega.
	WeightCos
	// WeightSin is sin(Omega·x). SciPy weight="sin", wvar=omega.
	WeightSin
	// WeightAlg is (x-a)^Alpha·(b-x)^Beta (QAWS). SciPy weight="alg",
	// wvar=(alpha, beta).
	WeightAlg
	// WeightAlgLogA is (x-a)^Alpha·(b-x)^Beta·log(x-a). SciPy "alg-loga".
	WeightAlgLogA
	// WeightAlgLogB is (x-a)^Alpha·(b-x)^Beta·log(b-x). SciPy "alg-logb".
	WeightAlgLogB
	// WeightAlgLog is (x-a)^Alpha·(b-x)^Beta·log(x-a)·log(b-x). SciPy
	// "alg-log".
	WeightAlgLog
	// WeightCauchy is 1/(x-C): the Cauchy principal value (QAWC). SciPy
	// weight="cauchy", wvar=c.
	WeightCauchy
)

Source Files

View all Source files

Jump to

Keyboard shortcuts

? : This menu
/ : Search site
f or F : Jump to
y or Y : Canonical URL
go.dev uses cookies from Google to deliver and enhance the quality of its services and to analyze traffic. Learn more.