百一 岩手県大船渡市猪川町 田茂山町神明杜 田茂山神社（奉納年不明）
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円6個，円弧，正三角形
#Julia #SymPy #算額 #和算 #数学
円弧の中に正三角形 1 個，甲円 2 個，乙円 4 個を容れる。甲円の直径が 1 寸のとき，弦の長さを求めよ。
正三角形の一辺の長さを \(2a\)
円弧の半径と中心座標を \(R,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (x_1,\ y_0 + r_1)\)
乙円の半径と中心座標を \(r_2,\ (x_{21},\ y_0 + r_2),\ (x_{22},\ y_{22})\)
弦と \(y\) 軸の交点座標を \(y_0\)とおき，以下の連立方程式を解く。
弦の長さは \(2\sqrt{R^2 - y_0^2}\) である。
甲円の中心と正三角形の底辺の右側の頂点を結ぶ直線は円弧の中心を通り，\(x\) 軸となす角は 60° なので，弦と \(y\) 軸の交点座標は \(y_0 = R/2\) である。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positice, y0::positive, R::positive,
r1::positive, x1::positive,
r2::positive, x21::positive, x22::positive, y22::positive
x1 = r1 /√Sym(3)+ a
y0 = √Sym(3)a
R = 2y0
eq1 = x1^2 + (y0 + r1)^2 - (R - r1)^2 |> expand
eq2 = x21^2 + (y0 + r2)^2 - (R - r2)^2 |> expand
eq3 = x22^2 + y22^2 - (R - r2)^2 |> expand
eq4 = (x21 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand
eq5 = dist2(a, y0, 0, R, x22, y22, r2);
println(eq1, ", # eq1")
println(eq2, ", # eq2")
println(eq3, ", # eq3")
println(eq4, ", # eq4")
println(eq5, ", # eq5")
-8*a^2 + 20*sqrt(3)*a*r1/3 + r1^2/3, # eq1
-9*a^2 + 6*sqrt(3)*a*r2 + x21^2, # eq2
-12*a^2 + 4*sqrt(3)*a*r2 - r2^2 + x22^2 + y22^2, # eq3
a^2 + 2*sqrt(3)*a*r1/3 - 2*a*x21 + r1^2/3 - 4*r1*r2 - 2*sqrt(3)*r1*x21/3 + x21^2, # eq4
3*a^2 - 3*a*x22 - sqrt(3)*a*y22 - r2^2 + 3*x22^2/4 + sqrt(3)*x22*y22/2 + y22^2/4, # eq5
eq1 を \(a\) について解くことができる。その \(2\sqrt{3}\) 倍が弦の長さである。
ans_a = solve(eq1, a)[2]
ans_a |> println
r1*(5*sqrt(3) + 9)/12
y0 = √Sym(3)ans_a
R = 2y0
x0 = sqrt(R^2 - y0^2)
弦 = 2x0 |> factor
#= 弦 =# 弦 |> println
弦(r1 => 1/2).evalf() |> println
r1*(5*sqrt(3) + 9)/2
4.41506350946110
弦の長さは，甲円の半径の \( (5\sqrt{3} + 9)/2\)，甲円の直径の \( (5\sqrt{3} + 9)/2 = (\sqrt{75} + 9)/4 = 4.415063509461097\) 倍である。
これは「術」と一致する。
---
以下は，図を描くためにほかのパラメータを求める手順である。
まず，eq2, eq3, eq4, eq5 の \(a\) に上で求めた \(ans_a\) を代入する。
eq12 = eq2(a => ans_a) |> simplify
eq13 = eq3(a => ans_a) |> simplify
eq14 = eq4(a => ans_a) |> simplify
eq15 = eq5(a => ans_a) |> simplify
eq12 |> println
eq13 |> println
eq14 |> println
eq15 |> println
-39*r1^2/4 - 45*sqrt(3)*r1^2/8 + 15*r1*r2/2 + 9*sqrt(3)*r1*r2/2 + x21^2
-13*r1^2 - 15*sqrt(3)*r1^2/2 + 5*r1*r2 + 3*sqrt(3)*r1*r2 - r2^2 + x22^2 + y22^2
9*sqrt(3)*r1^2/8 + 9*r1^2/4 - 4*r1*r2 - 3*sqrt(3)*r1*x21/2 - 3*r1*x21/2 + x21^2
r1^2*(5*sqrt(3) + 9)^2/48 - r1*x22*(5*sqrt(3) + 9)/4 - sqrt(3)*r1*y22*(5*sqrt(3) + 9)/12 - r2^2 + 3*x22^2/4 + sqrt(3)*x22*y22/2 + y22^2/4
eq14, eq15 を解いて \(x_{21},\ x_{22}\) を求める。
(ans_x21, ans_x22) = solve([eq14, eq15], (x21, x22))[4]
#= x21 =# ans_x21 |> println
#= x22 =# ans_x22 |> println
2*sqrt(r1)*sqrt(r2) + 3*r1*(1 + sqrt(3))/4
5*sqrt(3)*r1/6 + 3*r1/2 + 2*sqrt(3)*r2/3 - sqrt(3)*y22/3
eq12 に\( ans_{x21},\ ans_{x22}\) を代入し，\(r_2\) を求める。
eq22 = eq12(x21 => ans_x21, x22 => ans_x22) |> simplify
ans_r2 = solve(eq22, r2)[1] |> simplify |> sympy.sqrtdenest |> simplify
#= r2 =# ans_r2 |> println
ans_r2(r1 => 1/2).evalf() |> println
3*r1*(4*sqrt(3) + 13)/121
0.247043841697630
eq13 に \(ans_{x21},\ ans_{x22},\ ans_{r2}\) を代入し \(y_{22}\) を求める。
eq23 = eq13(x21 => ans_x21, x22 => ans_x22) |> simplify
eq33 = eq23(r2 => ans_r2) |> simplify
ans_y22 = solve(eq33, y22)[1] |> sympy.sqrtdenest |> simplify
#= y22 =# ans_y22 |> println
ans_y22(r1 => 1/2).evalf() |> println
7*r1*(151 + 93*sqrt(3))/484
2.25678210302411
以上をまとめると，\(r_1\) が与えられたとき各パラメータの値を得ることができる。
r1 = 1/2
a = r1*(5√3 + 9)/12
r2 = 3*r1*(4√3 + 13)/121
y22 = 7*r1*(151 + 93√3)/484
x21 = 2*sqrt(r1)*sqrt(r2) + 3r1*(1 + √3)/4
x22 = 5√3*r1/6 + 3*r1/2 + 2√3*r2/3 - √3y22/3
(a, r2, x21, x22, y22)
その他のパラメータは以下の通りである。
\(r_1 = 0.5;\ R = 2.54904;\ x_0 = 2.20753;\ a = 0.735844;\ y_0 = 1.27452\)
\(x_1 = 1.02452;\ r_2 = 0.247044;\ x_{21} = 1.72743\)
\(x_{22} = 0.453996;\ y_{22} = 2.25678\)
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 1/2
a = r1*(5√3 + 9)/12
r2 = 3*r1*(4√3 + 13)/121
y22 = 7*r1*(151 + 93√3)/484
x21 = 2*sqrt(r1)*sqrt(r2) + 3r1*(1 + √3)/4
x22 = 5√3*r1/6 + 3*r1/2 + 2√3*r2/3 - √3y22/3
x1 = r1 /√3 + a
y0 = √3a
R = 2y0
x0 = sqrt(R^2 - y0^2)
θ = atand(y0, x0)
@printf("甲円の直径が %g のとき，弦の長さは %g である。\n", 2r1, 2x0)
@printf("r1 = %g; R = %g; x0 = %g; a = %g; y0 = %g; x1 = %g; r2 = %g; x21 = %g; x22 = %g; y22 = %g\n", r1, R, x0, a, y0, x1, r2, x21, x22, y22)
plot([a, 0, -a, a], y0 .+ [0, √3a, 0, 0], color=:green, lw=0.5)
circle(0, 0, R, :green, beginangle=θ, endangle=180 - θ)
circle2(x1, y0 + r1, r1, :blue)
circle2(x21, y0 + r2, r2, :red)
circle2(x22, y22, r2, :red)
segment(-x0, y0, x0, y0)
segment(0, 0, x0, y0, :gray70)
segment(0, 0, -x0, y0, :gray70)
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(x1, y0 + r1, "甲円:r1\n(x1,y0+r1)", :blue, :center, :bottom, delta=delta/2)
point(x21, y0 + r2, "乙円:r2,(x21,y0+r2) ", :black, :right, :vcenter)
point(x22, y22, " 乙円:r2,(x22,y22)", :black, :left, :vcenter)
point(0, R, " R=y0+√3a", :black, :center, :bottom, delta=delta/2)
point(0, y0, "y0", :black, :center, delta=-delta/2)
point(a, y0, "(a,y0)", :black, :center, delta=-delta/2)
point(x0, y0, "(x0,y0)", :black, :right, delta=-delta/2)
end
end;
以下のアイコンをクリックして応援してください
