算額（その476）
宮城県丸森町小斎日向 鹿島神社 弘化5年
徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf
鈎股弦があり，鈎，股，弦を直径としその中点を中心とする 3 個の円，天円，地円，人円がある。
鈎股弦に内接する全円と，同じ直径を持つもう 1 つの全円は人円，地円，甲円に外接し，天円に内接している。
甲円は全円，地円と外接し，天円と内接している。
乙円は人円と地円の交わる部分にあり，人円と地円に内接する。
乙円の直径が与えられたとき，甲円の直径を求めよ。
鈎股弦は算額の写真を私が測定したところ，各辺が 3:4:5 の直角三角形と推測できる。
鈎, 股, 弦 の長さを 6, 8, 10 とする（偶数になるように設定する）
甲円の半径と中心座標を (r1, x1, y1)
乙円の半径と中心座標を (r2, x2, y2)
天円の半径と中心座標を (r3, x3, y3) = (弦/2, 股/2, 鈎/2)
地円の半径と中心座標を (r4, x4, |y4) = (股/2, 股/2, 0)
人円の半径と中心座標を (r5, x5, y5) = (鈎/2, 0, 鈎/2)
全円の半径と中心座標を (r6, x6, y6) = ((鈎 + 股 - 弦)/2, 鈎 + 股 - 弦)/2, 鈎 + 股 - 弦)/2)
とおいて以下の連立方程式を解く。
以下の連立方程式は条件が 7 個，未知数が 8 個なので，乙円の直径は一意に定まらない。
そこで，乙円の直径が最大になるときの解を求めるため，r1, x1, y1, x, y, r2, x2 について求めるが，それぞれは y3 を変数として含む式で表される。
include("julia-source.txt")
using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
r1::positive, x1::positive, y1::positive,
r2::positive, x2::positive, y2::positive,
r3::positive, x3::positive, y3::positive,
r4::positive, x4::positive, y4::positive,
r5::positive, x5::positive, y5::positive,
r6::positive, x6::positive, y6::positive,
x::positive, y::positive;
(鈎, 股, 弦) = (3, 4, 5).*2
(r3, x3, y3) = (弦//2, 股//2, 鈎//2) # 天円
(r4, x4, y4) = (股//2, 股//2, 0) # 地円
(r5, x5, y5) = (鈎//2, 0, 鈎//2) # 人円
(r6, x6, y6) = fill((鈎 + 股 - 弦)//2, 3) # 全円
eq1 = (x1 - x3)^2 + (y1 - y3)^2 - (r3 - r1)^2
eq2 = (x1 - x4)^2 + y1^2 - (r1 + r4)^2
eq3 = (x1 - x)^2 + (y1 - y)^2 - (r1 + r6)^2
eq4 = (x - x4)^2 + y^2 - (r6 + r4)^2
eq5 = (x - x3)^2 + (y - y3)^2 - (r3 - r6)^2
eq6 = (x4 - x2)^2 + y2^2 - (r4 - r2)^2
eq7 = x2^2 + (y5 - y2)^2 - (r5 - r2)^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, x1, y1, x, y, r2, x2))
4-element Vector{NTuple{7, Sym}}:
(3/2, 4 - sqrt(10), 9/2, 4, 6, y2/5 - 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5, 4*y2/5 - 2*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 4/5)
(3/2, 4 - sqrt(10), 9/2, 4, 6, y2/5 + 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5, 4*y2/5 + 2*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 4/5)
(3/2, sqrt(10) + 4, 9/2, 4, 6, y2/5 - 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5, 4*y2/5 - 2*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 4/5)
(3/2, sqrt(10) + 4, 9/2, 4, 6, y2/5 + 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5, 4*y2/5 + 2*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 4/5)
4 組の解が求まるが，3 番目のものが適解である。r2 はその 6 番目で，y2/5 - 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5 である。
res[3][6] |> println
y2/5 - 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5
この式は，y2 についての微分の式が 0 のとき，r2 が最大値を取ることを意味している。
pyplot(size=(400, 400), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(res[1][6], xlabel="y2", ylabel="r2", xlims=(0, 4))
r2 が最大値を取るときの y2 を求める。
solve(diff(res[3][6]))[1] |> println
9/5
y2 = 9/5 のとき，r2 が最大値を取る。そのときの実際の r2 は，以下で求められるように r2 = 1 である。
つまり，「鈎, 股, 弦 の長さを 6, 8, 10 とすると」乙円の半径は 1 である。
res[3][6](y2 => 9//5) |> println
1
甲円の半径は 3/2 である。つまり，甲円の直径は乙円の直径の 1.5 倍である。
res[3][1](y2 => 9//5) |> println
3/2
天円: r3 = 5; x3 = 4; y3 = 3
地円: r4 = 4; x4 = 4; y4 = 0
人円: r5 = 3; x5 = 0; y5 = 3
全円: r6 = 2; x6 = 2; y6 = 2 x = 4; y = 6
甲円: r1 = 1.5; x1 = 7.16228; y1 = 4.5
乙円: r2 = 1; x2 = 1.6; y2 = 1.8
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(鈎, 股, 弦) = (3, 4, 5).*2
(r3, x3, y3) = (弦/2, 股/2, 鈎/2) # 天円
(r4, x4, y4) = (股/2, 股/2, 0) # 地円
(r5, x5, y5) = (鈎/2, 0, 鈎/2) # 人円
(r6, x6, y6) = fill((鈎 + 股 - 弦)/2, 3) # 全円
y2 = 9/5
(r1, x1, y1, x, y, r2, x2) = (3/2, sqrt(10) + 4, 9/2, 4, 6, y2/5 - 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5, 4*y2/5 - 2*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 4/5)
@printf("天円: r3 = %g; x3 = %g; y3 = %g\n", r3, x3, y3)
@printf("地円: r4 = %g; x4 = %g; y4 = %g\n", r4, x4, y4)
@printf("人円: r5 = %g; x5 = %g; y5 = %g\n", r5, x5, y5)
@printf("全円: r6 = %g; x6 = %g; y6 = %g x = %g; y = %g\n", r6, x6, y6, x, y)
@printf("甲円: r1 = %g; x1 = %g; y1 = %g\n", r1, x1, y1)
@printf("乙円: r2 = %g; x2 = %g; y2 = %g\n", r2, x2, y2)
plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
circle(x3, y3, r3)
circle(x4, y4, r4, :blue)
circle(x5, y5, r5, :green)
circle(x6, y6, r6, :magenta)
circle(x1, y1, r1, :orange)
circle(x, y, r6, :magenta)
circle(x2, y2, r2, :black)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(x1, y1, " 甲円:r1,(x1,y1)", :orange, :center, :bottom, delta=delta)
#point(x2, y2, "乙円:r2,(x2,y2)", :orange, :center, :bottom, delta=delta)
point(x3, y3, " 天円:r3,(x3, y3)", :red, :left, :bottom, delta=delta)
point(x4, 0, "地円:r4,(x4, 0)", :blue, :center, :top, delta=-delta)
point(0, y5, "人円:r5,(0,y5) ", :green, :right, :bottom, delta=delta)
point(x6, y6, " 全円:r6,(x6,y6)", :magenta, :left, :bottom)
point(x, y, "全円:r6,(x,y)", :magenta, :center, :bottom, delta=delta)
point(x2, y2, " 乙円:r2,(x2,y2)", :black, :left, :top)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
end
end;
