Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Try using TaylorDiff.jl instead of TaylorIntegration.jl #253

Open
nathanaelbosch opened this issue Oct 13, 2023 · 1 comment
Open

Try using TaylorDiff.jl instead of TaylorIntegration.jl #253

nathanaelbosch opened this issue Oct 13, 2023 · 1 comment
Labels
help wanted Extra attention is needed

Comments

@nathanaelbosch
Copy link
Owner

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:

"""
Compute initial derivatives of an IIP ODEProblem with TaylorIntegration.jl
"""
function taylormode_get_derivatives(u, f::SciMLBase.AbstractODEFunction{true}, p, t, q)
tT = Taylor1(typeof(t), q)
tT[0] = t
uT = similar(u, Taylor1{eltype(u)})
@inbounds @simd ivdep for i in eachindex(u)
uT[i] = Taylor1(u[i], q)
end
duT = zero(uT)
uauxT = similar(uT)
TaylorIntegration.jetcoeffs!(f, tT, uT, duT, uauxT, p)
return [evaluate.(differentiate.(uT, i)) for i in 0:q]
end

Currently I'm not sure how to best get exactly this functionality from TaylorDiff.jl. If you do, please comment!

@nathanaelbosch nathanaelbosch added the help wanted Extra attention is needed label Oct 13, 2023
@nathanaelbosch
Copy link
Owner Author

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
help wanted Extra attention is needed
Projects
None yet
Development

No branches or pull requests

1 participant