# Way to find positive roots by treating complex numbers as negatives? \[SymEngine + Roots\]

**URL:** <https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688>\
**Category:** Numerics\
**Tags:** question, roots\
**Created:** [June 15, 2018, 7:25am UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688 "2018-06-15T07:25:41Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![djsegal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djsegal/32/13752_2.png) [@djsegal](https://discourse.julialang.org/u/djsegal)\
**Post date:** [June 15, 2018, 7:25am UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/1 "2018-06-15T07:25:41Z")

</div>

Let’s say you had some hideous nonlinear equation like so:

```julia
cur_equation = -1.0 + 32.5*(0.237334965250543*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) + 8.02921230974917e-05*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258) + 2.38498439733623e-10*((I_P^0.96)^9.6774193548387*(0.907736861673663*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)*(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)) + (0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 0.82398621004115*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)^2)*(1.2 + 994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 9.67240598811777e-05*(I_P^0.96)^3.2258064516129))^0.5)/(I_P*(0.771372182862621*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) - 0.000133194034060396*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258)))

```

* * *

However the value of I\_P stands for a real, positive quantity – i.e. a current.

Because of the funkiness in the equation, when you try to get roots from:

- `Roots.find_zeros(cur_equation, 0.1, 100.0)`

You get **complex solutions** and the root solver breaks (see the following error):

```julia
MethodError: no method matching isless(::Complex{Float64}, ::Int64)
Closest candidates are:
  isless(::Missings.Missing, ::Any) at /Users/dan/.julia/v0.6/Missings/src/Missings.jl:74
  isless(::AbstractFloat, ::Real) at operators.jl:98
  isless(::ForwardDiff.Dual{Tx,V,N} where N where V<:Real, ::Integer) where Tx at /Users/dan/.julia/v0.6/ForwardDiff/src/dual.jl:104

```

* * *

The way I’ve worked around this is doing some `lambdify`-foo a la:

```julia
cur_lambda = lambdify(cur_equation)

cur_func = function (work_I_P)
  cur_value = cur_lambda(complex(float(work_I_P)))
  iszero(imag(cur_value)) || return float(-100.0)
  return Real(cur_value)
end

cur_I_P_list = Roots.find_zeros(cur_equation, 0.1, 100.0)

```

And this has rewarded me with:

```julia
cur_I_P_list == [22.2564, 27.8827, 29.8268]

```

The problem is that this is really slow.

* * *

Has anyone had experience with a similar problem?

How could this be sped up?

// I don’t have a second per root solve

---

<div class="post-metadata">

**Author:** ![djsegal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djsegal/32/13752_2.png) [@djsegal](https://discourse.julialang.org/u/djsegal)\
**Post date:** [June 15, 2018, 7:43am UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/2 "2018-06-15T07:43:34Z")

</div>

For reference, here is an image of the function with the complex values replaced with the value `-5`

 ![51%20AM](https://global.discourse-cdn.com/julialang/original/3X/9/0/90186863b2f45977dc0d78f9b1a4d08f196fedb0.png)

The complex values happen at low values of current. (and are representative of impossible to achieve plasmas)

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 15, 2018, 8:28am UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/3 "2018-06-15T08:28:50Z")

</div>

I would consider introducing some

\zeta = (I\_P^{0.96})^{3.2258064516129}

then the equation seems to be in integer powers (particularly -1, 1, 2, there may be others I missed) of \zeta.

**EDIT:** fixed mistake pointed out by @mzaffalon.

---

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [June 15, 2018, 11:22am UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/4 "2018-06-15T11:22:20Z")

</div>

Not all I\_P have the numerical factor 0.96.

---

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [June 15, 2018, 11:24am UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/5 "2018-06-15T11:24:21Z")

</div>

Except for the numerical factor, @Tamas_Papp’ suggestion is a sound one in that a change of variables may help. Indeed there are recurrent factors (I\_P^{3.09677419354838})^{(0.3226I3.09677419354838P+6363.67829013947)}.

And do you really need to keep all the terms?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 15, 2018, 11:25am UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/6 "2018-06-15T11:25:28Z")

</div>

I now see that it was actually an exponent, corrected it, thanks for pointing it out.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [June 15, 2018, 12:52pm UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/7 "2018-06-15T12:52:30Z")

</div>

Here is a session using `IntervalRootFinding.jl` on the original function, although I agree that it would be better to simplify the expression.

```julia
f(I_P) = -1.0 + 32.5*(0.237334965250543*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) + 8.02921230974917e-05*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258) + 2.38498439733623e-10*((I_P^0.96)^9.6774193548387*(0.907736861673663*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)*(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)) + (0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 0.82398621004115*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)^2)*(1.2 + 994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 9.67240598811777e-05*(I_P^0.96)^3.2258064516129))^0.5)/(I_P*(0.771372182862621*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) - 0.000133194034060396*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258)))

using IntervalArithmetic, IntervalRootFinding

rts = roots(f, -100..101, Bisection, 1e-1)

julia> rts
19-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([29.8768, 29.9757], :unknown)
 Root([29.7779, 29.8769], :unknown)
 Root([29.6805, 29.778], :unknown)
 Root([27.9235, 28.0195], :unknown)
 Root([27.8277, 27.9236], :unknown)
 Root([27.7334, 27.8278], :unknown)
 Root([22.633, 22.732], :unknown)
 Root([22.4337, 22.5327], :unknown)
 Root([22.3348, 22.4338], :unknown)
 Root([22.2374, 22.3349], :unknown)
 Root([22.1868, 22.2375], :unknown)
 Root([22.137, 22.1869], :unknown)
 Root([22.0381, 22.1371], :unknown)
 Root([21.9392, 22.0382], :unknown)
 Root([21.8419, 21.9393], :unknown)
 Root([21.743, 21.842], :unknown)
 Root([19.2784, 19.3759], :unknown)
 Root([19.1811, 19.2785], :unknown)
 Root([19.0852, 19.1812], :unknown)

julia> rts = roots(f, rts, Bisection, 1e-2)
20-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([29.833, 29.8393], :unknown)
 Root([29.8269, 29.8331], :unknown)
 Root([29.8207, 29.827], :unknown)
 Root([29.8145, 29.8208], :unknown)
 Root([27.8871, 27.8932], :unknown)
 Root([27.8811, 27.8872], :unknown)
 Root([27.8752, 27.8812], :unknown)
 Root([27.8692, 27.8753], :unknown)
 Root([22.2796, 22.2858], :unknown)
 Root([22.2735, 22.2797], :unknown)
 Root([22.2674, 22.2736], :unknown)
 Root([22.2614, 22.2675], :unknown)
 Root([22.2553, 22.2615], :unknown)
 Root([22.2493, 22.2554], :unknown)
 Root([22.2433, 22.2494], :unknown)
 Root([22.2374, 22.2434], :unknown)
 Root([22.231, 22.2375], :unknown)
 Root([22.2246, 22.2311], :unknown)
 Root([19.1811, 19.1871], :unknown)
 Root([19.1749, 19.1812], :unknown)

julia> roots(f, rts, Newton)
4-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([29.8267, 29.8269], :unique)
 Root([27.8826, 27.8828], :unique)
 Root([22.2562, 22.2566], :unique)
 Root([19.1825, 19.1833], :unknown)

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [June 15, 2018, 12:55pm UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/8 "2018-06-15T12:55:16Z")

</div>

```julia
julia> roots(f, rts, Newton, 1e-10)
3-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([29.8267, 29.8268], :unique)
 Root([27.8826, 27.8827], :unique)
 Root([22.2564, 22.2565], :unique)

```

So I agree with the `Roots` package that these are the only 3 real roots of the equation in the interval -100…100.

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [June 15, 2018, 2:03pm UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/9 "2018-06-15T14:03:07Z")

</div>

If you replace `log` with `safer_log(x) = x <= 0 ? NaN : log(x)` and make the whole thing a function of `I_P` and not a symbolic expression it runs through `find_zeros` just fine and is reasonably fast. However, converting with `lambdify` makes this trick harder. If you had to, you could get the expression to lambdify with `convert(Expr, cur_equation)`, splice in such a substitute, then call `find_zeros`. Using `MacroTools`, something like this seems to work, which might help in general, but isn’t very friendly or speedy:

```julia
u = convert(Expr, cur_equation)
ex = MacroTools.postwalk(x -> @capture(x, log(xs__)) ? :(safer_log($(xs...))) : x, u) # slowish
fn = eval(Expr(:function, Expr(:call, gensym(), :(I_P)), ex))
find_zeros(fn, 1.0, 100.0)

```

---

<div class="post-metadata">

**Author:** ![djsegal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djsegal/32/13752_2.png) [@djsegal](https://discourse.julialang.org/u/djsegal)\
**Post date:** [June 15, 2018, 11:19pm UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/10 "2018-06-15T23:19:14Z")

</div>

> [@dpsanders](#):
>
> Here is a session using `IntervalRootFinding.jl` on the original function

After trying to do,

```julia
using IntervalArithmetic, IntervalRootFinding

function test_sanders()
    f(I_P) = -1.0 + 32.5*(0.237334965250543*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) + 8.02921230974917e-05*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258) + 2.38498439733623e-10*((I_P^0.96)^9.6774193548387*(0.907736861673663*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)*(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)) + (0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 0.82398621004115*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)^2)*(1.2 + 994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 9.67240598811777e-05*(I_P^0.96)^3.2258064516129))^0.5)/(I_P*(0.771372182862621*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) - 0.000133194034060396*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258)))

    @time rts = roots(f, -100..101, Bisection, 1e-1)
    @time rts = roots(f, rts, Bisection, 1e-2)
    @time roots(f, rts, Newton)
end

test_sanders()

```

I get,

```julia
  6.549414 seconds (21.20 M allocations: 797.128 MiB, 6.75% gc time)
  5.895622 seconds (20.01 M allocations: 752.655 MiB, 6.61% gc time)
DomainError:
log will only return a complex result if called with a complex argument. Try log(complex(x)).

Stacktrace:
 [1] nan_dom_err at ./math.jl:300 [inlined]
 [2] log at ./math.jl:419 [inlined]
 [3] (::#f#1)(::Float64) at ./In[1]:4
 [4] guarded_mid(::#f#1, ::IntervalArithmetic.Interval{Float64}) at /Users/dan/.julia/v0.6/IntervalRootFinding/src/IntervalRootFinding.jl:48
 [5] N at /Users/dan/.julia/v0.6/IntervalRootFinding/src/newton.jl:37 [inlined]

```

1. Is this the fastest I can make it to catch all three roots?
2. Why does your code sample not trigger the DomainError for `log` calls?

---

<div class="post-metadata">

**Author:** ![djsegal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djsegal/32/13752_2.png) [@djsegal](https://discourse.julialang.org/u/djsegal)\
**Post date:** [June 15, 2018, 11:30pm UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/11 "2018-06-15T23:30:58Z")

</div>

> [@j\_verzani](#):
>
> If you had to, you could get the expression to lambdify with `convert(Expr, cur_equation)` , splice in such a substitute, then call `find_zeros`

Why would you do this if it was slower?

* * *

Just rehashing this quick, this is simple setup:

```julia
using SymEngine, Roots, MacroTools

I_P = SymEngine.symbols("I_P")
cur_equation = -1.0 + 32.5*(0.237334965250543*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) + 8.02921230974917e-05*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258) + 2.38498439733623e-10*((I_P^0.96)^9.6774193548387*(0.907736861673663*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)*(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)) + (0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 0.82398621004115*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129)^2)*(1.2 + 994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 9.67240598811777e-05*(I_P^0.96)^3.2258064516129))^0.5)/(I_P*(0.771372182862621*(0.0155717178747169 + 1.42646372678281e-11*(I_P^0.96)^6.4516129032258) - 0.000133194034060396*(-(0.210126719741548 + 1.0*(-1.41012671974155 - (994.155975738568*(I_P^0.96)^(-3.2258064516129)/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + 18861.5885198702*(I_P^0.96)^(-3.2258064516129)*(21212.2609671316*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)/(1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129)) + log((1 + 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))))/(1 - 5303.06524178289*(1.2 + 6.08327420636338e-05*(I_P^0.96)^3.2258064516129)*(I_P^0.96)^(-3.2258064516129))) + 0.00012773744412246*(I_P^0.96)^3.2258064516129))^2 + 3.5558715118258e-08*(I_P^0.96)^6.4516129032258))) 

```

This is your suggestion:

```julia
safer_log(x) = x <= 0 ? NaN : log(x)
u = convert(Expr, cur_equation)

@time ex = MacroTools.postwalk(x -> @capture(x, log(xs__)) ? :(safer_log($(xs...))) : x, u) # slowish
@time fn = eval(Expr(:function, Expr(:call, gensym(), :(I_P)), ex)) 
@time Roots.find_zeros(fn, 0.01, 100.0)

-------
  2.028164 seconds (102.16 k allocations: 4.982 MiB)
  0.010968 seconds (1.15 k allocations: 86.956 KiB)
  0.153697 seconds (172.94 k allocations: 7.850 MiB)

```

And this is the existing code from the initial question,

```julia
@time cur_lambda = lambdify(cur_equation)
@time cur_func = function (work_I_P)
    cur_value = cur_lambda(complex(float(work_I_P)))
    iszero(imag(cur_value)) || return float(-100.0)
    return Real(cur_value)
end
@time Roots.find_zeros(cur_func, 0.01, 100.0)

-------
  0.020340 seconds (22.95 k allocations: 831.798 KiB)
  0.000018 seconds (19 allocations: 1.483 KiB)
  1.060586 seconds (1.24 M allocations: 40.453 MiB, 1.03% gc time)

```

---

<div class="post-metadata">

**Author:** ![djsegal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djsegal/32/13752_2.png) [@djsegal](https://discourse.julialang.org/u/djsegal)\
**Post date:** [June 15, 2018, 11:39pm UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/12 "2018-06-15T23:39:30Z")

</div>

> [@Tamas\_Papp](#):
>
> I would consider introducing some
> 
> ζ=(I0.96P)3.2258064516129 \zeta = (I\_P^{0.96})^{3.2258064516129}

The reason for this factor is because we substitute expressions for radius and magnetic field into the current equation (through the `G_blank` composite functions).

![14%20PM](https://global.discourse-cdn.com/julialang/original/3X/0/8/0888902fa340b28c8c779e597c8a15d876c9d13c.png)

Where R and B equations come in pairs for various constraints:

![20%20PM](https://global.discourse-cdn.com/julialang/original/3X/5/a/5af17b3b97b1b8455d84152b69e458fdaaecf5dd.png)

* * *

Therefore simplification would have to be done either manually for every constraint or digitally – which can take a few second.

Also, if you made your substitution, there is still a factor of I\_P floating around that will have a fractional exponent.

// i.e. the one that is currently sitting to a power of one:

```julia
... /(I_P*(0.771372182862621* ...

```

So I guess part of the question: is having one fractional exponent any better than having a bunch?

---

<div class="post-metadata">

**Author:** ![djsegal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djsegal/32/13752_2.png) [@djsegal](https://discourse.julialang.org/u/djsegal)\
**Post date:** [June 17, 2018, 6:19am UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/13 "2018-06-17T06:19:18Z")

</div>

Follow-up, is there a way to handle exponentiation with this trick to?

```julia
DomainError:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

```

And the next one after that would be negative exponents &nbsp; &nbsp; // to complete the set 😛

(basically can we throw away type stability for symbolic root solving functions)

* * *

**edit:** is `lambdify` what’s slow? (b/c of the `invokelatest` [call](https://github.com/symengine/SymEngine.jl/blob/9c393c39a80f972fd876b833b8fdee918df3980e/src/subs.jl#L128))

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [June 17, 2018, 1:11pm UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/14 "2018-06-17T13:11:19Z")

</div>

There is overhead with `lambdify`, about twice the time for smaller cases, but in my example the slowness comes from walking the expression tree. You can speed it up by converting to a string, replacing the calls as desired, then converting back to an expr. As for complex numbers, many of the algorithms in `find_zero` can use complex numbers, but none of those that use bracketing (as `find_zeros` does).

---

<div class="post-metadata">

**Author:** ![djsegal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/djsegal/32/13752_2.png) [@djsegal](https://discourse.julialang.org/u/djsegal)\
**Post date:** [June 19, 2018, 12:42pm UTC](https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/16 "2018-06-19T12:42:40Z")

</div>

For the interested, a robust `find_zeros` for complex numbers was proposed in:

[https://github.com/JuliaMath/Roots.jl/pull/108](https://github.com/JuliaMath/Roots.jl/pull/108)

// it’s able to find all the roots I need in around `1 second` !!!