Integrals with Caching Interface

Often, integral solvers allocate memory or reuse quadrature rules for solving different problems. For example, if one is going to solve the same integral for several parameters

using Integrals

prob = IntegralProblem((x, p) -> sin(x * p), (0, 1), 14.0)
alg = QuadGKJL()

solve(prob, alg)

prob = remake(prob, p = 15.0)
solve(prob, alg)

then it would be more efficient to allocate the heap used by quadgk across several calls, shown below by directly calling the library

using QuadGK
segbuf = QuadGK.alloc_segbuf()
quadgk(x -> sin(14x), 0, 1, segbuf = segbuf)
quadgk(x -> sin(15x), 0, 1, segbuf = segbuf)

Integrals.jl's caching interface automates this process to reuse resources if an algorithm supports it and if the necessary types to build the cache can be inferred from prob. To do this with Integrals.jl, you simply init a cache, solve!, replace p, and solve again. This uses the SciML init interface

using Integrals

prob = IntegralProblem((x, p) -> sin(x * p), (0, 1), 14.0)
alg = QuadGKJL()

cache = init(prob, alg)
sol1 = solve!(cache)

retcode: Success
u: 0.061661627270869046

cache.p = 15.0
sol2 = solve!(cache)

retcode: Success
u: 0.11731252752392136

The caching interface is intended for updating p and domain. Note that the types of these variables is not allowed to change. If it is necessary to change the integrand f instead of defining a new IntegralProblem, consider using FunctionWrappers.jl.

Caching for sampled integral problems

For sampled integral problems, it is possible to cache the weights and reuse them for multiple data sets.

using Integrals

x = 0.0:0.1:1.0
y = sin.(x)

prob = SampledIntegralProblem(y, x)
alg = TrapezoidalRule()

cache = init(prob, alg)
sol1 = solve!(cache)

retcode: Success
u: 0.4593145488579764

cache.y = cos.(x)   # use .= to update in-place
sol2 = solve!(cache)

retcode: Success
u: 0.8407696420884198

If the grid is modified, the weights are recomputed.

cache.x = 0.0:0.2:2.0
cache.y = sin.(cache.x)
sol3 = solve!(cache)

retcode: Success
u: 1.4114231970988789

For multi-dimensional datasets, the integration dimension can also be changed

using Integrals

x = 0.0:0.1:1.0
y = sin.(x) .* cos.(x')

prob = SampledIntegralProblem(y, x)
alg = TrapezoidalRule()

cache = init(prob, alg)
sol1 = solve!(cache)

retcode: Success
u: 11-element Vector{Float64}:
 0.0
 0.08393690598261781
 0.16703514214650947
 0.2484644183845503
 0.3274111202855099
 0.4030864385003036
 0.474734250264264
 0.5416386743258631
 0.6031312237955397
 0.6585974854457273
 0.7074832587247253

cache.dim = 1
sol2 = solve!(cache)

retcode: Success
u: 11-element Vector{Float64}:
 0.4593145488579764
 0.4570198892864838
 0.4501588380519307
 0.4387999485102908
 0.42305671493111446
 0.4030864385003037
 0.3790886556186594
 0.35130314420012265
 0.32000752789011855
 0.28551450214186125
 0.24816870986674897