裏 RjpWiki

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

Julia 1.8.4

2022年12月25日 | ブログラミング

Julia 1.8.4 が出ました

Current stable release: v1.8.4 (December 23, 2022)

 

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

算額(その74)

2022年12月25日 | Julia

算額(その74)

大阪府茨木市 井於神社 弘化3年(1846)
http://www.wasan.jp/osaka/iyo.html

鉤3寸,股4寸,弦5寸に甲円,乙円が入っている。それぞれの径を求めよ。

径は SymPy を使うまでもなく筆算(暗算)でも答えを出せる。

using SymPy
@syms AD::positive, DC::positive, BD::positive;
(AB, BC, CA) = (3, 4, 5)
eq1 = BD - AB * 4//5;
eq2 = AD - AB * 3//5;
eq3 = DC - (CA - AD);
result = solve([eq1, eq2, eq3], (BD, AD, DC))

   Dict{Any, Any} with 3 entries:
     BD => 12/5
     AD => 9/5
     DC => 16/5

r(鉤, 股, 弦) = 鉤 + 股 - 弦;  # 鉤股弦に内接する円の径(直径)

r(result[BD], result[DC], BC).evalf() |> println  # 甲円: 1.6 = 1寸6分

   1.60000000000000

r(result[AD], result[BD], AB).evalf() |> println  # 乙円: 1.2 = 1寸2分

   1.20000000000000

円の中心座標を求めて図を描く。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; fontsize=10, mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, fontsize, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))
end;

function centerofcircle(x1, y1, x2, y2, x3, y3)
   """ 鉤股弦に内接する円の中心座標を求める"""
   a = (x1 - x2)^2 + (y1 - y2)^2
   b = (x1 - x3)^2 + (y1 - y3)^2
   c = (x3 - x2)^2 + (y3 - y2)^2
   if isapprox(a + b, c) || isapprox(b + c, a) || isapprox(c + a, b)
       @syms Ox::positive, Oy::positive
       (a, b, c) = sqrt.((a, b, c))
       r = (a + b + c)/2 - max(a, b, c)
       eq1 = distance(x1, y1, x2, y2, Ox, Oy) - r
       eq2 = distance(x1, y1, x3, y3, Ox, Oy) - r
       res = solve([eq1, eq2])
       if typeof(res) == Dict{Any, Any}
           return (res[Ox], res[Oy], r)
       else
           for i in 1:length(res)
               x, y = res[i][Ox], res[i][Oy]
               if min(x1, x2, x3) <= x <= max(x1, x2, x3) &&
                   min(y1, y2, y3) <= y <= max(y1, y2, y3)
                   return (x, y, r)
               end
           end
       end
   else
       println("鉤股弦ではない")
       return
   end
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (AB, BC, CA) = (3, 4, 5)
   (BD, AD, DC) = (12/5, 9/5, 16/5)
   (Ax, Ay) = (0, AB)
   (Bx, By) = (0, 0)
   (Cx, Cy) = (BC, 0)
   gap = 0.15
   (sine, cosine) = (3, 4) ./ 5
    plot([Ax, Bx, Cx, Ax], [Ay, By, Cy, Ay], color=:black, lw=0.5,
        xlims=(-gap, Cx+gap), ylims=(-1.5gap, Ay+gap))
   (Dx, Dy) = (AD*cosine, Ay-AD*sine)
   segment(Bx, By, Dx, Dy, :black)
   (r1, r2) = (1.6, 1.2) ./ 2
   println("甲円の半径 = $r1;  乙円の半径 = $r2")
   O1x, O1y, r = centerofcircle(Bx, By, Dx, Dy, Cx, Cy)
   circle(O1x, O1y, r)
   O2x, O2y, r = centerofcircle(Ax, Ay, Bx, By, Dx, Dy)
   circle(O2x, O2y, r)
    if more
       point(Ax, Ay, "A ", :black, :right)
       point(Bx, By, "B ", :black, :right)
       point(Cx, Cy, " C", :black)
       point(Dx, Dy, "  D", :black)
       point(O2x, O2y, "乙円\n(O2x,O2y)", :red, :center)
       point(O1x, O1y, "甲円\n(O1x,O1y)", :red, :center)
    end
end;

   甲円の半径 = 0.8;  乙円の半径 = 0.6

 

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

算額(その73)

2022年12月25日 | Julia

算額(その73)

三九 加須市不動岡 総願寺不動堂 弘化5年(1848)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

埼玉県加須市不動岡 総願寺不動堂 弘化5年
http://www.wasan.jp/saitama/hudou.html

キーワード:円3個,直角三角形,正方形

鉤股弦の中に,1 個の正方形と大きさの異なる 3 個の円が互いに接して入っている。円の径を求めよ。

座標の割当の都合で,算額を左右反転させて記号を割り振り,方程式を解く。

円の径だけならば「径 = 鉤 + 股 - 弦」でよいが,図を描くためには円の中心座標を求めなければならない。

また,鉤股弦のいずれかが定数として与えられているのだろうが,写真では不鮮明なので,任意の定数を設定する。

鉤:AB, 股:BC, 弦:CA を 3, 4, 5 として各線分の長さを求めると 「長さ/37」 になるので,鉤股弦を 37 倍しておけば整数解が得られる。

using SymPy
@syms AD::positive, DE::positive, EA::positive,
     EB::positive, BF::positive,
     GC::positive, CF::positive;
(AB, BC, CA) = (3, 4, 5) .* 37
EF = DE
FG = DE
eq1 = AD - DE*(AB//BC);
eq2 = EB - DE*(AB//CA);
eq3 = FG - GC*(AB//BC);
eq4 = EA + EB - AB;
eq5 = BF + CF - BC;
eq6 = AD + EF + GC - CA;
eq7 = FG - GC * EB / BF
results = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7],
   (AD, EA, EB, BF, CF, DE, GC))

    (45, 75, 36, 48, 100, 60, 80)

小円,中円,大円の径は 12, 15, 20 である。

(AD, EA, EB, BF, CF, DE, GC) = (45, 75, 36, 48, 100, 60, 80)
EF = DE
FG = DE
(EB + BF - EF)/2 |> println
(AD + DE - EA)/2 |> println
(FG + GC - CF)/2 |> println

   12.0
   15.0
   20.0

function centerofcircle(x1, y1, x2, y2, x3, y3)
    """ 鉤股弦に内接する円の中心座標を求める"""
    a = (x1 - x2)^2 + (y1 - y2)^2
    b = (x1 - x3)^2 + (y1 - y3)^2
    c = (x3 - x2)^2 + (y3 - y2)^2
    if isapprox(a + b, c) || isapprox(b + c, a) || isapprox(c + a, b)
        @syms Ox::positive, Oy::positive
        (a, b, c) = sqrt.((a, b, c))
        r = (a + b + c)/2 - max(a, b, c)
        eq1 = dist2(x1, y1, x2, y2, Ox, Oy, r)
        eq2 = dist2(x1, y1, x3, y3, Ox, Oy, r)
        res = solve([eq1, eq2])
        for i = 1:length(res)
            (Ox1, Oy1) = (res[i][Ox], res[i][Oy])
            ((minimum([x1, x2, x3]) <= Ox1 <= maximum([x1, x2, x3])) &&
            (minimum([y1, y2, y3]) <= Oy1 <= maximum([y1, y2, y3]))) && return (float(Ox1), float(Oy1), r)
        end
    else
        println("鉤股弦ではない")
        return
    end
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    gap = 10
    (sine, cosine) = (3, 4) ./ 5
    (AD, EA, EB, BF, CF, DE, GC)= (45, 75, 36, 48, 100, 60, 80)
    (AB, BC, CA) = (3, 4, 5) .* 37
    EF = DE
    FG = DE
    (Ax, Ay) = (0, AB)
    (Bx, By) = (0, 0)
    (Cx, Cy) = (BC, 0)
    plot([Ax, Bx, Cx, Ax], [Ay, By, Cy, Ay], color=:black, lw=0.5,
         xlims=(-gap, Cx+gap), ylims=(-1.5gap, Ay+gap))
    (Dx, Dy) = (AD*cosine, Ay-AD*sine)
    (Ex, Ey) = (0, EB)
    (Fx, Fy) = (BF, 0)
    segment(Dx, Dy, Ex, Ey, :black)
    segment(Ex, Ey, Fx, Fy, :black)
    (Gx, Gy) = (Dx + BF, Dy - EB)
    segment(Fx, Fy, Gx, Gy)
    r3 = (EB + BF - EF)/2
    r2 = (AD + DE - EA)/2
    r1 = (FG + GC - CF)/2
    println("小円の半径 = $r3;  中円の半径 = $r2;  大円の半径 = $r1")
    O3x, O3y, r = centerofcircle(Bx, By, Ex, Ey, Fx, Fy)
    circle(O3x, O3y, r)
    O2x, O2y, r = centerofcircle(Ax, Ay, Ex, Ey, Dx, Dy)
    circle(O2x, O2y, r)
    O1x, O1y, r = centerofcircle(Fx, Fy, Cx, Cy, Gx, Gy)
    println((Fx, Fy, Cx, Cy, Gx, Gy))
    circle(O1x, O1y, r)
    if more
        point(Ax, Ay, "A ", :black, :right)
        point(Bx, By, "B ", :black, :right)
        point(Cx, Cy, " C", :black)
        point(Dx, Dy, "  D", :black)
        point(Ex, Ey, " E", :black)
        point(Fx, Fy, "\nF", :black)
        point(Gx, Gy, "  G", :black)
        point(O3x, O3y, "小円\n(O3x,O3y)", :red, :center)
        point(O2x, O2y, "中円\n(O2x,O2y)", :red, :center)
        point(O1x, O1y, "大円\n(O1x,O1y)", :red, :center)
    end
end;

   小円の半径 = 12.0;  中円の半径 = 15.0;  大円の半径 = 20.0

 

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

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

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