Lately I have been playing with the idea of writing a simple symbolic integrator based on the book by Bronstein “Symbolic Integration 1”, and have been more or less writing it from the ground up. So I have been writing things like a polynomial division algorithm.
This was all going well until I started getting rounding issues with the function that I wrote to divide polynomials.
At this point the solution that I decided to try was to convert the values of the coefficients of all of the equations into rationals. This seems to have worked for the issue that I was having however this brings up two new issues, both of witch seem to be part of symbolics.
The first issue, and so far this seems to be just a minor issue, is that if I do
using Symbolics
@variables x
A = 1//1 * x
B = Symbolics.coeff(A, x)
println(B)
typeof(B)
I would expect that this would return a rational{Int64} for the typeof(B) like it will if you use
A = 5//1 *x
however it seems to insist on returning a Int64. This seems to happen in all situations so it is not just a consequence of how I did it here. So far this has not seemed to give me any issues however it’s not the result that I expected.
The real issue that I am having is that if I do
using Symbolics
@variables t
B = (((t^2)^3)*(((800//1)*(t^3) / ((800//1)*(t^3)
- (14//1)*t))^2)) / ((( (60//1)*(t^2))^2)
* ( (14//1)*t)^3)
B = Symbolics.simplify(B)
The result is
ERROR: MethodError: no method matching div(::PolyForm{Real}, ::Rational{Int64}, ::RoundingMode{:ToZero})
Closest candidates are:
div(::MultivariatePolynomials.AbstractPolynomialLike, ::Number, ::Any...)
@ MultivariatePolynomials ~/.julia/packages/MultivariatePolynomials/ckbfK/src/division.jl:8
div(::P, ::Real, ::RoundingMode) where P<:Dates.Period
@ Dates /usr/share/julia/stdlib/v1.9/Dates/src/periods.jl:81
div(::Integer, ::Rational, ::RoundingMode)
@ Base rational.jl:449
...
Stacktrace:
[1] div(a::PolyForm{Real}, b::Rational{Int64})
@ Base ./div.jl:47
[2] rm_gcds(ns::Vector{Any}, ds::Vector{Any})
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/polyform.jl:538
[3] simplify_div(d::SymbolicUtils.BasicSymbolic{Real})
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/polyform.jl:275
[4] (::SymbolicUtils.var"#sdiv#126")(a::SymbolicUtils.BasicSymbolic{Real})
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/polyform.jl:328
[5] call_composed
@ ./operators.jl:1034 [inlined]
[6] (::ComposedFunction{SymbolicUtils.var"#sdiv#126", typeof(quick_cancel)})(x::SymbolicUtils.BasicSymbolic{Real}; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Base ./operators.jl:1031
[7] (::ComposedFunction{SymbolicUtils.var"#sdiv#126", typeof(quick_cancel)})(x::SymbolicUtils.BasicSymbolic{Real})
@ Base ./operators.jl:1031
[8] (::SymbolicUtils.Rewriters.Walk{:post, ComposedFunction{SymbolicUtils.var"#sdiv#126", typeof(quick_cancel)}, typeof(SymbolicUtils.frac_similarterm), false})(x::SymbolicUtils.BasicSymbolic{Real})
@ SymbolicUtils.Rewriters ~/.julia/packages/SymbolicUtils/Oyu8Z/src/rewriters.jl:200
[9] simplify_fractions(x::SymbolicUtils.BasicSymbolic{Real}; polyform::Bool)
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/polyform.jl:330
[10] simplify_fractions(x::SymbolicUtils.BasicSymbolic{Real})
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/polyform.jl:322
[11] simplify(x::SymbolicUtils.BasicSymbolic{Real}; expand::Bool, polynorm::Nothing, threaded::Bool, simplify_fractions::Bool, thread_subtree_cutoff::Int64, rewriter::Nothing)
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/simplify.jl:42
[12] simplify(x::SymbolicUtils.BasicSymbolic{Real})
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/simplify.jl:16
[13] simplify(n::Num; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Symbolics ~/.julia/packages/Symbolics/BQlmn/src/Symbolics.jl:146
[14] simplify(n::Num)
@ Symbolics ~/.julia/packages/Symbolics/BQlmn/src/Symbolics.jl:146
[15] top-level scope
@ Untitled-1:17
The Julia VERSION that is running is v"1.9.3" and the Symbolics version is Symbolics v5.5.1
This is not the equation that originally gave me this error as I simplified it as much as I could without losing the actual error message. However I don’t really know what I need to do to stop this from coming up or if I need to go back and reconsider how I am doing things.