算額（その1413）

算額（その1413）

八四 加須市不動岡 総願寺 明治12年(1879)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円3個，正方形，斜線2本
#Julia, #SymPy, #算額, #和算

問題は判読不能ということであるが，連立方程式を立ててそれを解くというやり方であれば，何を既知数として何を未知数として解くかさえ決めれば良いので，問題が判読不能であっても問題を設定しその解を求めることは容易である。

算額の図形を文章で記述すると以下のようになる。

正方形の中に一つの頂点から対辺へ 2 本の斜線を引く。区画された領域に大円と小円を入れる。

この図形を描くためにはパラメータが 4 個必要である。正方形の右下隅を原点として描く。

• 正方形の右下の頂点座標 (a, 0)
• 斜線と左側の辺の交点座標 (0, b)
• 大円の直径 r1
• 小円の直径 r2

4 個のパラメータの相互関係は以下の 2 本の連立方程式で記述できる。

1.  大円の中心 (r1, a - r1) と 2 点 (a, 0), (0, b) を結ぶ直線の距離が r1 である。
2.  小円の中心 (r2, r2) と 2 点 (a, 0), (0, b) を結ぶ直線の距離が r2 である。

その他にも何通りかの方程式が考えられることもあるが，見かけは違っても上の式と同じ意味を持つものもある。独立な式は 2 本である。

上の連立方程式には 4 個のパラメータが含まれるが，そのうちの 2 個が既知であれば，残りの2 個のパラメータは連立方程式を解けば得られる。
既知のパラメータに直接数値を代入しても良いが，変数名そのままで方程式を解けば，解の中にも変数名が残るので，一般解が得られる。その時点で変数名にいろいろな値を与えれば，それに応じて数値解が得られる。

4 個のパラメータのうち，どの 2 個のパラメータを既知とするかは 4C2 = 6 通りある。その中には方程式を解くうえで面白いものや難しいものがある。当然「問」として適切なのはいくつかに限られる。

今回の算額図からは，以下のような「問」が考えられる。

1.  正方形の一辺と小円の直径が与えられたとき，大円の直径を求める術（すべ）を述べよ。
2.  正方形の一辺と大円の直径が与えられたとき，小円の直径を求める術（すべ）を述べよ。
3.  大円と小円の直径がが与えられたとき，正方形の一辺を求める術（すべ）を述べよ。

1.　大円の直径を求める

問1. 正方形の一辺と小円の直径が与えられたとき，大円の直径を求める術（すべ）を述べよ。

include("julia-source.txt");
# # julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive
eq1 = dist2(a, 0, 0, b, r1, a - r1, r1)
eq2 = dist2(a, 0, 0, b, r2, r2, r2);
# eq3 = a + b - sqrt(a^2 + b^2) - 2r2;

res1 = solve([eq1, eq2], (r1, b))[2]  # 1 of 2

    ((a^3 - 4*a^2*r2 + 2*a*r2^2)/(2*a^2 - 6*a*r2 + 4*r2^2), 2*r2*(-a + r2)/(-a + 2*r2))

大円の半径は a*(a^2 - 4*a*r2 + 2*r2^2)/(2*(a^2 - 3*a*r2 + 2*r2^2)) で計算できる。

res1[1] |> simplify |> println

    a*(a^2 - 4*a*r2 + 2*r2^2)/(2*(a^2 - 3*a*r2 + 2*r2^2))

正方形の一辺の長さ a = 120，小円の直径 r2 = 48 のとき，大円の直径は 70 である。
具体的な数値が挙げられている場合には，与えられる数，得られる数が整数である場合が多い。

a = 120
r2 = 48/2
a*(a^2 - 4*a*r2 + 2*r2^2)/(2*(a^2 - 3*a*r2 + 2*r2^2))

    35.0

res1[1](a => 120, r2 => 24).evalf() |> println
res1[2](a => 120, r2 => 24).evalf() |> println

    (a^3 - 4.0*a^2*r2 + 2.0*a*r2^2)/(2.0*a^2 - 6.0*a*r2 + 4.0*r2^2)
    2.0*r2*(-a + r2)/(-a + 2.0*r2)

2.　小円の直径を求める

問2. 正方形の一辺と大円の直径が与えられたとき，小円の直径を求める術（すべ）を述べよ。

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive
eq1 = dist2(a, 0, 0, b, r1, a - r1, r1)
eq2 = dist2(a, 0, 0, b, r2, r2, r2);
# eq3 = a + b - sqrt(a^2 + b^2) - 2r2;

res2 = solve([eq1, eq2], (r2, b))[1]  # 1 of 4

    ((-r1*sqrt(2*a^2 - 4*a*r1 + r1^2) + (a - 2*r1)*(a - sqrt((a^2*(a - 2*r1)^2 + (r1*sqrt(2*a^2 - 4*a*r1 + r1^2) - (a - r1)^2)^2)/(a - 2*r1)^2)) + (a - r1)^2)/(2*(a - 2*r1)), (-r1*sqrt(2*a^2 - 4*a*r1 + r1^2) + (a - r1)^2)/(a - 2*r1))

小円の半径はかなり複雑な式になる。

a = 120
r1 = 70/2
(-r1*sqrt(2*a^2 - 4*a*r1 + r1^2) + (a - 2*r1)*(a - sqrt((a^2*(a - 2*r1)^2 + (r1*sqrt(2*a^2 - 4*a*r1 + r1^2) - (a - r1)^2)^2)/(a - 2*r1)^2)) + (a - r1)^2)/(2*(a - 2*r1))

    24.0

res2[1](a => 120, r1 => 70/2).evalf() |> println
res2[2](a => 120, r1 => 70/2).evalf() |> println

    0.5*(-2.0*r1*(0.5*a^2 - a*r1 + 0.25*r1^2)^0.5 + (a - 2.0*r1)*(a - 0.5*((4.0*a^2*(0.5*a - r1)^2 + 4.0*(r1*(0.5*a^2 - a*r1 + 0.25*r1^2)^0.5 - 0.5*(a - r1)^2)^2)/(0.5*a - r1)^2)^0.5) + (a - r1)^2)/(a - 2.0*r1)
    (-2.0*r1*(0.5*a^2 - a*r1 + 0.25*r1^2)^0.5 + (a - r1)^2)/(a - 2.0*r1)

3.　正方形の一辺の長さを求める

問3. 大円と小円の直径がが与えられたとき，正方形の一辺を求める術（すべ）を述べよ。

この場合は，SymPy の能力的に連立方程式を解くのは難しい。逐次的に解いていく。

iusing SymPy
@syms a::positive, b::positive, r1::positive, r2::positive
eq1 = dist2(a, 0, 0, b, r1, a - r1, r1) |> x -> x/a
eq2 = dist2(a, 0, 0, b, r2, r2, r2) |> x -> x/(a*b);

まず，eq2 を解いて b を求める。

ans_b = solve(eq2, b)[1];
ans_b |> println

    2*r2*(a - r2)/(a - 2*r2)

eq1 の b に ans_b を代入し，eq11 とする。eq11 を解いて a を求める。

eq11 = eq1(b => ans_b) |> simplify |> numerator;

ans_a = solve(eq11, a);

5 組の解が得られるが 5 番目のものが形式上は虚数解であるが，虚数部は実質的に 0 なので適解である。a = 120

ans_a[5](r1 => 35, r2 => 24).evalf() |> println

    120.0 - 8.470329472543e-22*I

ans_b に a, r1, r2 を代入すれば b が得られる。

ans_b(r1 => 35, r2 => 24, a => 120) |> println

    64

function draw(type, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho") 
    if type == 1
        (a, r2) = (120, 24)
        (r1, b) = ((a^3 - 4*a^2*r2 + 2*a*r2^2)/(2*a^2 - 6*a*r2 + 4*r2^2), 2*r2*(-a + r2)/(-a + 2*r2))
    elseif type == 2
        (a, r1) = (120, 35)
        (r2, b) = ((-r1*sqrt(2*a^2 - 4*a*r1 + r1^2) + (a - 2*r1)*(a - sqrt((a^2*(a - 2*r1)^2 + (r1*sqrt(2*a^2 - 4*a*r1 + r1^2) - (a - r1)^2)^2)/(a - 2*r1)^2)) + (a - r1)^2)/(2*(a - 2*r1)), (-r1*sqrt(2*a^2 - 4*a*r1 + r1^2) + (a - r1)^2)/(a - 2*r1))
    elseif type == 3
        (r1, r2) = (35, 24)
        a = 120
        b = 2*r2*(a - r2)/(a - 2*r2)
    end
    @printf("a = %g;  b = %g;  r1 = %g;  r2 = %g\n", a, b, r1, r2)
    plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
    circle(r1, a - r1, r1)
    circle(r2, r2, r2, :blue)
    circle(a - r2, a - r2, r2, :blue)
    segment(0, b, a, 0, :magenta)
    segment(a - b, 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)
        point(r1, a - r1, "大円:(r1,a-r1)", :red, :center, delta=-delta)
        point(r2, r2, "小円:(r2,r2)", :blue, :center, delta=-delta)
        point(0, b, "b ", :magenta, :right, :vcenter)
        point(a, 0, "a", :green, :center, delta=-delta/2)
        plot!(xlims=(-4delta, a + 5delta), ylims=(-4delta, a + 5delta))
    end
end;

draw(3, true)