Removing singularities from callable functions built from symbolic expressions

Specific Domains Numerics
,
1

I have a basic question regarding building functions from symbolic expressions for which I cannot find the answer on the documentation (Symbolic Calculations and Building Callable Functions · Symbolics.jl)

Suppose that I want to build a callable function from a reasonable symbolic expression that happens to have a singularity at a particular point. In ordinary non-symbolic programming I would be able to put a special exception to remove the singularity artificially, so that NaNs did not propagate through the code. I would like to do something similar to the callable function built from a symbolic expression.

For example, let us consider the function f = x/sqrt(x^2 + y^2). There is an undefined point at x = y = 0, and the value of the function there depends on how you approach x = y = 0. If one takes y->0 first, then x->0, the result is f =1, whereas if one takes x->0 then y-> 0 the result is f= 0. I would like to force a choice of limit on this point rather than evaluating 0/0 = NaN.

Output from the Julia REPL below. I would like the final call f(0.0,0.0) = 0.0 rather than NaN. Is there an easy way to build special exceptions into the build_function call? I can wrap if statements around calls to f_callable, but this seems like a dumb solution.

using Symbolics
julia> @variables x y
2-element Vector{Num}:
x
y

julia> f = x/sqrt(x^2 + y^2)
x / sqrt(x^2 + y^2)

julia> f_callable = build_function(f, x, y, expression=Val{false})

julia> f_callable(1.0,1.0)
0.7071067811865475

julia> f_callable(0.0,1.0)
0.0

julia> f_callable(0.0,0.0)
NaN

1 Like
2

I don’t know of an automatic way to do this, but you can write your expression using ifelse.

julia> using IfElse

julia> same_type_zero(x) = zero(x)
same_type_zero (generic function with 1 method)

julia> @register_symbolic same_type_zero(x)

julia> f = IfElse.ifelse(x==0.0, same_type_zero(x), x/sqrt(x^2+y^2))
IfElse.ifelse(x == 0.0, same_type_zero(x), x / sqrt(x^2 + y^2))

julia> f_callable = build_function(f, x, y, expression=Val{false});

julia> f_callable(0.0,0.0)
0.0
3

I think you are asking for limit(f, x=>xval, y=>yval). We currently don’t have limit computation in Symbolics yet, unfortunately.