Solving Equations with Julia-Defined Types

One of the nice things about DifferentialEquations.jl is that it is designed with Julia's type system in mind. What this means is, if you have properly defined a Number type, you can use this number type in DifferentialEquations.jl's algorithms! [Note that this is restricted to the native algorithms of OrdinaryDiffEq.jl. The other solvers such as ODE.jl, Sundials.jl, and ODEInterface.jl are not compatible with some number systems.]

DifferentialEquations.jl determines the numbers to use in its solvers via the types that are designated by Δt and the initial condition of the problem. It will keep the time values in the same type as Δt, and the solution values in the same type as the initial condition. [Note that adaptive timestepping requires that Δt be compaible with sqrt and ^ functions. Thus Δt cannot be Integer or numbers like that if adaptive timestepping is chosen].

Let's solve the linear ODE first define an easy way to get ODEProblems for the linear ODE:

In [1]:
using DifferentialEquations
f = (t,u) -> (1.01*u)
prob_ode_linear = ODEProblem(f,1/2,(0.0,1.0));

First let's solve it using Float64s. To do so, we just need to set u₀ to a Float64 (which is done by the default) and dt should be a float as well.

In [2]:
prob = prob_ode_linear 
sol =solve(prob,RK4(),dt=1/2^(6))
println(sol)
Retcode: Success
Interpolation: 3rd order Hermite
t: 65-element Array{Float64,1}
 0.0  0.015625  0.03125  0.046875  …  0.953125  0.96875  0.984375  1.0
u: 65-element Array{Float64,1}
 0.5     
 0.507953
 0.516033
 0.524241
 0.53258 
 0.541051
 0.549658
 0.558401
 0.567283
 0.576306
 0.585473
 0.594786
 0.604247
 ⋮       
 1.15403 
 1.17239 
 1.19103 
 1.20998 
 1.22923 
 1.24878 
 1.26864 
 1.28882 
 1.30932 
 1.33015 
 1.35131 
 1.3728  

Notice that both the times and the solutions were saved as Float64. Let's change the time to use rational values:

In [4]:
prob = ODEProblem(f,1/2,(0//1,1//1));
sol = solve(prob,RK4(),dt=1//2^(6))
println(sol)
Retcode: Success
Interpolation: 3rd order Hermite
t: 65-element Array{Rational{Int64},1}
 0//1  1//64  1//32  3//64  1//16  5//64  …  61//64  31//32  63//64  1//1
u: 65-element Array{Float64,1}
 0.5     
 0.507953
 0.516033
 0.524241
 0.53258 
 0.541051
 0.549658
 0.558401
 0.567283
 0.576306
 0.585473
 0.594786
 0.604247
 ⋮       
 1.15403 
 1.17239 
 1.19103 
 1.20998 
 1.22923 
 1.24878 
 1.26864 
 1.28882 
 1.30932 
 1.33015 
 1.35131 
 1.3728  

Now let's do something fun. Let's change the solution to use Rational{BigInt} and print out the value at the end of the simulation. To do so, simply change the definition of the initial condition.

In [5]:
const rationalα = 101//100
f = (t,u) -> (rationalα*u)
prob = ODEProblem(f,BigInt(1)//BigInt(2),(0//1,1//1));
sol =solve(prob,RK4(),dt=1//2^(6))
println(sol[end])
4154032919386558883432944248380343483762044089219885824293861963690668280133800624271545564442460641100421478069957127705133139131053171319939289915624722195403241736871340745589519387833493153871994750550507166424767604170338332253959630697516305444248796250106488696552824425774652891031781638156634640665726706553562695794716367646798636566490125595141712720380867485868916531456648814528917577693417533965049279568879801863167212171389128029079788394889712773514836798543384276326561054294342851708282050876790968869065128360584151770000714515194551497614161342119347668187950856166437783338125107242946094385126468080818490755092469614835748767521966870937090173768929887202086899128132689201712566935821453568568851761907310360889009454819233203019261511646422045122043461427963067831419822632761257565485308244276118163333934078610669354885645888806741789229076806586507072844471249752898840782835318816592414922484506856439857852070928805249944302969170900303083044962139908567605824428891872081720287044135359380045755621121//3025955263570019164018502277869853398058543745963126397283707470775892712704232437030043920740033026198847216426264951289188498307633591122471111874163926157374989814610878574225506571713008520940845805558579429855707382314196875257835647882856218717417250856125102284683546912020709544155188247379716859572950811281937944702307676679453365814328593305957854274867553594143460475201489987084725797475032257007739929467758191052369579260681352907875927458926484892315482757871323905647524505025315981027903769053444125491200000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

That's one huge fraction!

Other Compatible Number Types

ArbFloats

Let's test a bunch of other number types. First I'm going to test Jeffrey Sarnoff's ArbFloats. These high precision numbers which are much faster than Bigs for less than 500-800 bits of accuracy. Having already installed Nemo and ArbFloats, I can use them in DifferentialEquations.jl via:

In [6]:
using ArbFloats
const arbalpha = ArbFloat(1.01)
f = (t,u) -> (arbalpha*u)
prob_ode_arbfloatlinear = ODEProblem(f,ArbFloat(1)/ArbFloat(2),(0//1,1//1))
sol =solve(prob_ode_arbfloatlinear,RK4(),dt=1//2^(6))
println(sol[end])
2.00341833848949843198821e-166..
WARNING: Method definition string(ArbFloats.ArbFloat{P}) in module ArbFloats at /home/crackauc/.julia/v0.6/ArbFloats/src/basics/string.jl:69 overwritten at /home/crackauc/.julia/v0.6/ArbFloats/src/basics/string.jl:220.

Bingo! ArbFloats work with DifferentialEquations.jl.

DecFP.jl

Next let's try DecFP. DecFP is a fixed-precision decimals library which is made to give both performance but known decimals of accuracy. Having already installed DecFP with Pkg.add("DecFP"), I can run the following:

In [4]:
using DecFP
const decalpha = Dec128(1.01) # Set the constant to a Dec128 as well
f = (t,u) -> (decalpha*u)
prob_ode_decfplinear = ODEProblem(f,Dec128(1)/Dec128(2),(0.0,1.0))
sol =solve(prob_ode_decfplinear,RK4(),dt=1//2^(6))
println(sol[end]); println(typeof(sol[end]))
+1372800506801159568896177661311276E-33
DecFP.Dec128

Bingo! DecFP works with DifferentialEquations.jl

Incompatible Number Systems

Decimals.jl

Install with Pkg.add("Decimals").

In [6]:
using Decimals
const decalpha2 = decimal(1.01) # Set the constant to a Dec128 as well
f = (t,u) -> (decalpha2*u)
prob_ode_decimallinear = ODEProblem(f,[decimal("1.0")]./[decimal("2.0")],(0//1,1//1))
sol =solve(prob_ode_decimallinear,RK4(),dt=1/2^(6)) #Fails
println(sol[end]); println(typeof(sol[end]))
WARNING: redefining constant decalpha2
MethodError: no method matching recursivecopy(::Array{Decimals.Decimal,1})
Closest candidates are:
  recursivecopy(::AbstractArray{T<:Number,N}) where {T<:Number, N} at /home/crackauc/.julia/v0.6/RecursiveArrayTools/src/utils.jl:2
  recursivecopy(::AbstractArray{T<:AbstractArray,N}) where {T<:AbstractArray, N} at /home/crackauc/.julia/v0.6/RecursiveArrayTools/src/utils.jl:6

Stacktrace:
 [1] #init#134(::Int64, ::Array{Rational{Int64},1}, ::Array{Rational{Int64},1}, ::Array{Rational{Int64},1}, ::Void, ::Bool, ::Void, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::Rational{Int64}, ::Void, ::Void, ::Int64, ::Rational{Int64}, ::Rational{Int64}, ::Bool, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Rational{Int64}, ::Rational{Int64}, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN, ::DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE, ::Void, ::Void, ::Bool, ::Bool, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Array{Decimals.Decimal,1},Rational{Int64},false,##11#12,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.RK4, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/crackauc/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:95
 [2] (::DiffEqBase.#kw##init)(::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Array{Decimals.Decimal,1},Rational{Int64},false,##11#12,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.RK4, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0
 [3] #solve#133(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Array{Decimals.Decimal,1},Rational{Int64},false,##11#12,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.RK4, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/crackauc/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [4] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODEProblem{Array{Decimals.Decimal,1},Rational{Int64},false,##11#12,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.RK4, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)
 [5] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/crackauc/.julia/v0.6/IJulia/src/execute_request.jl:160
 [6] eventloop(::ZMQ.Socket) at /home/crackauc/.julia/v0.6/IJulia/src/eventloop.jl:8
 [7] (::IJulia.##11#14)() at ./task.jl:335

At the time of writing this, Decimals are not compatible. This is not on DifferentialEquations.jl's end, it's on partly on Decimal's end since it is not a subtype of Number. Thus it's not recommended you use Decimals with DifferentialEquations.jl

DoubleDouble.jl

Install via Pkg.add("DoubleDouble")

In [9]:
using DoubleDouble
const doublealpha = Double(1.01) # Set the constant to a Dec128 as well
f = (t,u) -> (doublealpha*u)
prob_ode_doublelinear = ODEProblem(f,Double(1.0)/Double(2.0),(0.0,1.0))
sol =solve(prob_ode_doublelinear,RK4(),dt=1/2^(6))
println(sol[end]); println(typeof(sol[end]))
no promotion exists for DoubleDouble.Double{Float64} and Rational{Int64}

Stacktrace:
 [1] promote_type(::Type{DoubleDouble.Double{Float64}}, ::Type{Rational{Int64}}) at ./promotion.jl:161
 [2] promote at ./promotion.jl:175 [inlined]
 [3] * at ./promotion.jl:250 [inlined]
 [4] #init#134(::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Void, ::Bool, ::Void, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::Rational{Int64}, ::Void, ::Void, ::Int64, ::Rational{Int64}, ::Rational{Int64}, ::Bool, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN, ::DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE, ::Void, ::Void, ::Bool, ::Bool, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{DoubleDouble.Double{Float64},Float64,false,##17#18,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.RK4, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/crackauc/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:118
 [5] (::DiffEqBase.#kw##init)(::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{DoubleDouble.Double{Float64},Float64,false,##17#18,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.RK4, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0
 [6] #solve#133(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{DoubleDouble.Double{Float64},Float64,false,##17#18,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.RK4, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/crackauc/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [7] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODEProblem{DoubleDouble.Double{Float64},Float64,false,##17#18,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.RK4, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)
 [8] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/crackauc/.julia/v0.6/IJulia/src/execute_request.jl:160
 [9] eventloop(::ZMQ.Socket) at /home/crackauc/.julia/v0.6/IJulia/src/eventloop.jl:8
 [10] (::IJulia.##11#14)() at ./task.jl:335

DoubleDouble errors for many reason. For one, it cannot convert the default parameters form rationals to Double. But even if you set all of the default parameters to Double, it will still error because DoubleDoubles cannot multiply Ints! An issue has been filed to the DoubleDouble.jl repo for this case. If you checkout the branch from the Issue, you will see that it will still error because DoubleDouble isn't compatible with exp. There's another issue for that..

Conclusion

As you can see, DifferentialEquations.jl can use arbitrary Julia-defined number systems in its arithmetic. The only limitations are the limitations inherent in the number systems themselves. ArbFloats and ArbReals are the most feature-complete and give great performance compared to BigFloats, and thus I recommend their use when high-precision (less than ~512-800 bits) is required. DecFP is a great library for high-performance decimal numbers and works well as well. Other number systems could use some modernization.