算額（その1119）改訂版
参考文献
1). 十 岩手県胆沢町若柳市野々 個人宅 安政2年(1855)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
2). 今有如図 03046
https://w.atwiki.jp/sangaku/pages/190.html
キーワード：円3個，正方形，斜線2本
#Julia, #SymPy, #算額, #和算
はじめに
本稿は，山村の図が不適切なために，算額に記載されている「答」，「術」に一致する答えが得られなかったので，稿を改めたものである。
問は以下のとおりである。
今有方内如図設甲円二个及大小斜容乙円一个其乙円径一寸問小斜幾何
補足も加えて現代語訳すると以下のようにもなるであろう。
正方形の中に甲円を 2 個と，大小（長短）の斜線を引き，区画された領域に乙円を容れる。乙円の直径が 1 寸のとき，短い方の斜線の長さはいかほどか。
この問に添付された図は以下のようになっている。
この図を山村が彼の解釈にしたがって描いたのか，元々の算額の図を正確に写したのかは定かではない。
別の算額についての資料を検索中に，この算額を記述した文書があることを知った。。
ページの先頭部分に示した参考文献の 2 番目の「今有如図 03046」である。
そこには，山村の図とは全く異なる図が描かれていた。
2 個の甲円は正方形の左辺に接している。長短の斜線は正方形の頂点を端点とする甲円と乙円の共通接線で，短い方の斜線は下側の甲円の円周上に端点を持つ。
この図に基づいて，短い方の斜線（オレンジ色）の長さを求めてみよう。
正方形の一辺の長さを a
長い方の斜線と正方形の辺との交点を (x02, a)
短い方の斜線の甲円の円周上の端点を (x01, y01)
甲円の半径と中心座標を r1, (r1, r1), (r1, a - r1)
乙円の半径と中心座標を r2, (a - r2, y2)
とおき，以下の連立方程式を解く。
短い方の斜線の長さは sqrt((a - x01)^2 + (a - y01)^2) である。
include("julia-source.txt")
using SymPy
@syms a::positive, r1::positive, r2::positive,
y2::positive, x01::positive,
y01::positive, x02::positive
r1 = a/4
eq1 = dist2(x01, y01, a, a, r1, a - r1, r1)
eq2 = dist2(x01, y01, a, a, a - r2, y2, r2)
eq3 = dist2(x02, a, a, 0, r1, a - r1, r1)
eq4 = dist2(x02, a, a, 0, a - r2, y2, r2)
eq5 = (x01 - r1)^2 + (y01 - r1)^2 - r1^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (a, y2, x01, y01, x02))
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)
(a, y2, x01, y01, x02) = u
return [
a^2*(a^2 + 3*a*x01 - 5*a*y01 - 3*x01*y01 + 4*y01^2), # eq1
a^4 - 2*a^3*r2 - 2*a^3*x01 - 2*a^3*y2 - a^2*r2^2 + 2*a^2*r2*x01 + 2*a^2*r2*y01 + 2*a^2*r2*y2 + a^2*x01^2 + 4*a^2*x01*y2 + a^2*y2^2 + 2*a*r2^2*x01 - 2*a*r2*x01*y01 - 2*a*r2*x01*y2 - 2*a*r2*y01*y2 - 2*a*x01^2*y2 - 2*a*x01*y2^2 - r2^2*x01^2 + 2*r2*x01*y01*y2 + x01^2*y2^2, # eq2
a^2*(-a^2 + a*x02 + 4*x02^2), # eq3
-a^2*r2^2 - 2*a^2*r2*y2 + a^2*y2^2 + 2*a*r2^2*x02 + 2*a*r2*x02*y2 - 2*a*x02*y2^2 - r2^2*x02^2 + x02^2*y2^2, # eq4
-a^2/16 + (-a/4 + x01)^2 + (-a/4 + y01)^2, # eq5
]
end;
r2 = 1/2
iniv = BigFloat[10, 6.40388, 3.2, 4.9, 3.90388]
res = nls(H, ini=iniv)
([2.780776406404415, 1.7807764064044151, 0.8898484500494128, 1.3625804391381635, 1.0855823048033113], true)
乙円の直径が 1 のとき，小斜の長さは 2.36365994544375 という結果になった。
術は「乙円の直径を (√49.13 + 11.9)/8 倍すれば小斜」ということで，(√49.13 + 11.9)/8 = 2.363659945443753 となり，見事なまでに一致した。
全てのパラメータは以下のとおりである。
r2 = 0.5; a = 2.78078; x01 = 0.889848; y01 = 1.36258; x02 = 1.08558
他にも，同じような難点をはらんだ問題があるようだ。時間を無駄にした。
function draw(r2, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, y2, x01, y01, x02) = res[1]
r1 = a/4
小斜 = sqrt((a - x01)^2 + (a - y01)^2)
@printf("乙円の直径が %g のとき，小斜の長さは %.15g である。\n", 2r2, 小斜)
@printf("r2 = %g; a = %g; x01 = %g; y01 = %g; x02 = %g\n", r2, a, x01, y01, x02)
plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
circle(r1, r1, r1)
circle(r1, a - r1, r1)
circle(a - r2, y2, r2, :blue)
segment(x01, y01, a, a, :orange, lw=1)
segment(x02, a, a, 0, :magenta)
if more
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)
point(0, a, "a ", :green, :right, :vcenter)
point(a, 0, " a", :green, :center, delta=-delta/2)
point(a, a, "(a,a)", :orange, :center, :bottom, delta=delta/2)
point(x02, a, "(x02,a)", :magenta, :center, :bottom, delta=delta/2)
point(x01, y01, "(x01,y01)", :orange, :right, delta=-delta/2)
point(r1, r1, "甲円:r1,(r1,r1)", :red, :center, delta=-delta/2)
point(r1, a - r1, "甲円:r1,(r1,a-r1)", :red, :center, delta=-delta/2)
point(a - r2, y2, "乙円:r2,(a-r2,y2)", :blue, :center, delta=-delta/2)
plot!(xlims=(-6delta, a+3delta), ylims=(-4delta, a+3delta))
end
end;
draw(1/2, true)
