裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

算額(その1531)

2025年01月10日 | Julia

算額(その1531)

四十五 岩手県一関市真滝 熊野白山滝神社 明治13年(1880)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03061
https://w.atwiki.jp/sangaku/pages/327.html

キーワード:円3個,長方形,斜線2本
#Julia, #SymPy, #算額, #和算

長方形の中に大円と 2 本の斜線を設け,小円 2 個を容れる。小斜が 1 寸のとき,大斜はいかほどか。

山村の図は不正確であるが,算額の図はこうあるべし(大円は 4 つの接点を持つ)と推測すれば,解くのに支障はない。

長方形の長辺,短辺を a, b
斜線と長方形の長辺との交点座標を (c, b)
大円の半径と中心座標を r1, (a - r1, r1); b = 2r1
小円の半径と中心座標を r2, (x2, r2), (r2, b - r2)
とおき,以下の連立方程式を解く。

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

using SymPy
@syms a::positive, b::positive, c::positive,
      r1::positive, r2::positive, x2::positive, K::positive
b = 2r1
eq1 = b + c - sqrt(b^2 + c^2) - 2r2
eq2 = (a - r1 - x2)^2 + (r1 -  r2)^2 - (r1 + r2)^2
eq3 = dist(0, 0, c, b, x2, r2) - r2^2
eq4 = K - sqrt(b^2 + c^2)
eq5 = dist(0, 0, c, b, a - r1, r1) - r1^2;

function H(u)
    (a, c, r1, r2, x2) = u
    b = 2r1
    return [
        b + c - sqrt(b^2 + c^2) - 2r2,
        (a - r1 - x2)^2 + (r1 -  r2)^2 - (r1 + r2)^2,
        dist(0, 0, c, b, x2, r2) - r2^2,
        K - sqrt(b^2 + c^2),
        dist(0, 0, c, b, a - r1, r1) - r1^2
    ]
end;
K = 1
iniv = BigFloat[1.1, 0.3, 0.5, 0.1, 0.2]
res = nls(H, ini=iniv)

    ([1.12, 0.28, 0.48, 0.12, 0.16], true)

大斜は sqrt((a - c)^2 + 4r1^2) = 1.275617497528158 である。

「術」は「1.6275 の平方根に小斜を掛ける」ということなのであるが,パラメータは正確に小数点以下 2 桁なので(100倍すれば整数)なので,(a - c)^2 + 4r1^2 =1.6272 である。
末尾の五は二の誤記または誤読ではないか?

(a, c, r1, r2, x2) = res[1]
length = sqrt((a - c)^2 + 4r1^2)

    1.275617497528158

(a - c)^2 + 4r1^2

    1.6272000000000002

function draw(r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (a, c, r1, r2, x2) = res[1] #[2.32023, 0.56338, 1, 0.24517, 0.32809]
    b = 2r1
    plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5)
    circle(a - r1, r1, r1)
    circle(x2, r2, r2, :blue)
    circle(r2, b - r2, r2, :blue)
    segment(0, 0, c, b, :magenta)
    segment(c, b, 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)
        vline!([0], color=:gray80, lw=0.5)
        point(a - r1, r1, "大円:r1,(a-r1,r1) ", :red, :right, :vcenter)
        point(x2, r2, "小円:r2\n(x2,r2)", :blue, :center, delta=-delta)
        point(r2, b - r2, "小円:r2\n(r2,b-r2)", :blue, :center, delta=-delta)
        point(c, b, "c", :magenta, :center, :bottom, delta=delta/2)
        point(a, b, "(a,b)", :green, :right, :bottom, delta=delta/2)
    end
end;

draw(1/2, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1129)改訂版

2025年01月10日 | Julia

算額(その1129)改訂版

算額(その1129)は依拠した図が間違っているので,改訂版を書いた。

四十五 岩手県一関市真滝 熊野白山滝神社 明治13年(1880)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03061
https://w.atwiki.jp/sangaku/pages/327.html

キーワード:円4個,長方形,斜線2本
#Julia, #SymPy, #算額, #和算

山村の図では,条件不足で解が一意に定まらない。「今有如図」の図に基づいて得られた解は「答」,「術」と一致する。

長方形の中に 2 本の斜線を引き,区画された領域に大円 1 個,中円 1 個,小円 2 個を容れる。小円の直径が 1 寸のとき,大円の直径はいかほどか。

長方形の長辺,短辺を a, b
斜線と長方形の辺の交点座標を (c, 0), (a, d)
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (a - r2, r2)
小円の半径と中心座標を r3, (r3, y31), (a - r3, y32)
とおき,以下の連立方程式を解く。

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

using SymPy
@syms a::positive, b::positive, c::positive, d::positive,
      r1::positive, r2::positive,
      r3::positive, y31::positive, y32::positive
eq1 = dist(0, b, a, d, r1, r1) - r1^2
eq2 = dist(0, b, a, d, a - r2, r2) - r2^2
eq3 = dist(0, b, a, d, r3, y31) - r3^2
eq4 = dist(0, b, a, d, a - r3, y32) - r3^2
eq5 = dist(c, 0, a, b, r1, r1) - r1^2
eq6 = dist(c, 0, a, b, a - r2, r2) - r2^2
eq7 = dist(c, 0, a, b, a - r3, y32) - r3^2
eq8 = (r1 - r3)^2 + (y31 - r1)^2 - (r1 + r3)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (a, b, c, d, r1, r2, y31, y32))

function H(u)
    (a, b, c, d, r1, r2, y31, y32) = u
    return [
        dist(0, b, a, d, r1, r1) - r1^2,
        dist(0, b, a, d, a - r2, r2) - r2^2,
        dist(0, b, a, d, r3, y31) - r3^2,
        dist(0, b, a, d, a - r3, y32) - r3^2,
        dist(c, 0, a, b, r1, r1) - r1^2,
        dist(c, 0, a, b, a - r2, r2) - r2^2,
        dist(c, 0, a, b, a - r3, y32) - r3^2,
        (r1 - r3)^2 + (y31 - r1)^2 - (r1 + r3)^2
    ]
end;

r3 = 1/2
iniv = BigFloat[4.6, 4.1, 2.4, 1.4, 1.5, 0.8, 3.2, 2.2]
res = nls(H, ini=iniv)

    ([4.732050807568878, 4.098076211353316, 2.366025403784439, 1.3660254037844386, 1.5, 0.8660254037844386, 3.232050807568877, 2.232050807568877], true)

小円の直径が 1 のとき,大円の直径は 3, 中円の直径は 1.73205 である。

全てのパラメータは,以下のとおり。

r3 = 0.5;  a = 4.73205;  b = 4.09808;  c = 2.36603;  d = 1.36603;  r1 = 1.5;  r2 = 0.866025;  y31 = 3.23205;  y32 = 2.23205

function draw(r3, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (a, b, c, d, r1, r2, y31, y32) = res[1]
    @printf("小円の直径が %g のとき,大円の直径は %g, 中円の直径は %g である。\n", 2r3, 2r1, 2r2)
    @printf("r3 = %g;  a = %g;  b = %g;  c = %g;  d = %g;  r1 = %g;  r2 = %g;  y31 = %g;  y32 = %g\n", r3, a, b, c, d, r1, r2, y31, y32)
    plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5)
    circle(r1, r1, r1)
    circle(a - r2, r2, r2, :blue)
    circle(r3, y31, r3, :magenta)
    circle(a - r3, y32, r3, :magenta)
    segment(0, b, a, d, :orange)
    segment(c, 0, a, b, :orange)
    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, b, "(a,b)", :green, :right, :bottom, delta=delta/2)
        point(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta)
        point(a - r2, r2, "中円:r2,(a-r2,r2)", :blue, :center, delta=-delta)
        point(r3, y31, "小円:r3\n(r3,y31)", :magenta, :center, delta=-delta)
        point(a - r3, y32, "小円:r3\n(a-r3,y32)", :magenta, :center, delta=-delta)
    end
end;

draw(1/2, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村