続神壁算法 〇四 中山道碓氷峠熊野社 享和元年(1801)
藤田貞資門人 上州安中驛 角田左源司親信
藤田嘉言(1807)：続神壁算法
http://www.wasan.jp/jinpeki/jinpeki.html
十一 群馬県碓氷郡松井田町峠 熊野神社 享和元年(1801)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円多数，菱形
#Julia #SymPy #算額 #和算 #数学
菱形の中に等円を数個描く（仮に11個描く）。菱形の対角線の長い方（菱長），短い方（菱平）の和と積，等円の直径が与えられたとき，等円の個数を求める術を述べよ。
菱形の中心を原点とする
菱長，菱平を \(2a,\ 2b\)
菱長と菱平の和を \(K_1\), 菱長と菱平の積を \(K_2\)
等円の半径を \(r\)
等円の中心座標を \( (x,\ y)\) として，\(x ≧ 0\) かつ \(y = 0\) の等円の個数を \(m\)，\(x = 0\) かつ \(y ≧ 0\) の等円の個数を \(n\) とする。
使う使わないはともかくとして，以下の 4 本の方程式を立てる。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, b::positive, r::positive,
m, n,
K1::positive, K2::positive
eq1 = r/(a - 2(m - 1)r) - b/sqrt(a^2 + b^2)
eq2 = r/(b - 2(n - 1)r) - a/sqrt(a^2 + b^2)
eq3 = 2a + 2b - K1
eq4 = (2a)*(2b) - K2;
1. \(r,\ m,\ n\) を与えて \(a,\ b\) を求める場合
算額問題としては最も多いパターンである。 eq1, eq2 を連立させて \(a,\ b\) を求める。SymPy の能力上，\(r,\ m,\ n\) を未知数のままとして解くことができないので，\(r,\ m,\ n\) が既知として解く。
図の場合を例として， \(r = 1/2,\ m = 4,\ n = 3\) として解く。
using SymPy
@syms a::positive, b::positive, r::positive,
m, n,
K1::positive, K2::positive
(r, m, n) = (1/2, 4, 3)
eq1 = r/(a - 2(m - 1)r) - b/sqrt(a^2 + b^2)
eq2 = r/(b - 2(n - 1)r) - a/sqrt(a^2 + b^2)
solve([eq1, eq2], (a, b))[1]
(3.90138781886600, 2.60092521257733)
2 .* (3.90138781886600, 2.60092521257733)
(7.802775637732, 5.20185042515466)
7.802775637732 + 5.20185042515466 |> println
7.802775637732 * 5.20185042515466 |> println
13.004626062886661
40.58887176852263
菱長 = 7.802775637732，菱平 = 5.20185042515466 である。
これを基礎として，図を描いたものが最初に示したものである。
後で使うことになる，菱長と菱平の和と積も 13.004626062886661, 40.58887176852263 である。
2. \(r,\ a,\ b\) を与えられて \(m,\ n\) を求める場合
eq1, eq2 は独立なので，それぞれを \(m,\ n\) について解けばよい。
using SymPy
@syms a::positive, b::positive, r::positive,
m, n,
K1::positive, K2::positive
eq1 = r/(a - 2(m - 1)r) - b/sqrt(a^2 + b^2)
eq2 = r/(b - 2(n - 1)r) - a/sqrt(a^2 + b^2);
まずは \(m\) を求める。
ans_m = solve(eq1, m)[1]
ans_m |> println
a/(2*r) + 1 - sqrt(a^2 + b^2)/(2*b)
前項のようにして別途計算しておいた結果，\(r = 1/2,\ a = 3.90138781886600,\ b = 2.60092521257733\) を与えると，\(m = 4\) が得られる。
ans_m(r => 1/2, a => 3.90138781886600, b => 2.60092521257733) |> println
4.00000000000000
次に \(n\) を求め，\(r,\ a,\ b\) を代入し実際の値を求める。
\(n = 3\) が得られる。
ans_n = solve(eq2, n)[1]
ans_n |> println
ans_n(r => 1/2, a => 3.90138781886600, b => 2.60092521257733) |> println
b/(2*r) + 1 - sqrt(a^2 + b^2)/(2*a)
3.00000000000000
\(a,\ b\) に近似値を与えると，\(m,\ n\) は整数から外れるが，四捨五入すればよい。
ans_m(r => 1/2, a => 3.9, b => 2.6) |> println
ans_n(r => 1/2, a => 3.9, b => 2.6) |> println
3.99861218113400
2.99907478742267
3. \(r,\ K_1,\ K_2\) を与えられて \(m,\ n\) を求める場合
eq1, eq2, eq3, eq4 を連立方程式として解けばよい。
using SymPy
@syms a::positive, b::positive, r::positive,
m, n,
K1::positive, K2::positive
eq1 = r/(a - 2(m - 1)r) - b/sqrt(a^2 + b^2)
eq2 = r/(b - 2(n - 1)r) - a/sqrt(a^2 + b^2)
eq3 = 2a + 2b - K1
eq4 = (2a)*(2b) - K2
res = solve([eq1, eq2, eq3, eq4], (m, n, a, b))[2] # 2 of 2
( (K1*r + K2/2 - r*sqrt(K1^2 - 4*K2) - r*sqrt(K1^2 - 2*K2))/(r*(K1 - sqrt(K1^2 - 4*K2))), (K1*r + K2/2 + r*sqrt(K1^2 - 4*K2) - r*sqrt(K1^2 - 2*K2))/(r*(K1 + sqrt(K1^2 - 4*K2))), K1/4 + sqrt(K1^2 - 4*K2)/4, K1/4 - sqrt(K1^2 - 4*K2)/4)
\(m,\ n\) はちょっと複雑な式になる。また，問では「長平の和と積」，この解答例では 長，平を \(2a,\ 2b\) としている点に注意が必要である。
長平の和と積は 13.004626062886661, 40.58887176852263 である。
\(r,\ K_1,\ K_2\) を代入して \(m,\ n\) は以下のようになる。
res[1](r => 1/2, K1 => 13.004626062886661, K2 => 40.58887176852263) |> println
res[2](r => 1/2, K1 => 13.004626062886661, K2 => 40.58887176852263) |> println
4.00000000000000
3.00000000000000
第 3，第 4 要素として \(a,\ b\) も求まる。
最初の項で求めた (3.90138781886600, 2.60092521257733) と一致する（当たり前だが）。
res[3](r => 1/2, K1 => 13.004626062886661, K2 => 40.58887176852263) |> println
res[4](r => 1/2, K1 => 13.004626062886661, K2 => 40.58887176852263) |> println
3.90138781886600
2.60092521257733
二数の和 \(u\) と積 \(v\) が与えられたとき，元の二数 \(x,\ y\) は以下の方程式を解けばわかる。
@syms u, v, x, y
eq01 = x + y ⩵ u
eq02 = x * y ⩵ v
res2 = solve( (eq01, eq02), (x, y))
2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
(u/2 - sqrt(u^2 - 4*v)/2, u/2 + sqrt(u^2 - 4*v)/2)
(u/2 + sqrt(u^2 - 4*v)/2, u/2 - sqrt(u^2 - 4*v)/2)
長平の和と積を与えて，長，平を得る。それを 1/2 すれば，\(a,\ b\) である。
事前に \(a,\ b\) を求めて，前項に従えば同じ結果になる。
res2[2][1](u => 13.004626062886661, v => 40.58887176852263) |> println
res2[2][2](u => 13.004626062886661, v => 40.58887176852263) |> println
7.80277563773200
5.20185042515466
res2[2][1](u => 13.004626062886661, v => 40.58887176852263)/2 |> println
res2[2][2](u => 13.004626062886661, v => 40.58887176852263)/2 |> println
3.90138781886600
2.60092521257733
4. 閑話休題
問題で要求されているのは，等円の個数である。
求めた \(m,\ n\) から，\(2(m - 1) + 2(n - 1) + 1 = 2m + 2n - 3\) で求めることができる。
\(m = 4,\ n = 3\) のとき，個数は 11 個である。
描画関数プログラムのソースを見る
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r, m, n) = (1/2, 4, 3)
(a, b) = (3.90138781886600, 2.60092521257733)
@printf("菱長 = %g，菱平 = %g\n", 2a, 2b)
@printf("長平和 = %g，長平積 = %g\n", 2a + 2b, 2a*2b)
plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:blue, lw=0.5)
circle(0, 0, r)
for i = 1:m - 1
circle2(i, 0, r)
end
for j = 1:n - 1
circle22(0, j, r)
end
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(a, 0, "a", :blue, :left, :bottom, delta=delta)
point(0, b, "b", :blue, :center, :bottom, delta=delta)
end
end;
draw(true)
以下のアイコンをクリックして応援してください
