No, I did not.
AFAIU, the Hermite interpolation in SciML here interpolates along the independent variable of the ODE (in my case, time t).
I was talking about a Chebyshev expansion in terms of a parameter of the problem (called q).
Anyway, I’ve now finally come up with a MWE. The direct method of ODE integration driven by numerical quadrature proposed by @stevengj would look like this:
using OrdinaryDiffEq
using Integrals
using BenchmarkTools
# Goal: Integrate a function with a some weight function with a peak.
wt(q) = sin(π*q)^151 + 0.1
# The function f to be integrated is an ODE solution, evaluated at a set of time
# points collected from some graph upstream.
const tpoints = rand(Float64, 24) * 5.
# (In real life, the peak occurs in f, not in wt.)
# right-hand side of the ODE
rhs(u) = - u^2 + 0.8 * u
rhs_alternate(u) = exp(-0.1 - cos(u)^2)
# Approach using direct ODE solution driven by integration
function direct(rhs)
# ODE initial condition depends on q
prob1 = ODEProblem((u, _, _) -> rhs(u), 1., (0., 5.))
prob(q) = remake(prob1, u0=q) # trying to save on problem construction
sol(q) = solve(prob(q), Tsit5())
f(q) = sum(sol(q).(tpoints))
# now evaluate the integral
iprob = IntegralProblem((q, _) -> wt(q)*f(q), (0.,1.))
function ival(;reltol=nothing)
solve(iprob, QuadGKJL(); reltol)
end
return ival
end
# and benchmark
a = direct(rhs);
a(;reltol=1e-9); a(;reltol=1e-3);
@bprofile a(;reltol=1e-9) # this gives around 7ms and 8.5MB allocated on my computer
@bprofile a(;reltol=1e-3) # this gives around 300μs and 500KB
aa = direct(rhs_alternate);
@bprofile aa(;reltol=1e-3) # this gives around 600μs and 500KB
From this I’ve learned that the time depends strongly on the desired precision and less strongly on the rhs of the ODE. So it can indeed be quite fast for modest precision.
My current implementation is less nice and compact:
(code continued from above):
using ApproxFun
using ApproxFunBase: ClenshawPlan, chebyshev_clenshaw
# Alternative approach based on expansion
function expanded(rhs)
qspace = Chebyshev(0.0..1.0)
qfun = Fun(qspace)
# k will be the fixed interpolation order
qpts(k) = collect(points(qspace, k))
# machinery for expansion and evaluation
can_qpts(k) = tocanonical.(qspace, qpts(k))
plan(k) = ClenshawPlan(Float64, qspace, k, k)
# ODE in coefficient space
Dop(g::Fun) = rhs(g)
function rhs!(dcoeff, coeff::AbstractVector, _, _)
ufun = Fun(qspace, coeff)
dcoeff .= @view coefficients(Dop(ufun))[1:length(coeff)]
end
# initial condition
function ginit(k::Integer)
cid = coefficients(qfun)
g0 = zeros(k)
g0[1:length(cid)] .= cid
g0
end
# a coupled ODE system of size k
aprob(k) = ODEProblem(rhs!, ginit(k), (0., 5.))
asol(k) = solve(aprob(k), Tsit5())
# collect terms over the time point list. This is where combination of
# expansions comes into play. We now do it in value space
function afpts(k)
s = asol(k) # one-time ODE solution for the coefficients
cp = can_qpts(k)
pl = plan(k)
workspace = zeros(k)
for t in tpoints
workspace .+= chebyshev_clenshaw(s(t), cp, pl)
end
return workspace
# sum(evaluate.((s(t),), qspace, p) for t in tpoints)
end
# expand also the weight, truncated to k, evaluated on the points
wtpts(k) = let
w = Fun(wt, qspace)
c = coefficients(w)[1:k]
return evaluate.((c,), qspace, qpts(k))
end
# and solve the integral
function aival(k)
vals = afpts(k) .* wtpts(k)
return sum(Fun(qspace, transform(qspace, vals)))
end
return aival
end
b = expanded(); b(20); b(50);
@bprofile b(20) # around 1ms and 770KB allocated
@bprofile b(50) # around 5ms and 3.2MB
ba = expanded(rhs_alternate); ba(20);
@bprofile ba(20) # around 3s and 266MB allocated (!)
I found it interesting that it can be faster than the direct version for a simple polynomial form of the ODE, but it really crawls for an exponential form. My guess is that the polynomial rhs is a narrowly banded operator in coefficient space? In my real problem, I in fact have such a simple rhs.
Also, I don’t quite know how to make a fair comparison in terms of precision achieved, as I don’t know the relation between expansion order k and reltol.
In my previous attempts, the direct approach was a lot slower than the expanded one; this may have been due to just requesting the default high precision, and having a simple ODE that favors expansion.
Any comments are very welcome! The first version is certainly much nicer to read. Is there something more I can do to speed it up?