裏 RjpWiki

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

算額(その219)

2023年05月01日 | Julia

算額(その219)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2(180) 長野県 渋薬師堂 上水内牟礼村牟礼大久保善賢氏保管 嘉永2年(1849)

外円に内接する 3 個の大円および小円,大円が交わる領域に 4 個の中円が入っている。外円の径が 14 寸のとき,大円,中円,小円の径を求めよ。

算額(その152) 岐阜県大垣市西外側町 大垣八幡神社 天保年間 を一段階複雑にしたもの

図のように記号を定め,連立方程式を解く。
外円の半径を r0 とする
右上にある小円の中心座標,半径を (x1, y1), r1
右上にある中円の中心座標,半径を (x2, y2), r2
右下にある大円の中心座標,半径を (x3, y3), r3

using SymPy

@syms r0::positive, r1::positive, r2::positive, r3::positive

r0 = 14
x1 = (r0 - r1) * cos(PI/6)
y1 = (r0 - r1) * sin(PI/6)
x2 = 2r2 * cos(PI/6)
y2 = 2r2 * sin(PI/6)
x3 = (r0 - r3) * cos(PI/6)
y3 = (r3 - r0) * sin(PI/6)
eq1 = 2r3 - r0 - r2  # 半径間の関係
eq3 = x2^2 + (r0 - r3 - y2)^2 - (r3 - r2)^2;  # 中円が大円に内接
eq5 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2;  # 小円と大円が外接

r0 を与えて solve() すると,数値解が得られるが,r0 との関係がわかりにくいので,r0 を未知数のまま解くと,それぞれの解に r0 が含まれる。

res = solve([eq1, eq3, eq5], (r3, r2, r1))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (3*r0/5, r0/5, r0/7)

つまり,大円,中円,小円の半径(直径)は外円の半径(直径)の 3/5, 1/5, 1/7 である。

   外円の半径 = 14.000,  大円の半径 = 8.400,  中円の半径 = 2.800;  小円の半径 = 2.000

外円の径が 14 寸のとき,大円,中円,小円の径は 8寸4分,2寸8分,2寸 である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
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 draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 14
   (r1, r2, r3) = (r0/7, r0/5, 3*r0/5)
   x1 = (r0 - r1) * cos(PI/6)
   y1 = (r0 - r1) * sin(PI/6)
   x2 = 2r2 * cos(PI/6)
   y2 = 2r2 * sin(PI/6)
   x3 = (r0 - r3) * cos(PI/6)
   y3 = (r3 - r0) * sin(PI/6)
   @printf("外円の半径 = %.3f,  大円の半径 = %.3f,  中円の半径 = %.3f;  小円の半径 = %.3f\n", r0, r3, r2, r1)
   plot()
   circle(0, 0, r0, :blue)
   circle(x1, y1, r1, :magenta)
   circle(-x1, y1, r1, :magenta)
   circle(0, r1 - r0, r1, :magenta)
   circle(0, r0 - r3, r3, :green)
   circle(x3, y3, r3, :green)
   circle(-x3, y3, r3, :green)
   circle(x2, y2, r2)
   circle(-x2, y2, r2)
   circle(0, -2r2, r2)
   circle(0, 0, r2)

   if more == true
       point(x3, y3, " 大円(x3,y3,r3)")
       point(x2, y2, " 中円(x2,y2,r2)", :red)
       point(x1, y1, " 小円(x1,y1,r1)", :magenta)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その218)

2023年05月01日 | Julia

算額(その218)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2(163)
長野県小諸市大久保 釈尊寺布引観音堂 天保14年(1843)

4 つの弧が交わり,中央部に大円と他に 8 個の等円を置く。大円の径が 1 寸のとき,等円の径はいかほどか。

大円,等円の半径をそれぞれ r1, r2 とおく。弧は半径 r0 の円の一部とする。
原点を中心として点対称なので,第1象限のみを考える。左下の等円の中心を (x2, x2) とする。右上の等円の中心を (x3, x3) とするが,x3 = (r1 + r0)/2 である。
以下の方程式を解く。

using SymPy

@syms r0::positive, r1::positive, r2::positive, x2::positive, x3::positive, x::positive, y::positive;

r1 = 1
eq1 = (r1 + r0 - x3)^2 + x3^2 - (r0 - r2)^2 |> expand |> simplify
eq2 = (r1 + r0 - x2)^2 + x2^2 - (r0 + r2)^2 |> expand |> simplify
eq3 = 2x2^2 - (r1 + r2)^2 |> expand |> simplify
eq4 = (r1 + r0) / 2 - x3
res = solve([eq1, eq2, eq3, eq4], (r0, r2, x2, x3))

   2-element Vector{NTuple{4, Sym}}:
    (-sqrt(2) + sqrt(14 - 8*sqrt(2)) + 3, -3*sqrt(2) - sqrt(7 - 4*sqrt(2)) + sqrt(14 - 8*sqrt(2)) + 4, -3 - sqrt(14 - 8*sqrt(2))/2 + sqrt(7 - 4*sqrt(2)) + 5*sqrt(2)/2, -sqrt(2)/2 + sqrt(2)*sqrt(7 - 4*sqrt(2))/2 + 2)
    (-sqrt(40*sqrt(6) + 98) + 7 + 3*sqrt(6), -sqrt(60*sqrt(6) + 147) + 2 + sqrt(6) + sqrt(40*sqrt(6) + 98), -sqrt(60*sqrt(6) + 147) - sqrt(6)/2 - 1 + 3*sqrt(40*sqrt(6) + 98)/2, -sqrt(10*sqrt(6) + 49/2) + 3*sqrt(6)/2 + 4)

最初の組が適切解である。

   r0 = 3.22478;  r1 = 1.00000;  r2 = 0.23741;  x2 = 0.87498;  x3 = 2.11239
   r2 = 0.23740866273917272

r2 = -3*sqrt(2) - sqrt(7 - 4*sqrt(2)) + sqrt(14 - 8*sqrt(2)) + 4
   = 0.23740866273917272

すなわち,0寸2分3厘7毛408663 である。

using Plots
using Printf

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 circle4(x, y, r, color=:red)
  circle(x, y, r, color)
  circle(x, -y, r, color)
  circle(-x, y, r, color)
  circle(-x, -y, r, color)
end;

function point(x, y, string="", color=:blue, 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 draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1
   (r0, r2, x2, x3) = (-sqrt(2) + sqrt(14 - 8*sqrt(2)) + 3,
       -3*sqrt(2) - sqrt(7 - 4*sqrt(2)) + sqrt(14 - 8*sqrt(2)) + 4,
       -3 - sqrt(14 - 8*sqrt(2))/2 + sqrt(7 - 4*sqrt(2)) + 5*sqrt(2)/2,
       -sqrt(2)/2 + sqrt(2)*sqrt(7 - 4*sqrt(2))/2 + 2)
   @printf("r0 = %.5f;  r1 = %.5f;  r2 = %.5f;  x2 = %.5f;  x3 = %.5f\n", r0, r1, r2, x2, x3)
   println("r2 = $r2")
   lim = r1 + r0 + 0.1
   plot()
   circle(r0 + r1, 0, r0, :black)
   circle(-r0 - r1, 0, r0, :black)
   circle(0, r0 + r1, r0, :black)
   circle(0, -r0 - r1, r0, :black)
   circle(0, 0, r1)
   circle4(x2, x2, r2, :green)
   circle4(x3, x3, r2, :green)
   plot!(xlims=(-lim, lim), ylims=(-lim, lim))
   if more == true
       point(0, r1, " r1", :blue, :right, :bottom)
       point(r1, 0, "r1 ", :blue, :right)
       point(0, r1 + r0, " r1+r0")
       point(r1 + r0, 0, "r1+r0 ", :blue, :right)
       point(x3, x3, "   (x3,x3)")
       point(x2, x2, "   (x2,x2)")
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

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

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

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