Fixed-Point Numbers

[gekaremi]

Good day to all!

I began to study the problem of a spinning pendulum (or, more precisely, a rotating asteroid). Most articles and papers use a floating point type variable to store the rotation angle. From the outside, this does not seem to be the best solution, since all angles are equivalent and limited to a range of 2pi. So I decided to try using fixed-point numbers for this task.

To begin, I tried to create a minimally simple example:

using DifferentialEquations, FixedPointNumbers

f(u,p,t) = u
u0 = Fixed{Int32,7}(0.5)
tspan = (Fixed{Int32,7}(0.0),Fixed{Int32,7}(1.0))
prob = ODEProblem(f,u0,tspan)
sol = solve(prob) 

which result to:

ERROR: LoadError: InexactError: trunc(Int32, 1280000000000)
Stacktrace:
 [1] throw_inexacterror(::Symbol, ::Any, ::Int64) at ./boot.jl:567
 [2] checked_trunc_sint at ./boot.jl:589 [inlined]
 [3] toInt32 at ./boot.jl:626 [inlined]
 [4] Type at ./boot.jl:716 [inlined]
 [5] convert at ./number.jl:7 [inlined]
 [6] round at ./int.jl:549 [inlined]
 [7] Type at /home/q/.julia/packages/FixedPointNumbers/3QAEg/src/fixed.jl:13 [inlined]
 [8] Type at /home/q/.julia/packages/FixedPointNumbers/3QAEg/src/fixed.jl:15 [inlined]
 [9] (::getfield(DiffEqBase, Symbol("#kw##__init")))(::NamedTuple{(:default_set,),Tuple{Bool}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Fixed{Int32,7},Tuple{Fixed{Int32,7},Fixed{Int32,7}},false,Nothing,ODEFunction{false,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Vern9,Rodas5{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}},AutoSwitch{Vern9,Rodas5{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./none:0
 [10] #__solve#257(::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:default_set,),Tuple{Bool}}}, ::Function, ::ODEProblem{Fixed{Int32,7},Tuple{Fixed{Int32,7},Fixed{Int32,7}},false,Nothing,ODEFunction{false,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Vern9,Rodas5{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}},AutoSwitch{Vern9,Rodas5{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/q/.julia/packages/OrdinaryDiffEq/miOSH/src/solve.jl:6
 [11] #__solve at ./none:0 [inlined] (repeats 5 times)
 [12] #__solve#2(::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::ODEProblem{Fixed{Int32,7},Tuple{Fixed{Int32,7},Fixed{Int32,7}},false,Nothing,ODEFunction{false,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Nothing) at /home/q/.julia/packages/DifferentialEquations/UZIcy/src/default_solve.jl:15
 [13] #__solve at ./none:0 [inlined]
 [14] #__solve#1 at /home/q/.julia/packages/DifferentialEquations/UZIcy/src/default_solve.jl:5 [inlined]
 [15] __solve at /home/q/.julia/packages/DifferentialEquations/UZIcy/src/default_solve.jl:2 [inlined]
 [16] #solve#429 at /home/q/.julia/packages/DiffEqBase/ZQVwI/src/solve.jl:41 [inlined]
 [17] solve(::ODEProblem{Fixed{Int32,7},Tuple{Fixed{Int32,7},Fixed{Int32,7}},false,Nothing,ODEFunction{false,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}) at /home/q/.julia/packages/DiffEqBase/ZQVwI/src/solve.jl:27
 [18] top-level scope at none:0
 [19] include at ./boot.jl:317 [inlined]
 [20] include_relative(::Module, ::String) at ./loading.jl:1041
 [21] include(::Module, ::String) at ./sysimg.jl:29
 [22] exec_options(::Base.JLOptions) at ./client.jl:229
 [23] _start() at ./client.jl:421
in expression starting at /home/q/FG/bagreport/fixednumbers.jl:7

and with other library:

using DifferentialEquations, FixedPointDecimals

f(u,p,t) = u
u0 = FixedDecimal{Int,4}(0.5)
tspan = (FixedDecimal{Int,4}(0.0),FixedDecimal{Int,4}(1.0))
prob = ODEProblem(f,u0,tspan)
sol = solve(prob) 

result to:

ERROR: LoadError: InexactError: Int64(Int64, 1//1000000)
Stacktrace:
 [1] Type at ./rational.jl:85 [inlined]
 [2] convert(::Type{FixedDecimal{Int64,4}}, ::Rational{Int64}) at /home/q/.julia/packages/FixedPointDecimals/Z4Eui/src/FixedPointDecimals.jl:279
 [3] Type at /home/q/.julia/packages/FixedPointDecimals/Z4Eui/src/FixedPointDecimals.jl:109 [inlined]
 [4] (::getfield(DiffEqBase, Symbol("#kw##__init")))(::NamedTuple{(:default_set,),Tuple{Bool}}, ::typeof(DiffEqBase.__init), ::ODEProblem{FixedDecimal{Int64,4},Tuple{FixedDecimal{Int64,4},FixedDecimal{Int64,4}},false,Nothing,ODEFunction{false,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Vern9,Rodas5{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}},AutoSwitch{Vern9,Rodas5{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./none:0
 [5] #__solve#257(::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:default_set,),Tuple{Bool}}}, ::Function, ::ODEProblem{FixedDecimal{Int64,4},Tuple{FixedDecimal{Int64,4},FixedDecimal{Int64,4}},false,Nothing,ODEFunction{false,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Vern9,Rodas5{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}},AutoSwitch{Vern9,Rodas5{0,false,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType},Rational{Int64},Float64}}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/q/.julia/packages/OrdinaryDiffEq/miOSH/src/solve.jl:6
 [6] #__solve at ./none:0 [inlined] (repeats 5 times)
 [7] #__solve#2(::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::ODEProblem{FixedDecimal{Int64,4},Tuple{FixedDecimal{Int64,4},FixedDecimal{Int64,4}},false,Nothing,ODEFunction{false,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Nothing) at /home/q/.julia/packages/DifferentialEquations/UZIcy/src/default_solve.jl:15
 [8] #__solve at ./none:0 [inlined]
 [9] #__solve#1 at /home/q/.julia/packages/DifferentialEquations/UZIcy/src/default_solve.jl:5 [inlined]
 [10] __solve at /home/q/.julia/packages/DifferentialEquations/UZIcy/src/default_solve.jl:2 [inlined]
 [11] #solve#429 at /home/q/.julia/packages/DiffEqBase/ZQVwI/src/solve.jl:41 [inlined]
 [12] solve(::ODEProblem{FixedDecimal{Int64,4},Tuple{FixedDecimal{Int64,4},FixedDecimal{Int64,4}},false,Nothing,ODEFunction{false,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}) at /home/q/.julia/packages/DiffEqBase/ZQVwI/src/solve.jl:27
 [13] top-level scope at none:0
 [14] include at ./boot.jl:317 [inlined]
 [15] include_relative(::Module, ::String) at ./loading.jl:1041
 [16] include(::Module, ::String) at ./sysimg.jl:29
 [17] exec_options(::Base.JLOptions) at ./client.jl:229
 [18] _start() at ./client.jl:421
in expression starting at /home/q/FG/bagreport/fixeddecim.jl:7

System info:

packages:

julia> pkgs = Pkg.installed();
julia> pkgs


Dict{String,Union{Nothing, VersionNumber}} with 8 entries:
  "DynamicalSystems"      => v"1.2.0"
  "DifferentialEquations" => v"6.3.0"
  "IJulia"                => v"1.17.0"
  "OrdinaryDiffEq"        => v"5.3.0"
  "FixedPointDecimals"    => v"0.3.0"
  "Plots"                 => v"0.23.1"
  "PyPlot"                => v"2.8.0"
  "FixedPointNumbers"     => v"0.5.3"

Julia:

Version 1.0.1
Ubuntu ⛬  julia/1.0.1-0ubuntu2

(k)ubuntu:

lsb_release -a
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 18.10
Release:        18.10
Codename:       cosmic

It is interesting that in other sources and literature, fixed-point numbers also seem to be quite rarely used for differential equations. This can be understood, but for some problems, like that, fixed-point numbers have some advantages.

If you are interested in a physical problem, then there are many articles on this topic, I can offer this as I studied it: The rotation states predominant among the planetary satellites https://arxiv.org/abs/1312.5236v1

Thanks!

[ChrisRackauckas]

New number types always find fun things. This is as far as I could get without hitting an impasse:

using OrdinaryDiffEq, FixedPointNumbers

f(u,p,t) = u
u0 = Fixed{Int32,7}(0.5)
tspan = (Fixed{Int32,7}(0.0),Fixed{Int32,7}(1.0))
prob = ODEProblem(f,u0,tspan)
sol = solve(prob,Tsit5(),dtmin = eltype(prob.tspan)(1//10^(5)),
                         abstol = Fixed{Int32,7}(0.1),
                         reltol = Fixed{Int32,7}(0.1))

Here's the difficulties of the number type that were uncovered. First, eltype(prob.tspan)(1//10^(10)) is how the dtmin parameter defaults which gives an error.

InexactError: trunc(Int32, 1280000000000)
throw_inexacterror(::Symbol, ::Any, ::Int64) at boot.jl:583
checked_trunc_sint at boot.jl:605 [inlined]
toInt32 at boot.jl:642 [inlined]
Type at boot.jl:732 [inlined]
convert at number.jl:7 [inlined]
round at int.jl:549 [inlined]
Type at fixed.jl:13 [inlined]
Fixed{Int32,7}(::Rational{Int64}) at fixed.jl:15
top-level scope at none:0

Making it larger is a good start. Then something funky was happening and I found out that abstol and reltol were defaulting to zero, since convert(Fixed{Int32,7},oneunit(Fixed{Int32,7})*1//10^6) == 0 since the defaults are lower than eps(Fixed{Int32,7}). Fantastic. So I just raised the defaults for now to see if it would work. Then I got an error which has the simple MWE:

Fixed{Int32,7}(0.1)^Fixed{Int32,7}(0.1)

^ not defined for Fixed{Int32,7}
error(::String, ::String, ::Type) at error.jl:42
no_op_err(::String, ::Type) at promotion.jl:388
^(::Fixed{Int32,7}, ::Fixed{Int32,7}) at promotion.jl:393
top-level scope at none:0

So these numbers won't work with adaptive time stepping until they know how to exponentiate. But fixed time steps have

sol = solve(prob,Tsit5(),dtmin = eltype(prob.tspan)(1//10^(5)),
                         adaptive = false, dt = Fixed{Int32,7}(0.1))

which has this error and MWE:

ceil(Int64,Fixed{Int32,7}(0.1))

MethodError: no method matching ceil(::Type{Int64}, ::Fixed{Int32,7})
Closest candidates are:
  ceil(::Type{T<:Integer}, !Matched::Integer) where T<:Integer at int.jl:552
  ceil(::Type{T<:Integer}, !Matched::Float16) where T<:Integer at float.jl:360
  ceil(::Type{T<:Union{Signed, Unsigned}}, !Matched::BigFloat) where T<:Union{Signed, Unsigned} at mpfr.jl:303
  ...
top-level scope at none:0

So it won't work in DiffEq, but it's not a me problem it's a you problem :P. Hopefully those MWEs make it pretty clear what methods need to be added to the FixedPointNumbers.jl.

[gekaremi]

I agree that this is my problem :b Many thanks for these examples, they helped to figure it out!