#==========
Julia の修行をするときに,いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
二分法による 1 変数方程式の解
http://aoki2.si.gunma-u.ac.jp/R/bisection.html
ファイル名: bisection.jl 関数名: bisection, bisection2
翻訳するときに書いたメモ
さすがに,引数の個数・型が同じだと,同じ関数名では作れない。
==========#
function bisection(func, lower, upper; ndiv=1, epsilon = 1e-14, maxiteration = 500)
if ndiv == 1
root = bisection2(func, lower, upper)
!isnan(root) && println("answer = $root")
return root
else
x = range(lower, upper, length = ndiv)
root = [bisection2(func, x[i], x[i + 1]; epsilon, maxiteration) for i = 1:ndiv - 1]
n = length(root)
if n == 0
println("no solution")
else
root = root[.!isnan.(root)]
for i in root
println("answer = $i")
end
return root
end
end
end
function bisection2(func, lower, upper; epsilon = 1e-14, maxiteration=500)
yl = func(lower)
yu = func(upper)
yl * yu > 0 && return NaN
for i in 1:maxiteration
mid = (lower + upper) / 2
ym = func(mid)
if abs(ym) < epsilon
return mid
elseif yu * ym > 0
upper = mid
yu = ym
else
lower = mid
yl = ym
end
abs(upper - lower) < epsilon && return (lower + upper) / 2
end
NaN
end
func(x) = (x + 6.7) * (x - 3.4) * (x - sqrt(2))
bisection(func, -10.0, 10.0, ndiv=500)
bisection(x -> (x + 6.7) * (x - 3.4) * (x - sqrt(2)), -10, 10, ndiv=10)
# answer = -6.7
# answer = 1.4142135623730976
# answer = 3.4000000000000012
bisection(func, -10, 10)
# answer = -6.700000000000004
roots = bisection(x -> sin(x) / x, 0.0000001, 50, ndiv=100)
# answer = 3.1415926535897847
# answer = 6.2831853071795445
# answer = 9.424777960769305
# answer = 12.566370614359183
# answer = 15.707963267949054
# answer = 18.849555921538926
# answer = 21.991148575128346
# answer = 25.13274122871822
# answer = 28.274333882308085
# answer = 31.415926535897967
# answer = 34.557519189487856
# answer = 37.69911184307725
# answer = 40.84070449666758
# answer = 43.98229715025701
# answer = 47.123889803847334
using Plots
x = 0:0.01:50
plot(x, sin.(x) ./ x, grid=false, ylabel="\$\\frac{\\sin x}{x}\$", title="\$\\frac{\\sin x}{x}\$", label="");
x = 0:15
xticks!(π .* x, string.(x) .* "π")
vline!(roots, lw=0.5, label="");
hline!([0], lw=0.5, label="");
savefig("fig1.png")