How to find the intercept of two functions?

How do I find the intercept of two functions? I thought it would be something simple like findall(in(a), b) but this gives me Vector{Int64} with 0 elements

MWE

attack_max2 = 0.90
attack_min2 = 0.52
β_V2 = 2.0
b0V2 = 16.0 
Bi = collect(0.1:0.1:200.0)

# clearance scenarios
demand = @. 1.0/((1.0/(attack_max2 + attack_min2))*Bi^β_V2/(b0V2^β_V2+Bi^β_V2)+1/attack_max2)
boom = @. (attack_max2 - attack_min2)*Bi^β_V2/(b0V2^β_V2+Bi^β_V2)+attack_min2

findall(in(demand), boom)

I know the overlap because the plot looks like

_, i = findmin(abs.(demand-boom))
[Bi[i] demand[i] boom[i]]

With the approach in this code, it will be very improbable to find an exact match between the curves. Mattrik’s reply works but the approach is very inefficient. You need a root finding algorithm (noting that the difference between the curves is zero iff they have the same value). Bisection or Newton’s method are good starting points.

1 Like

I am open to suggestions on code as I am not hard set on it. I just did it that way because I thought it would work.

It’s basically function demand(Bi*) = function boom(Bi*) where Bi* is the intercept, but I don’t know how to do that in Julia.


attack_max2 = 0.90
attack_min2 = 0.52
β_V2 = 2.0
b0V2 = 16.0

# clearance scenarios
demand_func(x) = 1.0/((1.0/(attack_max2 + attack_min2))*x^β_V2/(b0V2^β_V2+x^β_V2)+1/attack_max2)
boom_func(x) = (attack_max2 - attack_min2)*x^β_V2/(b0V2^β_V2+x^β_V2)+attack_min2

myfunc(x) = demand_func(x) - boom_func(x)

# julia> [myfunc(0),myfunc(50)]
# 2-element Array{Float64,1}:
#   0.3799999999999999
#  -0.29324849299165867

SideStep(num) = num >= 0 ? 1 : -1

function BiSector(func::Function,x0,x1,maxdepth::Int64=1024,debugflag::Bool=false)
      fx0 = func(x0)
      fx1 = func(x1)
      if SideStep(func(x0)) * SideStep(func(x1)) == 1
          return "Error: f(x0) and f(x1) must be opposite signs"
      end
      depth = 0
      lowerx = x0
      upperx = x1
      lowervalue = fx0
      uppervalue = fx1
      while true
          midpointx = (lowerx + upperx) * 0.5
          if depth >= maxdepth
             return midpointx
          end
          lowsign = SideStep(lowervalue)
          upsign = SideStep(uppervalue)
          if lowsign*upsign == 1 || lowerx == midpointx || midpointx == upperx
              return midpointx
          end
          midpointvalue = func(midpointx)
          midsign = SideStep(midpointvalue)
          if midsign == upsign
              upperx = midpointx
              uppervalue = midpointvalue
          else
              lowerx = midpointx
              lowervalue = midpointvalue
          end
          if debugflag
              println(" depth = ",depth," lowerx = ",lowerx," upperx = ",upperx)
          end
          depth += 1
      end
  end

  flush(stdout)
  println("solution is ",BiSector(myfunc,0.0,50.0))
Starting Julia...
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.4.2 (2020-05-23)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |
solution is 14.851477511590069

1 Like

You can use Roots.jl package.

using Roots
const attack_max2 = 0.90
const attack_min2 = 0.52
const β_V2 = 2.0
const b0V2 = 16.0 

demand(x) = 1.0/((1.0/(attack_max2 + attack_min2))*x^β_V2/(b0V2^β_V2+x^β_V2)+1/attack_max2)
boom(x) = (attack_max2 - attack_min2)*x^β_V2/(b0V2^β_V2+x^β_V2)+attack_min2

delta(x) = demand(x) - boom(x)

find_zero(delta, (0.1, 200.0), Bisection())
# 14.851477511590069

Or you can use Newton method

using ForwardDiff
D(f) = x -> ForwardDiff.derivative(f, float(x))
find_zero((delta,D(delta)), 10, Roots.Newton())
# 14.851477511590062

More details can be found in package itself.

5 Likes

Hello, I thought this was the answer but looking at the value, the intercept cannot be 14.851477511590069 or 14.851477511590062 since demand(x) and boom(x) never go over 1.0.

HA, that is the x intercept… I am a silly goose.