裏 RjpWiki

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

算額(その97)

2023年01月11日 | Julia

算額(その97)

東京都台東区下谷池端 東淵寺 明治2年(1869)
http://www.wasan.jp/tokyo/tobuti.html

正方形の中に半円,大円,小円が図のように配置されている。大円の径が二寸四分八厘のとき,小円の径はいくつか。

大円の半径を 248 として,図のように記号を定めて方程式を解く。

using SymPy
@vars R, r, x, y (positive=true, real=true);
@syms R::positive, r::positive, x::positive, y::positive;
R =248
eq1 = x^2 + (y - r)^2 - 4r^2
eq2 = x^2 + y^2 - R^2
eq3 = x^2 + (R + 2r - y)^2 - (R + r)^2;

(r, x, y) = solve([eq1, eq2, eq3], (r, x, y))[1]

   (-31*sqrt(17)*(109 + 27*sqrt(17))^(2/3)/12 - 961*(109 + 27*sqrt(17))^(1/3)/27 - 1240/27 + 31*sqrt(17)*(109 + 27*sqrt(17))^(1/3)/3 + 1457*(109 + 27*sqrt(17))^(2/3)/108, 248*sqrt(6)*sqrt(-(109 + 27*sqrt(17))^(5/3)*(909*sqrt(17) + 3755) - 2*(109 + 27*sqrt(17))^(4/3)*(4627 + 1125*sqrt(17)) + 32*(109 + 27*sqrt(17))*(2943*sqrt(17) + 12137))/(9*(109 + 27*sqrt(17))^(3/2)), -1984/(81*(109/729 + sqrt(17)/27)^(1/3)) + 248/9 + 248*(109/729 + sqrt(17)/27)^(1/3))

小円の半径は,100.00743003769321 である。元の単位で径に直せばほぼ 1 寸ということで,算額の答えにあっている。

-31*sqrt(17)*(109 + 27*sqrt(17))^(2/3)/12 - 961*(109 + 27*sqrt(17))^(1/3)/27 - 1240/27 + 31*sqrt(17)*(109 + 27*sqrt(17))^(1/3)/3 + 1457*(109 + 27*sqrt(17))^(2/3)/108

   100.00743003769321


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; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, 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 draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="")#, fontfamily="IPAMincho")
    (r, x, y) = (-31*sqrt(17)*(109 + 27*sqrt(17))^(2/3)/12 - 961*(109 + 27*sqrt(17))^(1/3)/27 - 1240/27 + 31*sqrt(17)*(109 + 27*sqrt(17))^(1/3)/3 + 1457*(109 + 27*sqrt(17))^(2/3)/108, 248*sqrt(6)*sqrt(-(109 + 27*sqrt(17))^(5/3)*(909*sqrt(17) + 3755) - 2*(109 + 27*sqrt(17))^(4/3)*(4627 + 1125*sqrt(17)) + 32*(109 + 27*sqrt(17))*(2943*sqrt(17) + 12137))/(9*(109 + 27*sqrt(17))^(3/2)), -1984/(81*(109/729 + sqrt(17)/27)^(1/3)) + 248/9 + 248*(109/729 + sqrt(17)/27)^(1/3))
    R = 248
    println("r = $r;  x = $x;  y = $y")
    plot()
    circle(0,2r + R, R, :green)
    circle(0, r, r)
    circle(x, y, r)
    circle(-x, y, r)
    circle(0, 0, R + r, :blue, beginangle=0, endangle=180)
    segment(-R-r, 0, R+r, 0, :blue)
    l = r + R
    plot!([-l, l, l, -l, -l], [0, 0, 2l, 2l, 0], color=:black, lw=0.5)
    if more
        point(0, 2r + R, "2r+R ", :green, :right)
        point(0, r, "r ", :red, :right)
        point(x, y, "(x,y)", :red, :top)
        hline!([0], color=:black, lw=0.5)
        vline!([0], color=:black, lw=0.5)
    end
end;

   r = 100.00743003769321;  x = 191.57807116330426;  y = 157.48600778909824

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

算額(その96)

2023年01月11日 | Julia

算額(その96)

東京都台東区下谷池端 東淵寺 明治2年(1869)
http://www.wasan.jp/tokyo/tobuti.html

大円2個,地円3個,人円3個が図のように配置されている。大円の径が 34 寸のとき,人円の径はいくつか。

大円の半径を 34 とし,図のように地円,人円の半径を r2, r3,人円の中心座標を (x, y) として方程式を解く。

using SymPy

@syms r1::positive, r2::positive, r3::positive, x::positive, y::positive;

r1 = 34
eq1 = x^2+ (y-2r2-r3)^2 - (r1+r3)^2 |> expand;
eq2 = (2r2 - x)^2 + (y - r2)^2 - (r1 + r2)^2 |> expand;
eq3 = r2^2 + (r2 - r3)^2 - (r2 + r3)^2 |> expand;
eq4 = x^2 + (y - r2)^2 - (r1 + r2)^2 |> expand;
a = solve([eq1, eq2, eq3, eq4], (r2, r3, x, y))

   1-element Vector{NTuple{4, Sym}}:
    (-68*sqrt(10)*sqrt(sqrt(5) + 3)/5 + 136*sqrt(5)/5 + 408/5, -17*sqrt(10)*sqrt(sqrt(5) + 3)/5 + 34*sqrt(5)/5 + 102/5, -68*sqrt(10*sqrt(5) + 30)/5 + 136*sqrt(5)/5 + 408/5, 238/5 + 136*sqrt(5)/5)

人円の半径は 11.002631123499281 である。

-17*sqrt(10)*sqrt(sqrt(5) + 3)/5 + 34*sqrt(5)/5 + 102/5

   11.002631123499281

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, vertical, position, color))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 34
   (r2, r3, x, y) = (-68*sqrt(10)*sqrt(sqrt(5) + 3)/5 + 136*sqrt(5)/5 + 408/5, -17*sqrt(10)*sqrt(sqrt(5) + 3)/5 + 34*sqrt(5)/5 + 102/5, -68*sqrt(10*sqrt(5) + 30)/5 + 136*sqrt(5)/5 + 408/5, 238/5 + 136*sqrt(5)/5)
   println("地 = $r2, 人 = $r3, x = $x, y= $y")
   plot()
   circle(x, y, r1, :red)
   circle(-x, y, r1, :red)
   circle(0, 2r2 + r3, r3, :blue)
   circle(r2, r3, r3, :blue)
   circle(-r2, r3, r3, :blue)
   circle(0, r2, r2, :magenta)
   circle(2r2, r2, r2, :magenta)
   circle(-2r2, r2, r2, :magenta)
   hline!([0], color=:black, lw=0.5)
   if more
       point(0, r2, "r2 ", :magenta, :right)
       point(x, y, "(x,y)", :red, :top)
       point(r2, r3, "(r2,r3)", :blue, :top)
       vline!([0], color=:black, lw=0.5)
   end
end;

   地 = 44.010524493997124, 人 = 11.002631123499281, x = 44.010524493997124, y= 108.42104898799428

 

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

算額(その95)

2023年01月11日 | Julia

算額(その95)

東京都台東区下谷池端 東淵寺 明治2年(1869)
http://www.wasan.jp/tokyo/tobuti.html

半円の中に等円 2 個と半円 1 個が入っている。外円の径を 2 寸 7 分 とすると等円の径はいくつか。

図のように記号を定め方程式を解く。

using SymPy

@syms R::positive, r::positive, x::positive, y::positive;

R = 27
eq1 = x^2 + (y - r)^2 - 4r^2 |> expand;
eq2 = x^2 + r^2 - (R - r)^2 |> expand;
eq3 = r^2 + y^2 - R^2 |> expand;
a = solve([eq1, eq2, eq3], (r, x, y));

方程式自体は簡単なものなのであるが,方程式の解 a をそのまま書き出すと,虚数単位 I が含まれたりする恐ろしい数式になる。それを .evalf() で数値になおすと以下の3組になる。
まだ I が残っているが 4.74210720465186e-23*I などで誤差の範囲で実数である。2 番目と 3 番目の解は不適切解である。
つまるところ r = 10.0929641600229, x = 13.5639203535985, y = 25.0426051852537 ということであろう。
よって等円の径は 1 寸 0 分 0 厘 9 毛である。

for i in 1:3
   println("r = $(a[i][1].evalf());  x = $(a[i][2].evalf());  y = $(a[i][3].evalf())")
end

   r = 10.0929641600229;  x = 13.5639203535985 + 4.74210720465186e-23*I;  y = 25.0426051852537 + 0.e-22*I
   r = -22.2346710141611 - 4.33680868994202e-19*I;  x = 43.928034724589 + 1.17139889977585e-22*I;  y = -15.3172910428713 - 0.e-21*I
   r = 17.5417068541382 + 4.33680868994202e-19*I;  x = 14.7733601500627*I;  y = -20.5253141423824 + 0.e-21*I

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, vertical, position, color))
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 draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r, x, y) = (27, 10.0929641600229, 13.5639203535985, 25.0426051852537)
   plot()
   circle(0, 0, R, :red, beginangle=0, endangle=180)
   circle(x, r, r, :blue)
   circle(-x, r, r, :blue)
   circle(0, y, r, :blue, beginangle=180, endangle=360)
   segment(-r, y, r, y, :blue)
   segment(-R, 0, R, 0, :red)
   if more
       point(0, y, "y ", :blue, :right)
       point(0, r, "r ", :blue, :right)
       point(x, r, "(x,r)", :blue, :top)
       point(x, 0, "  x", :blue, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

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

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

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