算額（その1511）

算額（その1511）

百四 岩手県大船渡市田茂山 根城八幡宮 天保12年(1841)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

今有如図 03018
https://w.atwiki.jp/sangaku/pages/206.html

キーワード：円9個，正三角形
#Julia, #SymPy, #算額, #和算

直角三角形の中に斜線を引き，大円 1 個，中円 1 個，小円 2 個を容れる。鈎，股が 3 寸，4 寸のとき，中円の直径はいかほどか。

注：あまりきれいではないが中円，小円の中心は水平線上にある。

正三角形の一辺の長さを 2a
大円の半径と中心座標を r1,(r1, r1)
中円の半径と中心座標を r2, (0, y2), (2r2, y2)
小円の半径と中心座標を r3, (r3, y3), (3r3, y3)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms a, r1, r2, y2, r3, y3
eq1 = (2r2 - r1)^2 + (y2 - r1)^2 - (r1 + r2)^2
eq2 = (3r3 - 2r2)^2 + (y3 - y2)^2 - (r2 + r3)^2
eq3 = dist(a, 0, 0, √Sym(3)*a, r1, r1) - r1^2
eq4 = dist(a, 0, 0, √Sym(3)*a, 2r2, y2) - r2^2
eq5 = dist(a, 0, 0, √Sym(3)*a, 3r3, y3) - r3^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (a, r1, r2, y2, y3))

function H(u)
    (a, r1, r2, y2, y3) = u
    return [
        (2r2 - r1)^2 + (y2 - r1)^2 - (r1 + r2)^2,
        (3r3 - 2r2)^2 + (y3 - y2)^2 - (r2 + r3)^2,
        dist(a, 0, 0, √3*a, r1, r1) - r1^2,
        dist(a, 0, 0, √3*a, 2r2, y2) - r2^2,
        dist(a, 0, 0, √3*a, 3r3, y3) - r3^2
    ]
end;
r3 = 4/2
iniv = BigFloat[23.5, 8.6, 3.6, 20.7, 26.2]
res = nls(H, ini=iniv)

    ([23.453792062852013, 8.58468371008167, 3.6391447593398145, 20.738502725786578, 26.23085463760209], true)

大円の直径は res[1][2]*2 = 17.16936742016334 である。

術は小円の直径の 4.25 倍とのみ書いているが，不明。4.25*4 = 17

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r3 = 4/2  # 0.08933
    (a, r1, r2, y2, y3) = [1, 0.36, 0.15239, 0.88, 1.12]
    (a, r1, r2, y2, y3) = res[1]
    plot([a, 0, -a, a], [0, √3a, 0, 0], color=:magenta, lw=0.5)
    circle2(r1, r1, r1)
    circle(0, y2, r2, :blue)
    circle2(2r2, y2, r2, :blue)
    circle2(r3, y3, r3, :green)
    circle2(3r3, y3, r3, :green)
    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(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta)
        point(2r2, y2, "中円:r2,(2r2,y2)", :blue, :center, delta=-delta)
        point(3r3, y3, "小円:r3,(3r3,y3)", :black, :center, delta=-delta)
        point(a, 0, "a", :magenta, :left, :bottom, deltax=delta/2)
        point(0, √3a, "√3a", :magenta, :left, :bottom, deltax=delta/2)
    end
end;

draw(true)

算額（その1510）

算額（その1510）

九十九 岩手県江刺市 雨宝堂 現在中善観音堂 文政10年(1827)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

今有如図 03007
https://w.atwiki.jp/sangaku/pages/205.html

キーワード：円4個，直角三角形，斜線1本
#Julia, #SymPy, #算額, #和算

直角三角形の中に斜線を引き，大円 1 個，中円 1 個，小円 2 個を容れる。鈎，股が 3 寸，4 寸のとき，中円の直径はいかほどか。

鈎，股をそのまま「鈎」，「股」
大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (股 - r2, y2)
小円の半径と中心座標を r3, (x3, r3), (x32, y32)
斜線と斜辺の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms 鈎, 股, r1, x1, r2, y2, r3, x3, x32, y32, x0, y0
eq1 = dist(0, 0, 股, 鈎, x1, r1) - r1^2
eq2 = dist(0, 0, 股, 鈎, 股 - r2, y2) - r2^2
eq3 = dist(0, 0, 股, 鈎, x32, y32) - r3^2
eq4 = dist(x0, y0, 股, 0, x1, r1) - r1^2
eq6 = dist(x0, y0, 股, 0, x32, y32) - r3^2
eq7 = dist(x0, y0, 股, 0, 股 - r2, y2) - r2^2
eq8 = x0^2 + y0^2 - 股^2
eq9 = (2股 + sqrt((股 - x0)^2 + y0^2))*r1 - 股*y0
eq11 = r3/(股 - x3) - r1/(股 - x1)
eq12 = (x3 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2;

function H(u)
    (r1, x1, r2, y2, r3, x3, x32, y32, x0, y0) = u
    return [
        dist(0, 0, 股, 鈎, x1, r1) - r1^2,
        dist(0, 0, 股, 鈎, 股 - r2, y2) - r2^2,
        dist(0, 0, 股, 鈎, x32, y32) - r3^2,
        dist(x0, y0, 股, 0, x1, r1) - r1^2,
        dist(x0, y0, 股, 0, x32, y32) - r3^2,
        dist(x0, y0, 股, 0, 股 - r2, y2) - r2^2,
        x0^2 + y0^2 - 股^2,
        (2股 + sqrt((股 - x0)^2 + y0^2))*r1 - 股*y0,
        r3/(股 - x3) - r1/(股 - x1),
        (x3 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2
    ]
end;
(鈎, 股) = (3, 4)
iniv = BigFloat[0.91753, 2.72942, 0.36849, 2.26201, 0.2366, 3.6581, 3.07768, 2.00873, 3.2, 2.4]
res = nls(H, ini=iniv)

    ([0.9116963119775494, 2.7350889359326485, 0.3675444679663241, 2.2649110640673515, 0.2389194451114356, 3.668516977010949, 3.0781652486756204, 2.0099746301174206, 3.2, 2.4], true)

中円の直径は res[1][3]*2 = 0.7350889359326482 である。

function draw(r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (鈎, 股) = (3, 4)
    (r1, x1, r2, y2, r3, x3, x32, y32, x0, y0) = [0.91753, 2.72942, 0.36849, 2.26201, 0.2366, 3.6581, 3.07768, 2.00873, 3.2, 2.4]
    (r1, x1, r2, y2, r3, x3, x32, y32, x0, y0) = res[1]
    plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:magenta, lw=0.5)
    circle(x1, r1, r1)
    circle(股 - r2, y2, r2, :blue)
    circle(x3, r3, r3, :green)
    circle(x32, y32, r3, :green)
    segment(股, 0, x0, y0)
    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, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta)
        point(股 - r2, y2, "中円:r2\n(股-r2,y2)", :blue, :center, :bottom, delta=delta)
        point(x3, r3, "小円:r3,(x3,r3)", :green, :right, :vcenter, deltax=-8delta)
        point(x32, y32, "小円:r3,(x32,y32)", :green, :right, :bottom, delta=2delta, deltax=-6delta)
        point(股, 0, " 股", :magenta, :left, :bottom, delta=delta, deltax=-delta/2)
        point(股, 鈎, "(股,鈎)", :magenta, :left, :bottom, delta=delta, deltax=-delta/2)
    end
end;

draw(2, true)

算額（その1509）

算額（その1509）

六十八 岩手県一関市川崎町薄衣諏訪前 浪分神社 文政5年(1822)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

今有如図 03004
https://w.atwiki.jp/sangaku/pages/202.html

キーワード：円6個，外円，弦4本
#Julia, #SymPy, #算額, #和算

外円の中に弦を 4 本引き，区画された領域に等円 5 個を容れる。等円の直径が与えられたとき，外円の直径を求める術を述べよ。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, R - r), (x1, y1), (x2, y2)
弦の端点座標を (x01, sqrt(R^2 - x01^2)), (x02, sqrt(R^2 - x02^2)), (0, -R)
とおき，以下の連立方程式を解く。

R, r1 を記号のままにして SymPy では簡単に解くことができない。数値を代入して解けば，数値解は求まる。

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

using SymPy
@syms R::positive, r::positive, x1::positive, y1::positive,
      x2::positive, y2::negative,  x01::positive, x02::positive
eq1 = x1^2 + y1^2 - (R - r)^2
eq2 = x2^2 + y2^2 - (R - r)^2
eq3 = dist(0, -R, x01, sqrt(R^2 - x01^2), 0, R - r) - r^2
eq4 = dist(0, -R, x01, sqrt(R^2 - x01^2), x1, y1) - r^2
eq5 = dist(0, -R, x02, sqrt(R^2 - x02^2), x1, y1) - r^2
eq6 = dist(0, -R, x02, sqrt(R^2 - x02^2), x2, y2) - r^2
eq7 = (sqrt(R^2 - x02^2) + R)/x02 * y2/x2 + 1;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (R, x1, y1, x2, y2, x01, x02))  

function H(u)
    (R, x1, y1, x2, y2, x01, x02) = u
    return [
        x1^2 + y1^2 - (R - r)^2,
        x2^2 + y2^2 - (R - r)^2,
        dist(0, -R, x01, sqrt(R^2 - x01^2), 0, R - r) - r^2,
        dist(0, -R, x01, sqrt(R^2 - x01^2), x1, y1) - r^2,
        dist(0, -R, x02, sqrt(R^2 - x02^2), x1, y1) - r^2,
        dist(0, -R, x02, sqrt(R^2 - x02^2), x2, y2) - r^2,
        (sqrt(R^2 - x02^2) + R)/x02 * y2/x2 + 1
    ]
end;
r = 2
iniv = BigFloat[10, 5, 5, 7, -3, 3, 8]
res = nls(H, ini=iniv)

    ([7.509090269827044, 3.8101684797908897, 3.9790315098916964, 4.870544338179145, -2.5744656631880116, 2.2798713379711115, 6.204719439604394], true)

等円の直径が 4 のとき，外円の直径は 2*7.509090269827044 = 15.018180539654088 である。

function draw(r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, x1, y1, x2, y2, x01, x02) = res[1]
    plot()
    circle(0, 0, R, :green)
    circle(0, R - r, r)
    circle2(x1, y1, r)
    circle2(x2, y2, r)
    plot!([-x01, 0, x01], [sqrt(R^2 - x01^2), -R, sqrt(R^2 - x01^2)], color=:blue, lw=0.5)
    plot!([-x02, 0, x02], [sqrt(R^2 - x02^2), -R, sqrt(R^2 - x02^2)], color=:blue, lw=0.5)
    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, R, "R", :green, :center, :bottom, delta=delta)
        point(x1, y1, "r,(x1,y1)", :red, :center, delta=-delta)
        point(x2, y2, "r,(x2,y2)", :red, :center, delta=-delta)
        point(x01, sqrt(R^2 - x01^2), "(x01,y01)", :blue, :left, :bottom, delta=delta)
        point(x02, sqrt(R^2 - x02^2), "(x02,y02)", :blue, :left, :bottom, delta=delta)
        point(0, R - r, " R-r", :blue, :left, :vcenter)
    end
end;

draw(2, true)

術は，「外円の直径を x，等円の直径を d として 以下の d の 5 次式を解くとしている。
-d^5 + 4*d^4*x - 14*d^3*x^2 + 52*d^2*x^3 - 73*d*x^4 + 16*x^5

using SymPy
@syms x, d
A = ((16x - 73d)*x + 52d^2)*x - 14d^3
eq = (A*x + 4d^4)*x - d^5 |> expand
eq |> println

    -d^5 + 4*d^4*x - 14*d^3*x^2 + 52*d^2*x^3 - 73*d*x^4 + 16*x^5

d = 4 のときの解を求める。5 個の解のうち，3 番目のものが適解 15.0181805396541 である。

res = solve(eq(d => 4));
res[3].evalf() |> println

    15.0181805396541

 

算額（その1508）

算額（その1508）

茨城県水戸市文京 笠原神社 明治10年(1877)
https://w.atwiki.jp/sangaku/pages/201.html

キーワード：円6個，外円，弦
#Julia, #SymPy, #算額, #和算

外円の中に弦を引き，甲円 1 個，乙円 4 個を容れる。外円の直径が 12.5 寸，甲円の直径が 10.5 寸のとき，乙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, y2), (x22, R - 2r1 + r2)
とおき，以下の連立方程式を解く。

R, r1 を記号のままにして SymPy では簡単に解くことができない。数値を代入して解けば，数値解は求まる。

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

using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive, y2::negative, x22::positive
R = 12.5/2
r1 = 10.5/2
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = x22^2 + (r1 - r2)^2 - (r1 + r2)^2
eq4 = (x2 - x22)^2 + (y2 - (R - 2r1 + r2))^2 - 4r2^2;
res = solve([eq1, eq2, eq3, eq4], (r2, x2, y2, x22))[1]  # 1 of 2

    (0.750000000000000, 4.96078370824611, -2.37500000000000, 3.96862696659689)

外円の直径が 12.5 寸，甲円の直径が 10.5 寸のとき，乙円の直径は 1.5 寸である。

function draw(R, r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r2, x2, y2, x22) = (0.750000000000000, 4.96078370824611, -2.37500000000000, 3.96862696659689)
    y = R - 2r1
    x = sqrt(R^2 - y^2)
    plot()
    circle(0, 0, R, :green)
    circle(0, R - r1, r1)
    circle2(x2, y2, r2, :blue)
    circle2(x22, y + r2, r2, :blue)
    segment(-x, y, x, y)
    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, 0, "外円:R,(0,0)", :green, :center, delta=-delta)
        point(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta)
        point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, delta=-delta)
        point(x22, y + r2, "乙円:r2,(x22,y+r2)", :blue, :center, delta=-delta)
        point(0, R, "R", :green, :center, :bottom, delta=delta)
        point(0, y, " y=R-2r", :black, :center, :bottom, delta=delta)
    end
end;

draw(12.5/2, 10.5/2, true)