二つの任意の曲線（関数）の共通接線を求めるプログラム

find_common_tangent.jl

# 二つの任意の曲線（関数）の、共通接線を求めるJuliaプログラム
# それほど精度が良くないので、使用する際は注意してください
# ライセンスは2-Clause BSDライセンス
using ForwardDiff
using Printf
using Roots

# 内側のfind_xerosで、根を探したときに見つからなかった時に適当に返す値
const FIND_G_ROOT_TMP = 1000.0
# 共通接線の探索範囲
const XRANGE = (-100.0, 100.0)

# 関数f(x)の接線tan_line(x)が、関数g(x)と交わるときの座標bを考え、h'(b)を返す（h(x) = f(x) - tan_line(x)）
# 共通接線を持つなら、h(x)は重根を持たなければならない→h'(b)も0になる必要がある
function find_g_root(f, g, a)
    # 関数f(x)のx = aでの接線
    tan_line = tangent_line(f, a, false)

    # 関数h(x)
    h = x -> g(x) - tan_line(x)

    # h(x)の根を探す
    res = find_zeros(h, XRANGE, atol = 0.0, rtol = 0.0)

    # 根の数をチェック
    if length(res) == 0
        # @printf "NG\n"
        # h(x)の根が見つからなかったので適当な値を返す
        return FIND_G_ROOT_TMP 
    end

    #dfx = ForwardDiff.derivative(f, a)
    #dg = x -> ForwardDiff.derivative(g, x)

    # h(x)の微分h'(x)
    dh = x -> ForwardDiff.derivative(h, x)

    #for i in 1:length(res)
    #    @printf "length(res) = %d x0 = %.7f x1(%d) = %.7f dh(a) = %.7f\n" length(res) a i res[i] dh(res[i])
    #    @printf "df(a) = %.7f dg(a) = %.7f\n" dfx dg(res[i])
    #end

    return dh(res[1])
end

# 関数f(x)の、x = aでの接線が、関数g(x)と接するときのその座標bを返す
function find_h_root(f, g, a)
    # 関数f(x)のx = aでの微分f'(a)
    dfa = ForwardDiff.derivative(f, a)

    # f'(a) - g'(b) = 0となる座標bを探す
    dh = x -> dfa - ForwardDiff.derivative(g, x) 
    res = find_zeros(dh, XRANGE, atol = 0.0, rtol = 0.0)

    # 座標bを返す
    return res[1]
end

# 関数f(x)と関数g(x)の共通接線を探す
function find_common_tangent(f, g)
    # 関数f(x)と関数g(x)が共通接線を持つときの、f(x)の方の座標を探す（複数ある）
    res = find_zeros(x -> find_g_root(f, g, x), XRANGE, atol = 0.0, rtol = 0.0)
    # @printf "length(res) = %d\n" length(res)

    return res
end

function tangent_line(f, a, printflag)
    # a における関数の値と微分係数を求める
    slope = ForwardDiff.derivative(f, a)
    y_intercept = f(a) - slope * a

    # 接線の方程式を表示
    if printflag
        @printf "y = %.7fx + (%.7f)" slope y_intercept
    end

    return x -> slope * x + y_intercept
end

# 関数f(x)と関数g(x)を定義
f(x) = x^3 - 3*x
g(x) = x^2

# 共通接線を求める
result = find_common_tangent(f, g)

# 共通接線を表示
for item in result
    @printf("共通接線:")
    tangent_line(f, item, true)
    println()
end