算額（その1190）
(12) 京都市伏見区御香宮門前町 御香宮神社（ごこうのみやじんじゃ）文久3年(1863)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円6個，直線上
直線上に，天円と地円が互いに接して載っている。天円，地円と直線の隙間に，午，未，...子の 11 円が互いに接し，かつ，天円，地円のいずれかとも接してて入っている。
子円，亥円の直径がそれぞれ 64 寸，225 寸のとき，午円の直径はいかほどか。
解が得られるまでの過程が長いので，結論を先に述べておく。
午円の直径は 1427.16049382716 寸である。
---
まず，デカルトの円定理を用いて，それぞれの円の半径の数式解を求める。
include("julia-source.txt");
using SymPy
@syms 天::positive, 地::positive, 午::positive
#(天, 地) = (1854.1323234573401, 4951.3384889946465)
k午 = 1/天 + 1/地 + 2sqrt(1/(天*地))
午 = 1/k午
k未 = 1/地 + 1/午 + 2sqrt(1/(地*午))
未 = 1/k未
k酉 = 1/地 + 1/未 + 2sqrt(1/(地*未))
酉 = 1/k酉
k戌 = 1/地 + 1/酉 + 2sqrt(1/(地*酉))
戌 = 1/k戌
k亥 = 1/地 + 1/戌 + 2sqrt(1/(地*戌))
亥 = 1/k亥;
天, 地, 午, 未, 酉, 戌, 亥
略1
k巳 = 1/天 + 1/午 + 2sqrt(1/(天*午))
巳 = 1/k巳
k辰 = 1/天 + 1/巳 + 2sqrt(1/(天*巳))
辰 = 1/k辰
k卯 = 1/天 + 1/辰 + 2sqrt(1/(天*辰))
卯 = 1/k卯
k寅 = 1/天 + 1/卯 + 2sqrt(1/(天*卯))
寅 = 1/k寅
k丑 = 1/天 + 1/寅 + 2sqrt(1/(天*寅))
丑 = 1/k丑
k子 = 1/天 + 1/丑 + 2sqrt(1/(天*丑))
子 = 1/k子;
巳,辰,卯,寅,丑,子
略2
天円と地円の大きさにより，末円（亥円と子円）の大きさが異なる。
そこで，問にあるような，子と亥の直径が 64 寸，225 寸となるのは，天と地がいくつのときか，数値解を求める。
@syms 子2, 亥2
eq1 = 子 - 64//2
eq2 = 亥 - 225//2;
# solve([eq1, eq2], (天, 地))
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 Float64.(v), r.f_converged
end;
function H(u)
(天, 地) = u
return [
略3
]
end;
iniv = BigFloat[13000, 7500]./2
res = nls(H, ini=iniv)
([1854.1323234573401, 4951.3384889946465], true)
天円の半径は r1 = 1854.1323234573401，地円の半径は r2 = 4951.3384889946465 であることがわかった。
最初に戻り，デカルトの円定理で各円の半径を数値として求めることはできるが，ここでは図を書いて解が正しいことを確認するために，2 円の間の外接関係（ピタゴラスの定理）に基づく連立方程式を解き，半径と位置座標（円の中心の x 座標）を求める。
以下，午，未，酉，戌，亥，および巳，辰，卯，寅，丑，子の半径と中心座標を逐次求める。
天円の半径と中心座標を r1, (x1, r1)
地円の半径と中心座標を r2, (0, r2)
午円の半径と中心座標を r3, (x3, r3)
:
亥円の半径と中心座標を r7, (x7, r7)
巳円の半径と中心座標を r14, (x14, r14)
:
子円の半径と中心座標を r19, (x19, r19)
とする。
---
午の半径と中心の x 座標
@syms r1::positive, r2::positive, r3::positive, x1::positive, x3::positive
(r1, r2) = (1854.1323234573401, 4951.3384889946465)
eq1 = x1^2 + (r2 - r1)^2 - (r2 + r1)^2
eq2 = x3^2 + (r2- r3)^2 - (r2 + r3)^2
eq3 = (x1- x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
res3 = solve([eq1, eq2, eq3], (r3, x1, x3))[1] # 1 of 2
(713.580246913579, 6059.84710593375, 3759.34959349594)
未の半径と中心の x 座標
@syms r4::positive, x4::positive
(r3, x1, x3) = res3
eq4 = x4^2 + (r2 - r4)^2 - (r2 + r4)^2
eq5 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
res4 = solve([eq4, eq5], (r4, x4))[1] # 1 of 2
(374.902031440026, 2724.89769192995)
酉の半径と中心の x 座標
@syms r5::positive, x5::positive
(r4, x4) = res4
eq6 = x5^2 + (r2 - r5)^2 - (r2 + r5)^2
eq7 = (x4 - x5)^2 + (r4 - r5)^2 - (r4 + r5)^2
res5 = solve([eq6, eq7], (r5, x5))[1] # 1 of 2
(230.559556786704, 2136.89345314506)
戌の半径と中心の x 座標
@syms r6::positive, x6::positive
(r5, x5) = res5
eq8 = x6^2 + (r2 - r6)^2 - (r2 + r6)^2
eq9 = (x5 - x6)^2 + (r5 - r6)^2 - (r5 + r6)^2
res6 = solve([eq8, eq9], (r6, x6))[1] # 1 of 2
(155.979085849216, 1757.61799176434)
亥の半径と中心の x 座標
問の通り，直径が 225 寸になった。
@syms r7::positive, x7::positive
(r6, x6) = res6
eq10 = x7^2 + (r2 - r7)^2 - (r2 + r7)^2
eq11 = (x6 - x7)^2 + (r6 - r7)^2 - (r6 + r7)^2
res7 = solve([eq10, eq11], (r7, x7))[1] # 1 of 2
(112.500000000000, 1492.68292682927)
---
巳の半径と中心の x 座標
@syms r14::positive, x14::positive
eq24 = (x1 - x14)^2 + (r1 - r14)^2 - (r1 + r14)^2
eq25 = (x3 - x14)^2 + (r3 - r14)^2 - (r3 + r14)^2
res24 = solve([eq24, eq25], (r14, x14))[1] # 1 of 2
(271.777959183670, 4640.11149825784)
辰の半径と中心の x 座標
@syms r15::positive, x15::positive
(r14, x14) = res24
eq26 = (x1 - x15)^2 + (r1 - r15)^2 - (r1 + r15)^2
eq27 = (x14 - x15)^2 + (r14 - r15)^2 - (r14 + r15)^2
res25 = solve([eq26, eq27], (r15, x15))[1] # 1 of 2
(142.121439792364, 5033.17879459785)
卯の半径と中心の x 座標
@syms r16::positive, x16::positive
(r15, x15) = res25
eq28 = (x1 - x16)^2 + (r1 - r16)^2 - (r1 + r16)^2
eq29 = (x15 - x16)^2 + (r15 - r16)^2 - (r15 + r16)^2
res26 = solve([eq28, eq29], (r16, x16))[1] # 1 of 2
(87.1712696766909, 5255.78972294575)
寅の半径と中心の x 座標
@syms r17::positive, x17::positive
(r16, x16) = res26
eq29 = (x1 - x17)^2 + (r1 - r17)^2 - (r1 + r17)^2
eq30 = (x16 - x17)^2 + (r16 - r17)^2 - (r16 + r17)^2
res27 = solve([eq29, eq30], (r17, x17))[1] # 1 of 2
(58.8727931190581, 5399.06590555266)
丑の半径と中心の x 座標
@syms r18::positive, x18::positive
(r17, x17) = res27
eq31 = (x1 - x18)^2 + (r1 - r18)^2 - (r1 + r18)^2
eq32 = (x17 - x18)^2 + (r17 - r18)^2 - (r17 + r18)^2
res28 = solve([eq31, eq32], (r18, x18))[1] # 1 of 2
(42.4114263002606, 5499.00346858998)
子の半径と中心の x 座標
問の通り，子の直径が 64 寸になった。
@syms r19::positive, x19::positive
(r18, x18) = res28
eq33 = (x1 - x19)^2 + (r1 - r19)^2 - (r1 + r19)^2
eq34 = (x18 - x19)^2 + (r18 - r19)^2 - (r18 + r19)^2
res29 = solve([eq33, eq34], (r19, x19))[1] # 1 of 2
(31.9999999999989, 5572.68292682926)
各円の直径は以下のとおりである。
天 = 3708.26; 地 = 9902.68; 午 = 1427.16; 未 = 749.804; 酉 = 461.119; 戌 = 311.958; 亥 = 225
巳 = 543.556; 辰 = 284.243; 卯 = 174.343; 寅 = 117.746; 丑 = 84.8229; 子 = 64
午円の直径は 1427.16049382716 寸である。
「術」では 1600 寸としている。そもそも，算額の図では天円のほうが地円より大きいという違いがあるのだが。
function draw(lens, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2) = (1854.1323234573401, 4951.3384889946465)
(r3, x1, x3) = (713.580246913579, 6059.84710593375, 3759.34959349594)
(r4, x4) = (374.902031440026, 2724.89769192994)
(r5, x5) = (230.559556786703, 2136.89345314506)
(r6, x6) = (155.979085849216, 1757.61799176433)
(r7, x7) = (112.499999999999, 1492.68292682927)
(r14, x14) = (271.777959183669, 4640.11149825784)
(r15, x15) = (142.121439792360, 5033.17879459787)
(r16, x16) = (87.1712696766871, 5255.78972294577)
(r17, x17) = (58.8727931190587, 5399.06590555266)
(r18, x18) = (42.4114263002605, 5499.00346858998)
(r19, x19) = (31.9999999999989, 5572.68292682926)
@printf("天 = %g; 地 = %g; 午 = %g; 未 = %g; 酉 = %g; 戌 = %g; 亥 = %g\n", 2r1, 2r2, 2r3, 2r4, 2r5, 2r6, 2r7)
@printf("巳 = %g; 辰 = %g; 卯 = %g; 寅 = %g; 丑 = %g; 子 = %g\n", 2r14, 2r15, 2r16, 2r17, 2r18, 2r19)
plot()
circlef(0, r2, r2)
circlef(x1, r1, r1, :blue)
circlef(x3, r3, r3, 1)
circlef(x4, r4, r4, 2)
circlef(x5, r5, r5, 3)
circlef(x6, r6, r6, 4)
circlef(x7, r7, r7, 5)
circlef(x14, r14, r14, 6)
circlef(x15, r15, r15, 7)
circlef(x16, r16, r16, 8)
circlef(x17, r17, r17, 9)
circlef(x18, r18, r18, 10)
circlef(x19, r19, r19, 11)
if more && lens
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
annotate!(x3, r3, text("午", :white, 16))
annotate!(x4, r4, text("未", :white, 16))
annotate!(x5, r5, text("酉", :white, 16))
annotate!(x6, r6, text("戌", :white, 16))
annotate!(x7, r7, text("亥", :white, 16))
annotate!(x14, r14, text("巳", :white, 16))
annotate!(x15, r15, text("辰", :white, 16))
annotate!(x16, r16, text("卯", :white, 16))
annotate!(x17, r17 + delta, text("寅", :white, 16))
annotate!(x18 + delta/4, r18 + delta, text("丑", :white, 16))
annotate!(x19 + delta/2, r19 + delta, text("子", :white, 16))
end
lens && plot!(xlims=(1000, 6000), ylims=(0, 1500), size=(1000, 1000))
end;
