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]]
hwjsnc
June 4, 2020, 11:11am
3
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.
