裏 RjpWiki

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

算額(その1528)

2025年01月09日 | Julia

算額(その1528)

三十三 岩手県一関市舞川相川 菅原神社 嘉永3年(1850)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

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

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

直角三角形の鋭角の頂点と対辺の中点を結ぶ斜線を 2 本引き,円を容れる。直角三角形の直角を挟む 2 辺の長さが 3 寸,4 寸のとき,円の直径はいかほどか。

鈎,股をそのまま「鈎」,「股」
円の半径と中心座標を r, (x, y)
とおき,以下の連立方程式を解く。
一度には解けないので,順次解いていく。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive,
      r::positive, x::positive, y::positive
eq1 = dist2(0, 鈎, 股, 0, x, y, r)
eq2 = dist2(0, 鈎, 股/2, 0, x, y, r)
eq3 = dist2(0, 鈎/2, 股, 0, x, y, r);
# res = solve([eq1, eq2, eq3], (r, x, y))

ans_x = solve(eq1, x)[1];  # 1 of 2

# y
eq12 = eq2(x => ans_x) |> simplify
ans_y = solve(eq12, y)[1]  # 1 of 2
ans_y |> println

    (-2*r*sqrt(股^2 + 鈎^2) - r*sqrt(股^2 + 4*鈎^2) + 股*鈎)/股

# x
ans_x = ans_x(y => ans_y)
ans_x |> println

    (-r*sqrt(股^2 + 鈎^2) - 股*(-鈎 + (-2*r*sqrt(股^2 + 鈎^2) - r*sqrt(股^2 + 4*鈎^2) + 股*鈎)/股))/鈎

# r
eq13 = eq3(x => ans_x)(y => ans_y)
ans_r = solve(eq13, r)[1]  # 1 of 2
ans_r |> println

    股*鈎*(3*sqrt(股^2 + 鈎^2) + sqrt(股^2 + 4*鈎^2) - sqrt(4*股^2 + 鈎^2))/(6*(股^2 + 2*鈎^2 + sqrt(股^2 + 鈎^2)*sqrt(股^2 + 4*鈎^2)))

2ans_r(鈎 => 3, 股 => 4).evalf() |> println

    0.780358219829331

鈎,股が 3 寸,4 寸のとき,円の直径は 0.780358219829331 寸である。

プログラム的には以下のようにすれば,簡潔に求めることができる。

function rxy(鈎, 股)
    s = sqrt(股^2 + 4*鈎^2)
    弦 = sqrt(股^2 + 鈎^2)
    r = 股*鈎*(3弦 + s - sqrt(4股^2 + 鈎^2))/(6股^2 + 12鈎^2 + 6弦*s)
    y = (股*鈎 - 2r*弦 - r*s)/股
    x = (股*(鈎 - y)- r*弦)/鈎
    return (r, x, y)
end;
rxy(3, 4)

    (0.3901791099146656, 1.5881723747992602, 1.3211468315072228)

function draw(鈎, 股, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r, x, y) = rxy(鈎, 股)
    @printf("r = %g;  x = %g;  y = %g\n", r, x, y)
    plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
    circle(x, y, r)
    segment(0, 鈎/2, 股, 0, :blue)
    segment(0, 鈎, 股/2, 0, :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(股, 0, "股", :green, :left, :bottom, delta=delta)
        point(0, 鈎, " 鈎", :green, :left, :bottom, delta=delta)
        point(股/2, 0, "股/2", :green, :left, :bottom, delta=delta)
        point(0, 鈎/2, " 鈎/2", :green, :left, :bottom, delta=delta)
        point(x, y, "(x,y)", :red, :center, delta=-delta)
    end
end;

draw(3, 4, true)

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

算額(その1527)

2025年01月09日 | Julia

算額(その1527)

三十三 岩手県一関市舞川相川 菅原神社 嘉永3年(1850)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

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

キーワード:円5個,二等辺三角形,圭,弦2本
#Julia, #SymPy, #算額, #和算

全円の中に,二等辺三角形 1 個,大円 2 個,小円 2 個,斜線 2 本(圭)を容れる。小円の直径が 1 寸のとき,大円の直径はいかほどか。

山村の図は不正確である。

全円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, R - r1), (0, r1 - R)
小円の半径と中心座標を r2, (x2, 2r1 - R + r2)
斜線と全円の交点座標を (x, -sqrt(R^2 - x^2))
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive,
      x2::positive, x::positive
eq1 = x2^2 + (2R - 3r1 - r2)^2 - (r1 + r2)^2
eq2 = dist(0, R, sqrt(R^2 - (2r1- R)^2), 2r1 - R, x2, 2r1 - R + r2) - r2^2
eq3 = dist(0, R, x, -sqrt(R^2 - x^2), x2, 2r1 - R + r2) - r2^2
eq4 = dist(0, R, x, -sqrt(R^2 - x^2), 0, r1 - R) - r1^2;
# res = solve([eq1, eq2, eq3, eq4], (R, r1, x2, x))

function H(u)
    (R, r1, x2, x) = u
    return [
        x2^2 + (2R - 3r1 - r2)^2 - (r1 + r2)^2,
        dist(0, R, sqrt(R^2 - (2r1- R)^2), 2r1 - R, x2, 2r1 - R + r2) - r2^2,
        dist(0, R, x, -sqrt(R^2 - x^2), x2, 2r1 - R + r2) - r2^2,
        dist(0, R, x, -sqrt(R^2 - x^2), 0, r1 - R) - r1^2
    ]
end;
r2 = 1/2
iniv = BigFloat[4.5, 2.0, 2.23, 2.46] .* r2
res = nls(H, ini=iniv)

    ([2.25, 1.0, 1.118033988749895, 1.2321190896427412], true)

大円の半径は小円の半径の 2 倍である。ちなみに全円の半径は 4.5 倍である。

function draw(r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, r1, x2, x) = r2.*[4.5, 2.0, 70711162//31622993, 510056828//206983575]
    @printf("小円の直径が %g のとき,大円の直径は %g,全円の直径は %g である。\n", 2r2, 2r1, 2R)
    plot([0, sqrt(R^2 - (2r1 - R)^2), -sqrt(R^2 - (2r1 - R)^2), 0], [R, 2r1 - R, 2r1 - R, R], color=:green, lw=0.5)
    plot!([-x, 0, x], [-sqrt(R^2 - x^2), R, -sqrt(R^2 - x^2)], color=:orange, lw=0.5)
    circle(0, 0, R, :magenta)
    circle22(0, R - r1, r1)
    circle2(x2, 2r1 - R + 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(0, R, "R", :magenta, :center, :bottom, delta=delta/2)
        point(0, R - r1, "大円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
        point(0, r1 - R, "大円:r1,(0,r1-R)", :red, :center, delta=-delta/2)
        point(x2, 2r1 - R + r2, "小円:r2\n(x2,2r1-R+r2)", :blue, :center, :bottom, delta=delta/2)
        point(x, -sqrt(R^2 - x^2), "(x,-sqrt(R^2-x^2))", :magenta, :left, delta=-2delta, deltax=-3delta)
        point(0, 2r1 - R, "2r1-R", :green, :center, :bottom, delta=delta)
    end
end;

draw(1/2, true)

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

算額(その1526)

2025年01月09日 | Julia

算額(その1526)

三十三 岩手県一関市舞川相川 菅原神社 嘉永3年(1850)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

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

キーワード:円1個,直角三角形,正方形,斜線
#Julia, #SymPy, #算額, #和算

直角三角形の中に直角三角形と斜線を 1 本容れ,隙間に円を容れる。直角三角形の直角を挟む二辺(鈎,股)が 3 寸,4 寸のとき,円の直径はいかほどか。

鈎,股をそのまま「鈎」,「股」
正方形の一辺の長さを a
円の半径と中心座標を r, (x, a + r)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, a::positive,
      r::positive, x::positive
eq1 = dist2(0, 鈎, a, 0, x, a + r, r)
eq2 = dist2(0, 鈎, 股, 0, x, a + r, r)
eq3 = (鈎 - a)/a - 鈎/股;

一度には解けないので,まず a, x を求め,次に r を求める。

res = solve([eq3, eq1], (a, x))[2]

    (股*鈎/(股 + 鈎), (-r*股^2 - r*股*鈎 + r*股*sqrt(2*股^2 + 2*股*鈎 + 鈎^2) + r*鈎*sqrt(2*股^2 + 2*股*鈎 + 鈎^2) + 股*鈎^2)/(股^2 + 2*股*鈎 + 鈎^2))

eq2 に a, x の解を代入し,その式を r について解く。

eq12 = eq2(a => res[1], x => res[2]) |> simplify |> numerator;
res = solve(eq12, r)[1];  # 1 of 2

@syms d
ans_r = apart(res, d) |> simplify
ans_r |> println

    鈎*(-股^3 + 2*股^2*鈎 + 股^2*sqrt(股^2 + 鈎^2) + 股^2*sqrt(2*股^2 + 2*股*鈎 + 鈎^2) + 2*股*鈎^2 + 股*鈎*sqrt(股^2 + 鈎^2) - 股*鈎*sqrt(2*股^2 + 2*股*鈎 + 鈎^2) - 股*sqrt(2*股^4 + 2*股^3*鈎 + 3*股^2*鈎^2 + 2*股*鈎^3 + 鈎^4) + 鈎^3 - 鈎*sqrt(2*股^4 + 2*股^3*鈎 + 3*股^2*鈎^2 + 2*股*鈎^3 + 鈎^4))/(2*(3*股^3 + 5*股^2*鈎 - 2*股^2*sqrt(2*股^2 + 2*股*鈎 + 鈎^2) + 3*股*鈎^2 - 2*股*鈎*sqrt(2*股^2 + 2*股*鈎 + 鈎^2) + 鈎^3))

円の半径を求める式は長いが,具体的な鈎,股の値を代入すれば,円の直径が求まる。

2ans_r(鈎 => 3, 股 => 4).evalf() |> println

    0.547208709287844

鈎 = 3, 股 = 4 のとき,円の直径は 0.547208709287844 である。

プログラム的には,弦,s を一時変数として若干簡単に書くことができる。

function get_r(鈎, 股)
    弦 = sqrt(股^2 + 鈎^2)
    s = sqrt(2股^2 + 2股*鈎 + 鈎^2)
    r = 鈎*(-股^3 + 2股^2*鈎 + 股^2*弦 + 股^2*s + 2股*鈎^2 + 股*鈎*弦 - 股*鈎*s - 股*弦*s + 鈎^3 - 鈎*弦*s)/
        (2*(3股^3 + 5股^2*鈎 - 2股^2*s + 3股*鈎^2 - 2股*鈎*s + 鈎^3))
    println("円の直径 = $(2r)")
    return r
end
get_r(3, 4);

    円の直径 = 0.5472087092878442

術は 以下のように,簡約化されたものである。

function 術(鈎, 股)
    弦 = sqrt(鈎^2 + 股^2)
    天 = (sqrt((鈎 + 股)^2 + 股^2)*鈎 + (鈎 + 股)*弦+股^2)*(鈎 + 股)
    r = 股^2*鈎^2/天
    println("円の直径 = $(2r)")
end;
術(3, 4)

    円の直径 = 0.5472087092878439

function draw(鈎, 股, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    弦 = sqrt(股^2 + 鈎^2)
    s = sqrt(2股^2 + 2股*鈎 + 鈎^2)
    r = get_r(鈎, 股)
    (a, x) = (股*鈎/(股 + 鈎), (-r*股^2 - r*股*鈎 + r*股*s + r*鈎*s + 股*鈎^2)/(股^2 + 2股*鈎 + 鈎^2))
    plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
    rect(0, 0, a, a, :blue)
    segment(0, 鈎, a, 0, :magenta)
    circle(x, a + r, r)
    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, "股", :green, :left, :bottom, delta=delta)
        point(0, 鈎, "鈎", :green, :left, :bottom, delta=delta)
        point(a, 0, " a", :blue, :left, :bottom, delta=delta)
        point(x, a + r, "(x,a+r)", :red, :center, delta=-delta)
    end
end;

draw(3, 4, true)

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

算額(その1525)

2025年01月09日 | Julia

算額(その1525)

三十一 一関市舞川 観福寺内地蔵堂前額 明治34年(1901)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

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

キーワード:正五角形,折り紙
#Julia, #SymPy, #算額, #和算

内部に正五角形の空白(DGFHI)ができるように細長い長方形の紙を折る。紙の幅 AD が与えられたとき,長さを求める術を述べよ。

山村の図は不正確極まりない。実際に紙を折ってみたり,算額を解いてみたりしていないというのが透けて見える。

「今有如図」も細かいところは不正確である。折り紙をしたらどうなるかを理解していないから,直角であるところが直角でないし,線を延長すると通るべき所を通っていない。

正確な図は以下のようになる。赤が表,ベージュは裏

B, E の座標を (bx, by), (ex, ey)
B, E, C を円周上に持つ円の半径を R
D, G, F を円周上に持つ円の半径を r
とおき,BD の長さを求める。この時点では BD は R の関数である。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive
s18 = Sym(18)
s36 = Sym(36)
(ex, ey) = (R*cosd(s18), R*sind(s18))
(bx, by) = (R*sind(s36), -R*cosd(s36))
r = R*sind(s18)/cosd(s36)
BD = sqrt(bx^2 + (-r - by)^2) |> factor
BD |> println

    R*sqrt(160 - 32*sqrt(5))/(4*(1 + sqrt(5)))

紙の長さは 2(AB + BC) である。

AB = BD*sind(s18) |> simplify;
BC = sqrt((R - by)^2 + bx^2) |> simplify;
長さ = 2(AB + BC) |> simplify
長さ |> println

    R*(-sqrt(50 - 10*sqrt(5)) + 3*sqrt(10 - 2*sqrt(5)) + 4*sqrt(2*sqrt(5) + 10))/4

長さが R の関数として求められた。問は,長さを AD で表せということなので, R と AD の関係式を求める。

紙の幅 AD は BD*cosd(s18) である。

@syms AD
eq1 = AD - BD*cosd(s18)
eq1 |> println

    AD - R*sqrt(160 - 32*sqrt(5))*sqrt(sqrt(5)/8 + 5/8)/(4*(1 + sqrt(5)))

方程式 eq1 を解いて R が AD でどのように表されるかを見る。
R = AD*(sqrt(5) + 5)/5 である。

ans_R = solve(eq1, R)[1]
ans_R |> println

    AD*(sqrt(5) + 5)/5

長さを AD の関数で表す。

長さAD = 長さ(R => ans_R) |> sympy.sqrtdenest |> factor
長さAD |> println

    -AD*(-10*sqrt(2)*sqrt(sqrt(5) + 5) - 2*sqrt(10)*sqrt(sqrt(5) + 5) - 5*sqrt(2)*sqrt(5 - sqrt(5)) + sqrt(10)*sqrt(5 - sqrt(5)))/10

紙の幅 AD = 10 のときの紙の長さは 61.5536707435051 である。

長さAD(AD => 10).evalf() |> println

    61.5536707435051

「術」は 「紙の長さ = 2sqrt(sqrt(20) + 5) * 紙の幅」である。
2sqrt(sqrt(20) + 5) * 10 = 61.55367074350507 と,上の解と一致した。

「長さAD」は SymPy では簡約化できないが,「長さAD/AD」と術の「2sqrt(sqrt(20) + 5)」の差を取ると,ほぼゼロで,見かけの数式では同じように見えないが,両者は数値的に等しいことがわかる。

(長さAD/AD - 2sqrt(sqrt(Sym(20)) + 5)).evalf() |> println

    -0.e-123

function draw(AD, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    s18 = 18
    s36 = 36
    R = AD*(√5 + 5)/5
    r = R*sind(s18)/cosd(s36)
    (cx, cy) = (0, R)
    (dx, dy) = (0, -r)
    (ex, ey) = (R*cosd(s18), R*sind(s18))
    (bx, by) = (R*sind(s36), -R*cosd(s36))
    (ax, ay) = (AD*sind(s36), -AD*cosd(s36) - r)
    (fx, fy) = (r*sind(s36), r*cosd(s36))
    (gx, gy) = (r*cosd(s18), -r*sind(s18))
    r = R*sind(s18)/cosd(s36)
    plot(showaxis=false)
    plot!([bx, ex, cx, bx], [by, ey, cy, by], seriestype=:shape, color=:red, fillcolor=:red)
    plot!(-[bx, ex, cx, bx], [by, ey, cy, by], seriestype=:shape, color=:red, fillcolor=:red)
    plot!([fx, cx, -fx, fx], [fy, cy, fy, fy], seriestype=:shape, color=:wheat, fillcolor=:wheat)
    plot!([ax, bx, gx, dx, ax], [ay, by, gy, dy, ay], seriestype=:shape, color=:wheat, fillcolor=:wheat)
    plot!(-[ax, bx, gx, dx, ax], [ay, by, gy, dy, ay], seriestype=:shape, color=:wheat, fillcolor=:wheat)
    plot!([gx, ex, fx], [gy, ey, fy], linestyle=:dot, color=:wheat)
    plot!(-[gx, ex, fx], [gy, ey, fy], linestyle=:dot, color=:wheat)
    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(ax, ay, "A ", :green, :right, :vcenter)
        point(bx, by, " B", :green, :left, :vcenter)
        point(cx, cy, "C", :green, :center, :bottom, delta=delta)
        point(dx, dy, "D", :green, :center, :bottom, delta=delta)
        point(ex, ey, " E", :green, :left, :vcenter)
        point(fx, fy, "F ", :green, :right, delta=-delta/2)
        point(gx, gy, "G ", :green, :right, :bottom)
        point(-fx, fy, " H", :green, :left, delta=-delta/2)
        point(-gx, gy, " I", :green, :left, :bottom)
    end
end;

draw(10, true)

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

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

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