裏 RjpWiki

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

算額(その1404)改訂版

2025年03月06日 | Julia

算額(その1404) 改訂版

算額(その1404)は,依拠した図がでたらめなものであったので,改訂版を書いた。

三十二 一関市舞川 観福寺内地蔵堂 明治43年(1901)

山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

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


キーワード:円2個,直角三角形4個,等脚台形
#Julia, #SymPy, #算額, #和算

問の原文は以下のとおりである。
「今有如図半円内設等円二個従其親所洩二斜容等円一個其等円径若干問得半円径術如何」

山村の図は正しいものではない。「今有如図」の図に従って解を求めると,術に一致する解が得られた。

円弧の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, r), (0 - r)
円弧と等円の接点座標を (x0, y0)
とおき,以下の連立方程式を解く。

斜線は,(x0, y0), (x, 0), (0, y01) を通る。

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

@syms R::positive, r::positive, x::positive, x0::positive, y0::positive
eq1 = dist2(x, 0, x0, y0, 0, -r, r)
eq2 = x^2 + r^2 - (R - r)^2
eq3 = (x0 - x)^2 + (y0 - r)^2 - r^2
eq4 = x0^2 + y0^2 - R^2
res = solve([eq1, eq2, eq3, eq4], (R, x, x0, y0))[1]

    (3*r, sqrt(3)*r, 3*sqrt(3)*r/2, 3*r/2)

円弧の半径 R は,等円の半径 r の 3 倍である。
等円の直径が 1 寸のとき,円弧直径は 3 寸である。
x = 0.866025,x0 = 1.29904, y0 = 0.75

function draw(r, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, x, x0, y0) = (3*r, sqrt(3)*r, 3*sqrt(3)*r/2, 3*r/2)
    @printf("等円の直径が %g のとき,円弧の直径は %g である。ただし,x = %g,x0 = %g, y0 = %g である。\n", 2r, 2R, x, x0, y0)
    plot()
    circle(0, 0, R, beginangle=0, endangle=180)
    circle2(x, r, r, :blue)
    circle(0, -r, r, :blue)
    (x01, y01) = float.(intersection(x0, y0, x, 0, -x0, y0, -x, 0))
    println((x01, y01))
    segment(x0, y0, 0, y01)
    segment(-x0, y0, 0, y01)
    segment(-R, 0, R, 0, :red)
    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(x0, y0, "(x0,y0)", :blue, :right, :vcenter, deltax=-delta)
        point(x, 0, "x", :blue, :center, delta=-delta)
        point(R, 0, "R", :red, :center, delta=-delta)
        point(x, r, "等円:r,(x,r)", :blue, :center, delta=-delta, deltax=-2delta)
        point(0, -r, "等円:r,(0,-r)", :blue, :center, delta=-delta, deltax=-2delta)
        point(0, y01, "y01", :black, :left, :vcenter, deltax=delta)
    end
end;

draw(1/2, true)

 

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

算額(その1403) 改訂版

2025年03月06日 | Julia

算額(その1403) 改訂版

算額(その1403)は,依拠した図がでたらめなものであったので,改訂版を書いた。

三十二 一関市舞川 観福寺内地蔵堂後額 明治34年(1901)

山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03086
https://w.atwiki.jp/sangaku/pages/158.html

キーワード:円2個,等脚台形,斜線
#Julia, #SymPy, #算額, #和算, #数学

合同な直角三角形を4個組み合わせて等脚台形を作り,大円 1 個,小円 1 個を容れる。小円の直径が与えられたときに大円の直径を得る術を述べよ。

山村の図は正しいものではない。「今有如図」の図に従って解を求めると,術に一致する解が得られた。

直角三角形の直角を挟む 2 辺の短い方,長い方の長さをそれぞれ a, h
大円の半径と中心座標を r1, (2a - r1, r1)
小円の半径と中心座標を r2, (a + r2, h - r2)
とおき,以下の連立方程式を解く。

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

using SymPy

@syms a::positive, h::positive,
      r1::positive, r2::positive
eq1 = dist2(0, 0, a, h, 2a - r1, r1, r1)
eq2 = dist2(3a, 0, 2a, h, a + r2, h - r2, r2)
eq3 = ((a + r2) - (2a - r1))^2 + (h - r2 - r1)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (r1, a, h))[5]

    (-5*r2 + sqrt(41)*r2, -3*r2/2 + sqrt(41)*r2/2, r2*(3 + sqrt(41))/2)

大円の半径 r1 は,小円の半径 r2 の (√41 - 5) = 1.4031242374328485 倍である。
小円の直径が 1 寸のとき,大円の直径は 1.4031242374328485 寸である。

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, a, h) = (-5*r2 + sqrt(41)*r2, -3*r2/2 + sqrt(41)*r2/2, r2*(3 + sqrt(41))/2)
    plot([0, 3a, 2a, a, 0], [0, 0, h, h, 0], color=:green, lw=0.5)
    plot!([a, a, 2a, 2a], [0, h, 0, h], color=:green, lw=0.5)
    circle(2a - r1, r1, r1)
    circle(a + r2, h - r2, r2, :blue)
    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(2a - r1, r1, "大円:r1\n(2a-r1,r1)", :red, :center, delta=-delta, deltax=2delta)
        point(a + r2, h - r2, "小円:r2\n(a+r2,h-r2)", :blue, :center, delta=-delta)
        point(a, h, "(a,h)", :green, :center, :bottom, delta=delta)
        point(2a, h, "(2a,h)", :green, :center, :bottom, delta=delta)
        point(3a, 0, "3a", :green, :right, :bottom, delta=delta/2, deltax=-delta)
    end
end;

draw(1/2, true)

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

算額(その1647)

2025年03月06日 | Julia

算額(その1647)

三四 武蔵国埼玉郡下忍村遍照院境内 金毘羅社(神楽堂) 天保11年(1840)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.
キーワード:円7個,楕円3個,正方形
#Julia, #Julia, #SymPy, #算額, #和算, #数学

正方形の中に,交差する大楕円 2 個を設け,四隅に丙円 4 個,中央部に小楕円,更にその中に甲円 1 個,乙円 2 個を容れる。乙円とは丙円はそれぞれの楕円の長径端において楕円と 1 点で接する最大の円である(曲率円)。小楕円の短径が 1 寸のとき,大楕円の短径はいかほどか。

大楕円の長半径,短半径,中心座標を a, b, (0, 0)
小楕円の長半径,短半径,中心座標を c, d, (0, 0)
甲円の半径と中心座標を r1, (0, 0); r1 = b
乙円の半径と中心座標を r2, (b - r2, 0)
丙円の半径と中心座標を r4, ((c - r4)/√2, (c - r4)/√2)
大楕円と小楕円の接点座標を A:(x0, y0)
接点から長軸へおろした垂線の脚を B,原点を O としたとき AB = y1, 0B = x1
とおき,以下の連立方程式を解く。

小楕円の中の乙円は曲率円なので,r2 = b^2/a
また,a = b + 2r2 である。

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

using SymPy

@syms a::positive, b::positive, r2::positive,
      c::positive, d::positive, r4::positive,
      x0::positive, y0::positive, x1::positive, y1::positive

eq01 = r2 - b^2/a
eq02 = a - (b + 2r2)
solve([eq01, eq02], (a, r2))[1]

    (2*b, b/2)

a = 2b, r2 = b/2 である。

eq03 = r4 - d^2/c
eq04 = c - (d + 2r4)
solve([eq03, eq04], (c, r4))[1]

    (2*d, d/2)

c = 2d, r4 = d/2 である。

大楕円,小楕円は相似である。

y1 = (x0 + y0)/√Sym(2)
eq0 = y0^2 - (4b^2 - x0^2)/4
eq4 = (x0^2 + y0^2) - ((x0 + y0)^2/2 + x1^2)
eq5 = 3√Sym(2)x0*y0/(4y0 + x0) - 3x1/4
eq6 = x1^2 + 4y1^2 - 4d^2
res = solve([eq0, eq4, eq5, eq6], (d, x0, y0, x1))[2]  # 2 of 2

    (b*sqrt(3*sqrt(41) + 25)/4, b*sqrt(205 - sqrt(41))*(41*sqrt(2) + 13*sqrt(82))/1312, b*sqrt(205 - sqrt(41))*(-3*sqrt(82) + 41*sqrt(2))/1312, b*sqrt(8405 - 41*sqrt(41))/82)

# d
res[1] |> println
res[1].evalf() |> println

    b*sqrt(3*sqrt(41) + 25)/4
    1.66225322815709*b

大楕円の短半径 d は,小楕円の短半径 b の sqrt(3√41+ 25)/4 = 1.66225322815709 倍」である。
小楕円の短径が 1 寸のとき,大楕円の短径は 1.66225322815709 寸である。

術は「大楕円の短径 d は,小楕円の短径 b の sqrt(√369 + 25)/4 倍」である。
sqrt(√369 + 25)/4 = sqrt(3√41 + 25)/4 なので,同じである。

function draw(b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (d, x0, y0, x1) = (b*sqrt(3*sqrt(41) + 25)/4,
        b*sqrt(205 - sqrt(41))*(41*sqrt(2) + 13*sqrt(82))/1312,
        b*sqrt(205 - sqrt(41))*(-3*sqrt(82) + 41*sqrt(2))/1312,
        b*sqrt(8405 - 41*sqrt(41))/82)
    f = d/b
    a = 2b
    c = 2d
    r1 = b
    @printf("factor = %g\n", f)
    @printf("a = %g;  b = %g;  c = %g;  d = %g\n", a, b, c, d)
    r2 = b^2/a  # 曲率円
    r3 = r2*f
    r4 = d^2/c  # 曲率円
    s = sqrt(2(c^2 + d^2))  # 算法助術−94
    plot(s/2 .*[1, 1, -1, -1, 1], s/2 .*[-1, 1, 1, -1, -1], color=:orange, lw=0.5)
    ellipse(0, 0, a, b, color=:red)
    ellipse(0, 0, a, b, color=:pink, φ=90)
    circle(0, 0, b, :green)
    circle2(a - r2, 0, r2, :brown)
    ellipse(0, 0, c, d, color=:blue, lw=0.5, φ=45)
    ellipse(0, 0, c, d, color=:blue, lw=0.5, φ=135)
    ellipse(0, 0, c, d, color=:gray80, lw=0.5, φ=0)
    circle(0, 0, d, :lightgreen)
    circle4((c - r4)/√2, (c - r4)/√2, r4, :deeppink)
    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, "甲円:r1,(0,0)", :green, :center, :bottom, delta=2delta)
        point(a - r2, 0, "乙円:r2\n(a-r2,0)", :brown, :center, delta=-delta/2)
        point((c - r4)/√2, (c - r4)/√2, "丙円:r4\n((c-r4)/√2,\n(c-r4)/√2)", :magenta, :center, delta=-delta/2)
        point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta, deltax=-delta/2)
        point(0, b, "b", :red, :center, :bottom, delta=delta/2)
        point(0, d, "d", :red, :center, :bottom, delta=delta/2)
        point(a, 0, "a", :red, :right, :vcenter, deltax=-delta/2)
        point(c, 0, "c", :red, :right, :vcenter, deltax=-delta/2)
        segment(0, 0, -s/2, -s/2, :gray80)
        point(0, 0, "O", :black, :left, delta=-delta)
        segment(0, 0, -x0, y0)
        point(-x0, y0, "A", :black, :right, :vcenter, deltax=-delta)
        segment(0, 0, -x1/√2, -x1/√2)
        point(-x1/√2, -x1/√2, "B", :black, :right, :vcenter, deltax=-delta)
        segment(-x1/√2, -x1/√2, -x0, y0)
    end
end;

draw(1, true)

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

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

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