算額(その479)
宮城県丸森町小斎日向 鹿島神社 明治9年
徳竹亜紀子,谷垣美保: 2021年度の算額調査,仙台高等専門学校名取キャンパス 研究紀要,第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf
外円内に甲,乙,丙,丁円が入っている。丁円の直径を知って,丙円の直径を求めよ。
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (0, r2), (0, r0 - r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (x4, r4)
とおいて以下の連立方程式の数値解を求める。
include("julia-source.txt")
using SymPy
@syms r1::positive, x1::positive, y1::positive, r2::positive,
r3::positive, x3::positive, r4::positive,
x4::positive, r0::positive;
r0 = 4r2
y1 = 2r2
r4 = 1//2
eq1 = x1^2 + y1^2 - (r0 - r1)^2
eq2 = x3^2 + r3^2 - (r0 - r3)^2
eq3 = x1^2 + (r0 - r2 - y1)^2 - (r1 + r2)^2
eq4 = (x3 - x1)^2 + (y1 - r3)^2 - (r1 + r3)^2
eq5 = (x1 - x4)^2 + (y1 - r4)^2 - (r1 + r4)^2
eq6 = x4^2 + (r2 - r4)^2 - (r2 + r4)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6])#, (r1, x1, r2, r3, x3, x4))
using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
v = r.zero
end
return v, r.f_converged
end;
function H(u)
(r1, x1, r2, r3, x3, x4) = u
return [
4*r2^2 + x1^2 - (-r1 + 4*r2)^2, # eq1
r3^2 + x3^2 - (4*r2 - r3)^2, # eq2
r2^2 + x1^2 - (r1 + r2)^2, # eq3
-(r1 + r3)^2 + (2*r2 - r3)^2 + (-x1 + x3)^2, # eq4
-(r1 + r4)^2 + (2*r2 - r4)^2 + (x1 - x4)^2, # eq5
x4^2 + (r2 - r4)^2 - (r2 + r4)^2, # eq6
]
end;
r0 = 4r2
y1 = 2r2
r4 = 1/2
iniv = [big"2.6", 4.3, 2.2, 1.5, 7.1, 3.0]
res = nls(H, ini=iniv)
(BigFloat[1.311396103067892771960759925894364135356352343919626632929505882095829615307969, 2.141500868793756520305245638521174546550359518581729934628195842934091980001108, 1.092830085889910643300633271578636779463626953266355527441254901746524679423308, 0.7285533905932737622004221810524245196424179688442370182941699344976831196155384, 3.569168114656260867175409397535290910917265864302883224380326404890153300001858, 1.478397839480233171313044189429409031462889497069357846136076631538706602818749], true)
r1 = 1.3114; x1 = 2.1415; r2 = 1.09283; r3 = 0.728553; x3 = 3.56917; x4 = 1.4784
丙円の直径は丁円の直径の 1.45711 倍
丙円の半径は丁円の半径の 1.4571067811865475 倍である。
この数値は術によれば sqrt(1/2)+3/4 である。
Float64(res[1][4]/r4) |> println
sqrt(1/2)+3/4 |> println
1.4571067811865475
1.4571067811865475
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, x1, r2, r3, x3, x4) = res[1]
r0 = 4r2
y1 = 2r2
r4 = 1/2
@printf("r1 = %g; x1 = %g; r2 = %g; r3 = %g; x3 = %g; x4 = %g\n", r1, x1, r2, r3, x3, x4)
@printf("丙円の直径は丁円の直径の %g 倍\n", r3/r4)
plot()
circle(0, 0, 4r2, :black)
circle4(x1, y1, r1, :magenta)
circle(0, r0 - r2, r2, :blue)
circle(0, r2 - r0, r2, :blue)
circle(0, r2, r2, :blue)
circle(0, -r2, r2, :blue)
circle4(x3, r3, r3, :green)
circle4(x4, r4, r4, :orange)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
point(x1, y1, "甲円:r1,(x1,y1)", :magenta, :center, :top, delta=-delta/2)
point(0, r0 - r2, " 乙円:r2,(0,r0-r2)", :blue, :left, :vcenter)
point(0, r2, " 乙円:r2,(0,r2)", :blue, :left, :vcenter)
point(x3, r3, "丙円::r3,(x3,r3) ", :green, :right, :vcenter)
point(x4, r4, " 丁円::r4,(x4,r4) ", :orange, :left, :top, delta=-delta/2)
point(0, r0, " r0", :black, :left, :bottom, delta=delta/3)
end
end;