裏 RjpWiki

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

算額(その1478)

2024年12月17日 | Julia

算額(その1478)

千葉県印西市 竜湖寺 文久元年(1861)
山根誠司:算法勝負!「江戸の数学」に挑戦,ブルーバックス,講談社,東京都,2015年。

キーワード:折り紙,折り鶴
#Julia, #SymPy, #算額, #和算

折り鶴の羽の横幅と長さの比を求めよ。

折り紙の一辺の長さを a とおく。
一旦折った鶴を開いて,折り目を図に示す。

左上が尾(頭),右上が羽の先端である。羽の横幅は DE,長さは CF である。
y1, y2 は左上の頂点を通り,傾きが -tan(π/4)*(3/8) と -tan(π/4)*(1/4) の直線の y 切片である。 

include("julia-source.txt");

using SymPy
@syms a, θ1, θ2, x, y
θ1 = 90*(Sym(3)/8)
eq1 = (y - a) + tand(θ1)*(x + a)
y1 = solve(eq1(x => 0), y)[1] |> simplify
y1 |> println

    a*(-2*sqrt(sqrt(2) + 2) + sqrt(2) + sqrt(2*sqrt(2) + 4))

幅 DE は y1*√2 である。

width = y1*sqrt(Sym(2))
width |> println
width(a => 10).evalf() |> println

    sqrt(2)*a*(-2*sqrt(sqrt(2) + 2) + sqrt(2) + sqrt(2*sqrt(2) + 4))
    4.69266270539641

長さは CF = OF - OA - 2AB である。

θ2 = 90*(Sym(1)/4)
eq2 = (y - a) + tand(θ2)*(x + a)
y2 = solve(eq2(x => 0), y)[1] |> simplify
y2 |> println

    a*(2 - sqrt(2))

AB = (y2/2 - y1/2)*sqrt(Sym(2));

OA = y1/2*sqrt(Sym(2));

length = sqrt(Sym(2))*a - OA - 2AB |> simplify
length |> println

    a*(-sqrt(2*sqrt(2) + 4) - sqrt(2) + sqrt(sqrt(2) + 2) + 3)

length(a => 10).evalf() |> println

    8.20419572896725

ratio = length/width;

@syms d
apart(ratio, d) |> simplify |> println

    -sqrt(2*sqrt(2) + 4)/2 + 1/2 + sqrt(2)/2 + sqrt(sqrt(2) + 2)

ratio.evalf() |> println

    1.74830288133274

羽根の幅は sqrt(2)*a*(-2*sqrt(sqrt(2) + 2) + sqrt(2) + sqrt(2*sqrt(2) + 4))
羽の長さは a*(-sqrt(2*sqrt(2) + 4) - sqrt(2) + sqrt(sqrt(2) + 2) + 3)
長さ/幅は -sqrt(2*sqrt(2) + 4)/2 + 1/2 + sqrt(2)/2 + sqrt(sqrt(2) + 2)

折り紙の一辺の長さを 10cm とすれば,幅は 4.69266270539641,長さは 8.20419572896725 で,長さは幅の 1.74830288133274 倍である。

function draw(a, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    θ = 90*(3/8)
    y1= a*(-2*sqrt(sqrt(2) + 2) + sqrt(2) + sqrt(2*sqrt(2) + 4))
    println("y1 = $y1")
    width = sqrt(2*y1^2)
    println("width = $width")
    y2 = a*(2 - sqrt(2))

    plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:blue, lw=0.5)
    segment(-a, a, 0, y1, :blue)
    segment(-a, a, 0, y2, :red)
    segment(0, 0, a, a, :gray70)
    segment(0, y1, y1, 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(0, y1, "y1 ", :blue, :right)
        point(y1, 0, "y1 ", :blue, :right)
        point(0, y2, "y2 ", :red, :right)
        point(y2, 0, "y2 ", :red, :right)
        point(y1/2, y1/2, "A")
        point(y2/2, y2/2, "B")
    end
end;

draw(10, true)

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

算額(その1477)

2024年12月17日 | Julia

算額(その1477)

四十三 岩手県一関市真滝 熊野白山滝神社 弘化3年(1846)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html

キーワード:円5個,直角三角形
#Julia, #SymPy, #算額, #和算

外円の中に,直角三角形と,その 3 辺の中点を中心とする 3 個の円を容れる。また,直角三角形に内接する等円を容れる。等円の直径が与えられたとき,外円の直径はいかほどか。

しかし,設問には難点がある。
1. 鈎円,股円,弦円は鈎,股,弦を直径とする円とされているが,図では股円は股を直径とする円には見えない。
2. 弦円の中にも等円が描かれている。しかし,弦円に内接しているように見えるが,弦に接しているわけでもなく,「そこにある」という意味しか持っていない。
3. これは重大な欠点であるが,等円の直径,外円の直径はともに鈎,股により決まるがその比は一定ではないということである。

本問は,算額(その1230)https://blog.goo.ne.jp/r-de-r/e/af96a9c20cd935b40c76e49b5c6a4826 の類題である。鈎,股が特定の値を取るとき外円の直径が,直角三角形に内接する等円の直径の何倍かと問うものである。

外円の半径と中心座標を R, (x, y)
鈎,股,弦を a, b, c
それぞれの中点を中心とする円の半径を r1, r2, r3
とおき,以下の連立方程式を解く。
直角三角形に内接する等円の半径は (a + b - c)/2 である。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive, R::positive, r1::positive, r2::positive, r3::positive, x::positive, y::positive
@syms a, b, c, R, r1, r2, r3, x, y
#c = sqrt(a^2 + b^2)
r1 = a/2
r2 = b/2
r3 = c/2
eq1 = x^2 + (r1 - y)^2 - (R - r1)^2
eq2 = (r2 - x)^2 + y^2 - (R - r2)^2
eq3 = (r2 - x)^2 + (r1 - y)^2 - (R - r3)^2
res = solve([eq1, eq2, eq3], (R, x, y))[1];  # 1 of 2

# R
ans_R = res[1](c^2 => a^2 + b^2)(sqrt(a^2*b^4) => a*b^2)
ans_R |> println
ans_R(a => 3, b => 4, c => 5) |> println

    (a^5 - a^4*c - 3*a^3*b^2 - a^3*(a^2 + b^2) - 3*a^2*b^3 + 2*a^2*b^2*c + a^2*c^3 + b^5 - b^4*c - b^3*(a^2 + b^2) + b^2*c^3)/(4*(a^4 - 2*a^3*c - a^2*b^2 + a^2*(a^2 + b^2) + b^4 - 2*b^3*c + b^2*(a^2 + b^2)))
    72/23

分子分母を別々に因子分解して再度分数にすることにより,若干簡潔になる。

num = numerator(ans_R) |> factor
den = denominator(ans_R) |> factor;
ans_R2 = num/den
ans_R2 |> println
ans_R2(a => 3, b => 4, c => 5) |> println

    (-a^4*c - 4*a^3*b^2 - 4*a^2*b^3 + 2*a^2*b^2*c + a^2*c^3 - b^4*c + b^2*c^3)/(4*(2*a^4 - 2*a^3*c + a^2*b^2 + 2*b^4 - 2*b^3*c))
    72/23

鈎,股が 3 寸,4 寸のとき,外円の直径は 2*72/23 = 6 + 6/23 である。
等円の直径は 3 + 4 - 5 = 2 で,外円の直径は等円の直径の 144/23/2 = 72/23 = 3.130434782608696 倍である。

術は「等円の直径を 105/34 倍すれば,全円の直径になる」とでたらめなことを書いている。そもそも鈎股の比率が異なる図形は相似ではないので,単純に等円の何倍ということはない。たとえば,鈎 = 3,股 = 6 のとき,外円の直径は 4.07571052090885,等円の直径は 1.1458980337503153 となり,比は 3.556782890681637 である。

@syms d
ratio = apart(2ans_R2/(a + b - c), d)
ratio |> println

    -(a^4*c + 4*a^3*b^2 + 4*a^2*b^3 - 2*a^2*b^2*c - a^2*c^3 + b^4*c - b^2*c^3)/(4*a^5 + 4*a^4*b - 8*a^4*c + 2*a^3*b^2 - 4*a^3*b*c + 4*a^3*c^2 + 2*a^2*b^3 - 2*a^2*b^2*c + 4*a*b^4 - 4*a*b^3*c + 4*b^5 - 8*b^4*c + 4*b^3*c^2)

ratio(a => 3, b => 4, c => 5).evalf() |> println
ratio(a => 3, b => 6, c => √45).evalf() |> println

    3.13043478260870
    3.55678289068168

ans_R2(a => 3, b => 6, c => √45) |> println
r4 = (3 + 6 - √45)/2 |> println
4.07571052090885/1.1458980337503153 |> println

    4.07571052090886
    1.1458980337503153
    3.556782890681637

x, y については,数値解のみ述べる。

res[2](a => 3, b => 4, c => 5) |> println 
res[3](a => 3, b => 4, c => 5) |> println 

    36/23
    24/23

function draw(a, b, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    c = sqrt(a^2 + b^2)
    r4 = (a + b - c)/2
    (R, x, y) = ((a^5 - a^4*c - a^3*b^2 - a^3*c^2 - a^2*b^3 + a^2*c^3 - sqrt(2)*a^2*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - sqrt(2)*a*b*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + sqrt(2)*a*c*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + b^5 - b^4*c - b^3*c^2 + b^2*c^3)/(4*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)), (a^4*b^2 + a^3*b^3 - 3*a^3*b^2*c + sqrt(2)*a^3*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - 3*a^2*b^3*c + 3*a^2*b^2*c^2 + sqrt(2)*a^2*b*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - 2*sqrt(2)*a^2*c*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) - a*b^5 + a*b^4*c + a*b^3*c^2 - a*b^2*c^3 - sqrt(2)*a*b*c*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + sqrt(2)*a*c^2*sqrt(a*b^3*(-a^2 + 2*a*b - b^2 + c^2)) + b^6 - b^5*c - b^4*c^2 + b^3*c^3)/(4*b*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)), a*(a^4 - a^3*b - a^3*c + a^2*b*c - a^2*c^2 + a*b^3 - 3*a*b^2*c + a*b*c^2 + a*c^3 + b^4 - 3*b^3*c + 3*b^2*c^2 - b*c^3)/(4*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)) + sqrt(2)*sqrt(-a*b^3*(a - b - c)*(a - b + c))*(b - c)*(a + b - c)/(4*(a^4 - 2*a^3*c - a^2*b^2 + a^2*c^2 + b^4 - 2*b^3*c + b^2*c^2)))
    @printf("外円は等円の %g 倍大きい\n", R/r4)
    @printf("a = %g;  b = %g;  c = %g;  R = %g;  x = %g;  y = %g\n", a, b, c, R, x, y)
    plot([0, b, 0, 0], [0, 0, a, 0], color=:green, lw=0.5)
    circle(x, y, R)
    circle(0, a/2, a/2, :blue)
    circle(b/2, 0, b/2, :magenta)
    circle(b/2, a/2, c/2, :orange)
    circle(r4, r4, r4, :tomato)
    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(x, y, " 外円:R,(x,y)", :red, :left, :vcenter)
        point(0, a/2, "鈎円:r1,(0,a/2)", :blue, :center, delta=-delta/2)
        point(b/2, 0, "股円:r2,(b/2,0)", :magenta, :center, delta=-delta/2)
        point(b/2, a/2, "弦円:r3,(b/2,a/2)", :orange, :center, delta=-delta/2)
        point(r4, r4, "等円:r4,(r4,r4)", :tomato, :center, delta=-delta/2)
        point(0, a, "a ", :green, :right, :bottom, delta=delta)
        point(b, 0, " b", :green, :left, delta=-delta)
    end
end;

draw(3, 4, true)

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

算額(その1476)

2024年12月17日 | Julia

算額(その1476)

三十九 一関市前堀 前堀熊野神社 明治37年(1904)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html
キーワード:円3個,半円5個
#Julia, #SymPy, #算額, #和算

半円の中に 4 個の半円を描き,隙間に甲円 2 個,乙円 1 個を容れる。甲円の直径が与えられたとき,乙円の直径を求める術を述べよ。

一番小さい半円の半径を r0
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (0, y2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, x1::positive, y1::positive, r2::positive, y2::positive
eq1 = r0^2 + y2^2 - (r0 + r2)^2
eq2 = (r0/2)^2 + y2^2 - (3r0/2 - r2)^2
eq3 = (r0 - x1)^2 + y1^2 - (r1 + r0)^2
eq4 = (x1 + r0/2)^2 + y1^2 - (3r0/2 + r1)^2
eq5 = (x1 - r0/2)^2 + y1^2 - (3r0/2 - r1)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r0, x1, y1, r2, y2))[1]

    (4*r1, 3*r1, 2*sqrt(6)*r1, 8*r1/5, 8*sqrt(6)*r1/5)

乙円の半径は甲円の半径の 8/5 倍である。

function draw(r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r0, x1, y1, r2, y2) = (4*r1, 3*r1, 2*sqrt(6)*r1, 8*r1/5, 8*sqrt(6)*r1/5)
    @printf("甲円の直径が %g のとき,乙円の直径は %g である。\n", 2r1, 2r2)
    plot()
    circle(0, 0, 2r0, :green, beginangle=0, endangle=180)
    circle(r0/2, 0, 3r0/2, beginangle=0, endangle=180)
    circle(-r0/2, 0, 3r0/2, beginangle=0, endangle=180)
    circle(r0, 0, r0, :blue, beginangle=0, endangle=180)
    circle(-r0, 0, r0, :blue, beginangle=0, endangle=180)
    circle2(x1, y1, r1, :orange)
    circle(0, y2, r2, :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(x1, y1, "甲円:r1,(x1,y1)", :black, :center, delta=-delta, deltax=8delta)
        point(0, y2, "乙円:r2\n(0,y2)", :magenta, :center, delta=-delta)
        point(r0, 0, "r0", :blue, :left, :bottom, delta=delta)
        point(r0/2, 0, "r0/2", :red, :left, :bottom, delta=delta)
    end
end;

draw(1/2, true)

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

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

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