Symbolics.jl
Symbolics.jl copied to clipboard
NaN and division by 0 in solve_for
I originally posted this to the Julia #helpdesk slack on October 9th, but didn't get an answer then. I'm feeling like it's a bug.
I'm having trouble using Symbolics.solve_for
.
I'm working with vectors in three dimensional space. Because @variables
is for real numbers, I've defined
eltname(name, i) = Symbol("($name)_$i")
function symbolic_vector(name::Symbol, n::Integer)
map(variable, map(i -> eltname(name, i), 1:n))
end
This gives me a vector of n
Symbolics variables.
I think this is the approach suggested to me by (I think) Sashi at a recent CAJUN.
As a convenience, when I want to substitute numeric values for one of these vectors, I define
function vsubs(v::Vector{Symbolics.Num}, values...)
Pair.(v, [1,2,3])
end
I should be able to wrap the results of multiple calls to vsubs
in a Dict
to pass to `Symbolics.substitute'.
I can define a function that will give me the parametric formula for points on a line if I know two points on the line:
function line(p1, p2, param)
p1 + param * (p2 - p1)
end
Suppose I want to find the point where two lines intersect:
let
P1 = symbolic_vector(:P1, 3)
P2 = symbolic_vector(:P2, 3)
@variables r
l1 = line(P1, P2, r)
P3 = symbolic_vector(:P3, 3)
P4 = symbolic_vector(:P4, 3)
@variables s
l2 = line(P3, P4, s)
eq = Equation(l1, l2)
solve_for(eq, s)
end
Yes, there is no guarantee that these lines intersect, but I don't see how that is a justification to divide by zero:
Symbolics.Num[
(1//0)*(var"(P1)_1"
- var"(P3)_1"
+ (-var"(P1)_1"
+ var"(P2)_1")*r
- (-var"(P3)_1"
+ var"(P4)_1")*s),
(1//0)*(var"(P1)_2"
- var"(P3)_2"
+ (-var"(P1)_2"
+ var"(P2)_2")*r
- (-var"(P3)_2"
+ var"(P4)_2")*s),
(1//0)*(var"(P1)_3"
- var"(P3)_3"
+ (-var"(P1)_3"
+ var"(P2)_3")*r
- (-var"(P3)_3"
+ var"(P4)_3")*s)]
[1:3][1:3]
If I substitute concrete values for my points
let
P1 = symbolic_vector(:P1, 3)
P2 = symbolic_vector(:P2, 3)
@variables r
l1 = line(P1, P2, r)
P3 = symbolic_vector(:P3, 3)
P4 = symbolic_vector(:P4, 3)
@variables s
l2 = line(P3, P4, s)
subs = Dict([
vsubs(P1, 0, 0, 0)...,
vsubs(P2, 10, 10, 0)...,
vsubs(P3, 4, 0, 0)...,
vsubs(P4, 0, 4, 0)...
])
eq = Equation(l1, l2)
# Expected solution: r = 0.2, s = 0.5.
eq2 = substitute(eq, subs)
solve_for(eq2, s)
end
I can just look at eq2
and verify that r = 0.2 and s = 0.5, but the result of solve_for
still has zero denominators:
Symbolics.Num[(1//0)*(-4 + 10r + 4s), (1//0)*(10r - 4s), NaN][1:3][1:3]
Am I doing something wrong?