Open nathanaelbosch opened 1 year ago
Here is a small example to get this forward, that uses the current implementation (linked above):
using ProbNumDiffEq
function lorenz(du, u, p, t)
du[1] = 10.0(u[2] - u[1])
du[2] = u[1] * (28.0 - u[3]) - u[2]
du[3] = u[1] * u[2] - (8 / 3) * u[3]
return nothing
end
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 1.0)
p = nothing
prob = ODEProblem(lorenz, u0, tspan)
order = 10
taylor_coeffs = ProbNumDiffEq.taylormode_get_derivatives(prob.u0, prob.f, prob.p, prob.tspan[1], order)
returns
11-element Vector{Vector{Float64}}:
[1.0, 0.0, 0.0]
[-10.0, 28.0, 0.0]
[380.0, -308.0, 28.0]
[-6880.0, 10920.0, -942.6666666666667]
[178000.00000000003, -201777.33333333337, 54593.77777777778]
[-3.797773333333334e6, 5.029636888888889e6, -2.256960740740741e6]
[8.827410222222221e7, -1.0087210725925928e8, 1.0874346553086421e8]
[-1.891462094814815e9, 1.8127303928395057e9, -4.805741615341564e9]
[3.7041924876543205e10, -3.093252391818924e9, 2.2031089507078738e11]
[-4.013517726836212e11, -2.5410226724751836e12, -9.768228414733674e12]
[-2.1396708997915617e13, 2.351267406708858e14, 4.275183771080162e14]
The goal of this issue would be to find a way to have exactly this functionality with TaylorDiff.jl. Ideally even with less allocations, e.g. writing directly into some pre-allocated output vector.
Not quite sure what to gain here, but at the very least the dependency could be more lightweight.
The part that would be changed is this one: https://github.com/nathanaelbosch/ProbNumDiffEq.jl/blob/70b6d57d719cdae709dc3f3e73d55570c1f8d84e/src/initialization/taylormode.jl#L35-L49
Currently I'm not sure how to best get exactly this functionality from TaylorDiff.jl. If you do, please comment!