裏 RjpWiki

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

算額(その890)

2024年09月02日 | Julia

算額(その890)

改訂版:2024/09/02

七〇 加須市大字外野 棘脱地蔵堂 明治6年(1873)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.
キーワード:円7個,外円,円弧

外円内に円弧 6 個と等円 6 個を入れる。外円の直径が 3 寸のとき,等円の直径はいかほどか。

円弧の半径は外円の半径と同じ。また,等円の中心は正六角形を構成する。
外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, y); 第1象限に中心がある等円
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive, y::positive
y = R/2  # R*sin(PI/6)
eq1 = x^2 + y^2 - (R - r)^2  # 等円が外円に内接する
eq2 = R*(1 - cos(PI/6)) +  r - x  # 等円の中心の x 座標
res = solve([eq1, eq2], (r, x))[1];

等円の半径 r は 外円の半径 R の (3√3 - 1)/13 倍である。

res[1] |> factor |> println

    R*(-1 + 3*sqrt(3))/13

外円の直径が 3 寸のとき,等円の直径は 3(3√3 - 1)/13 = 0.9683428667784535 である。

3*(3√3 - 1)/13

    0.9683428667784535

「術」には,「置二十四個開平方一個減以除十三個乗外円径」とあるが,数式は「外円直径 × (√24 - 1)/13」であるが,「置二十四個」ではなく「置二十七個」の誤り(誤記?)で,「外円直径 × (√27 - 1)/13 = 外円直径 × (3√3 - 1)/13」となり,前述の式と一致する。
なお,「答」は,この修正のあとでも一致しない。

等円の中心の y 座標は y = R\*sin(π/6) = R/2 であるが,x 座標は x = R\*(24 - 7√3)/26 である。

res[2] |> factor |> println

  -R*(-24 + 7*sqrt(3))/26

function draw(R)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    y = R/2  # R*sind(30)
    r = R*(3√3 - 1)/13
    x = R*(24 - 7√3)/26
    @printf("外円の直径 = %g;  等円の直径 = %g\n", 2R, 2r)
    plot()
    delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
    dimension_line(0, 0, x, y, "R-r", :gray70, deltax=15delta, delta=5delta)
    point(0,0)
    hline!([0], color=:gray80, lw=0.5)
    vline!([0], color=:gray80, lw=0.5)
    x0 = R*cosd(30)
    y0 = R*sind(30)
    segment(-x0, y0, x0, y0, :gray70, linestyle=:dashdot)
    circle(0, 0, R, :green)
    rotate(x, y, r, angle=60)
    for θ = 30:60:330
        circle(R*cosd(θ), R*sind(θ), R, :blue, beginangle=180 + θ, endangle=240 + θ)
    end
    point(x, y, "(x,y)", :red, :center, :bottom, delta=delta)
    point(0, R, "R", :green, :center, :bottom, delta=delta)
end;

draw(3/2)

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

算額(その1267)

2024年09月02日 | Julia

算額(その1267)

百三十 現存するが奉納場所不明 明治14年(1881)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円3個,正三角形,正方形

正三角形と正方形が交わっており,その区分領域に大円,中円,小円を容れる。正方形の一辺の長さが与えられたとき小円の直径を求めよ。

正三角形の一辺の長さを a
正方形の一辺の長さを b
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (x3, r3)
とおき,以下の連立方程式を解く。

1. a, r2, x2 を求める

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     r1::positive, r2::positive, x2::positive;
# r1 = b/2
b = 2r1
eq1 = b/(a - b) - √Sym(3)
eq2 = (a - x2)/r2 - √Sym(3)
eq3 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand;

ans_a = solve(eq1, a)[1]
ans_a |> println

   2*r1*(sqrt(3) + 3)/3

eq12 = eq2(a => ans_a) |> simplify;
eq13 = eq3(a => ans_a);
res1 = solve([eq12, eq13], (r2, x2))[1]

   (r1*(-2*sqrt(sqrt(3) + 3) + sqrt(3) + 4)/3, r1*(-2*sqrt(3) + 3 + 2*sqrt(3*sqrt(3) + 9))/3)

2. r3, x3 を求める

前節で求めた a, r2, x2 を既知として r3, x3 を求める。

@syms a::positive, b::positive,
     r1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive;
b = 2r1
# r1 = b/2
a = 2r1*(√Sym(3) + 3)/3
r2 = r1*(-2sqrt(√Sym(3) + 3) + √Sym(3) + 4)/3
x2 = r1*(-2√Sym(3) + 3 + 2sqrt(3√Sym(3) + 9))/3
eq4 = (x3 - r1)^2 + (r1 - r3)^2 - (r1 + r3)^2 |> expand
eq5 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2 |> expand;
res2 = solve([eq4, eq5], r3, x3)[2];  # 2 of 2

@syms d
ans_r3 = apart(res2[1]) |> simplify
ans_r3 |> println

   r1*(-204*sqrt(sqrt(3) + 3) - 92*sqrt(3*sqrt(3) + 9) + 229*sqrt(3) + 449)/338

r1 = b/2 として,r3 は r1 の (-204*sqrt(sqrt(3) + 3) - 92*sqrt(3*sqrt(3) + 9) + 229*sqrt(3) + 449)/338 倍である。

たとえば,b = 10 のとき,r3 = 0.8171180745896389,直径は 1.6342361491792778 である。

ans_x3 = apart(res2[2]) |> simplify
ans_x3 |> println

   r1*(-9*sqrt(sqrt(3) + 3) - sqrt(3*sqrt(3) + 9) + 8*sqrt(3) + 33)/13

b = 10 のとき,各パラメータは以下のとおりである。

   b = 10;  a = 15.7735;  r1 = 5;  r2 = 2.30233;  x2 = 11.7858;  r3 = 0.817118;  x3 = 9.04257

function draw(b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = b/2
   a  = 2r1*(√3 + 3)/3
   r2 = r1*(√3 + 4 - 2sqrt(√3 + 3))/3
   x2 = r1*(2sqrt(3√3 + 9) -2√3 + 3)/3
   r3 = r1*(229√3 + 449 - 204sqrt(√3 + 3) - 92sqrt(3√3 + 9))/338
   x3 = r1*(8√3 + 33 - 9sqrt(√3 + 3) - sqrt(3√3 + 9))/13
   @printf("b = %g;  a = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", b, a, r1, r2, x2, r3, x3)
   plot([0, a, a/2, 0], [0, 0, √3a/2, 0], color=:blue, lw=0.5)
   plot!([0, b, b, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5) 
   circle(r1, r1, r1)
   circle(x2, r2, r2, :orange)
   circle(x3, r3, r3, :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, 0, "a", :blue, :center, delta=-delta)
       point(0, b, "b ", :green, :right, :vcenter)
       point(b, 0, "b", :green, :left, delta=-delta)
       point(a/2, √3a/2, "(a/2,√3a/2)", :blue, :left, :bottom, delta=delta/2)
       point(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta)
       point(x2, r2, "中円:r2\n(x2,r2)", :orange, :center, delta=-delta)
       point(x3, r3, "小円:r3,(x3,r3)", :magenta, :right, :bottom, delta=2delta, deltax=-4delta)
       plot!(xlims=(-5delta, a + 4delta), ylims=(-5delta, √3a/2 + 5delta))
   end
end;

draw(10, true)

 

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

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

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