integrate

// 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.

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})

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)

package main

import (
	"fmt"

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

func init() { register("lotka", lotka) }

// The README's SolveIVP example: Lotka–Volterra with RK45, the default
// method, and the default tolerances. The markers are the solver's own
// accepted steps, sol.T: no interpolation.
func lotka() *plot.Figure {
	f := 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(f, [2]float64{0, 15}, []float64{10, 5}, nil)
	if err != nil {
		panic(err)
	}
	prey := make([]float64, len(sol.T))
	predators := make([]float64, len(sol.T))
	for k, y := range sol.Y {
		prey[k], predators[k] = y[0], y[1]
	}
	fig, axs := plot.Subplots(1, 1, plot.FigSize(6.4, 4))
	ax := axs[0][0]
	ax.Plot(sol.T, prey, plot.Marker("o"), plot.MarkerSize(3), plot.Label("prey"))
	ax.Plot(sol.T, predators, plot.Marker("s"), plot.MarkerSize(3), plot.Label("predators"))
	ax.SetTitle(fmt.Sprintf("SolveIVP, RK45: %d steps, %d evaluations", sol.NSteps, sol.NFev))
	ax.SetXLabel(`time $t$`)
	ax.Legend()
	return fig
}

package main

import (
	"fmt"

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

func init() { register("vanderpol", vanderpol) }

// The README's stiff example: van der Pol with mu = 1000 over [0, 3000],
// solved with Radau. An explicit method needs millions of steps here.
func vanderpol() *plot.Figure {
	const mu = 1000
	f := 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(f, [2]float64{0, 3000}, []float64{2, 0},
		&integrate.IVPOptions{Method: integrate.Radau})
	if err != nil {
		panic(err)
	}
	y0 := make([]float64, len(sol.T))
	for k, y := range sol.Y {
		y0[k] = y[0]
	}
	fig, axs := plot.Subplots(1, 1, plot.FigSize(6.4, 3.6))
	ax := axs[0][0]
	ax.Plot(sol.T, y0, plot.Marker("o"), plot.MarkerSize(2))
	ax.SetTitle(fmt.Sprintf(`Radau, $\mu = 1000$: %d steps, %d LU decompositions`, sol.NSteps, sol.NLU))
	ax.SetXLabel(`time $t$`)
	ax.SetYLabel(`$y_0$`)
	return fig
}

package main

import (
	"fmt"
	"math"

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

func init() { register("quad", quad) }

// The README's Quad example: an integrand with a singularity at 0, which
// QUADPACK's extrapolation handles. The title prints what Quad returned.
func quad() *plot.Figure {
	f := func(x float64) float64 { return math.Log(x) / math.Sqrt(x) }
	r, err := integrate.Quad(f, 0, 1, nil)
	if err != nil {
		panic(err)
	}
	x := linspace(0.002, 1, 300)
	y := make([]float64, len(x))
	for i, v := range x {
		y[i] = f(v)
	}
	fig, axs := plot.Subplots(1, 1, plot.FigSize(6.4, 4))
	ax := axs[0][0]
	ax.Plot(x, y, plot.Label(`$\log x / \sqrt{x}$`))
	ax.FillBetween(x, y, nil, plot.Alpha(0.3))
	ax.SetYLim(-12, 0.5)
	ax.SetTitle(fmt.Sprintf("Quad: %.12f, error estimate %.1e, %d evaluations", r.Value, r.AbsErr, r.NEval))
	ax.SetXLabel(`$x$`)
	ax.Legend(plot.Loc("lower right"))
	return fig
}