Rational Type support for DifferentialEquations

General Usage
, ,
1

Hello,

I am trying to use the DifferentialEquations functions with Rational{BigInt} as data type.

for simpler cases it works as it’s shown in https://tutorials.sciml.ai/html/type_handling/01-number_types.html

Even though I turn off the adaptive option in the solve command, i get an

InexactError: BigInt(1//1000000)

and after searching around I found out that DiffEq is trying to set a tolerance (1//1000000) and convert it to BigInt. In order to reproduce the issue:

function f(u, p, t)
       return 3//1*u
end

u0 = ones(Rational{BigInt}, 10, 3)
prob = ODEProblem(f, u0, (0//1,1//1))
sol = solve(prob, RK4(), dt=1//4, adaptive=false)

This will result in the following error message:

ERROR: InexactError: BigInt(1//1000000)
Stacktrace:
 [1] Integer
   @ ./rational.jl:110 [inlined]
 [2] convert
   @ ./number.jl:7 [inlined]
 [3] __init(prob::ODEProblem{Matrix{Rational{BigInt}}, Tuple{Rational{Int64}, Rational{Int64}}, false, SciMLBase.NullParameters, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(d), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::RK4, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Rational{Int64}, dtmin::Nothing, dtmax::Rational{Int64}, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/8rtul/src/solve.jl:166
 [4] __solve(::ODEProblem{Matrix{Rational{BigInt}}, Tuple{Rational{Int64}, Rational{Int64}}, false, SciMLBase.NullParameters, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(d), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::RK4; kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:dt, :adaptive), Tuple{Rational{Int64}, Bool}}})
   @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/8rtul/src/solve.jl:5
 [5] #solve_call#26
   @ ~/.julia/packages/DiffEqBase/Lq1gG/src/solve.jl:473 [inlined]
 [6] solve_up(prob::ODEProblem{Matrix{Rational{BigInt}}, Tuple{Rational{Int64}, Rational{Int64}}, false, SciMLBase.NullParameters, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(d), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Matrix{Rational{BigInt}}, p::SciMLBase.NullParameters, args::RK4; kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:dt, :adaptive), Tuple{Rational{Int64}, Bool}}})
   @ DiffEqBase ~/.julia/packages/DiffEqBase/Lq1gG/src/solve.jl:835
 [7] #solve#31
   @ ~/.julia/packages/DiffEqBase/Lq1gG/src/solve.jl:802 [inlined]
 [8] top-level scope
   @ REPL[68]:1

is it possible to overcome this and have something working purely with Rationals, even the solution?

Thanks in advance

2

It doesn’t error for me using Julia v1.8.2 and DifferentialEquations v7.6.0

julia> function f(u, p, t)
              return 3//1*u
       end
f (generic function with 1 method)

julia> u0 = ones(Rational{BigInt}, 10, 3)
10×3 Matrix{Rational{BigInt}}:
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1

julia> prob = ODEProblem(f, u0, (0//1,1//1))
ODEProblem with uType Matrix{Rational{BigInt}} and tType Rational{Int64}. In-place: false
timespan: (0//1, 1//1)
u0: 10×3 Matrix{Rational{BigInt}}:
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1
 1//1  1//1  1//1

julia> sol = solve(prob, RK4(), dt=1//4, adaptive=false)
retcode: Success
Interpolation: 3rd order Hermite
t: 5-element Vector{Rational{Int64}}:
 0//1
 1//4
 1//2
 3//4
 1//1
u: 5-element Vector{Matrix{Rational{BigInt}}}:
 [1//1 1//1 1//1; 1//1 1//1 1//1; … ; 1//1 1//1 1//1; 1//1 1//1 1//1]
 [4331//2048 4331//2048 4331//2048; 4331//2048 4331//2048 4331//2048; … ; 4331//2048 4331//2048 4331//2048; 4331//2048 4331//2048 4331//2048]
 [18757561//4194304 18757561//4194304 18757561//4194304; 18757561//4194304 18757561//4194304 18757561//4194304; … ; 18757561//4194304 18757561//4194304 18757561//4194304; 18757561//4194304 18757561//4194304 18757561//4194304]
 [81238996691//8589934592 81238996691//8589934592 81238996691//8589934592; 81238996691//8589934592 81238996691//8589934592 81238996691//8589934592; … ; 81238996691//8589934592 81238996691//8589934592 81238996691//8589934592; 81238996691//8589934592 81238996691//8589934592 81238996691//8589934592]
 [351846094668721//17592186044416 351846094668721//17592186044416 351846094668721//17592186044416; 351846094668721//17592186044416 351846094668721//17592186044416 351846094668721//17592186044416; … ; 351846094668721//17592186044416 351846094668721//17592186044416 351846094668721//17592186044416; 351846094668721//17592186044416 351846094668721//17592186044416 351846094668721//17592186044416]
1 Like
3

What’s the motivation to simulate ODEs using rationals? Those numerators and denominators are going to get large very fast, even when you cancel common factors.

Is it just to investigate or illustrate the time-stepping scheme exactly for a few steps?

4

I have some operators that need to satisfy a relation exactly. This would imply exact energy conservation for my simulation. I tried it with floats and the relation hold up to floating point error but the simulation gets unstable. I want to show exact energy conservation even for some time steps.

Now lower orders is fine since I can write a RK4 easily, but for an 8th order operator I need something of order 9 or higher to exclude errors from the time stepper.

It seems to be a problem with my version of the package, I am trying to resolve it and will come back to you.

Thanks for the help!

5

This sounds like you might want to try using floats and a symplectic solver.

1 Like
6

Even with 0 floating point error, methods like RK4 will still exhibit lack of energy conservation. It has nothing to do with numerical error if there’s discretization error.

4 Likes