裏 RjpWiki

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

算額(その230)

2023年05月13日 | Julia

算額(その230)

(12) 京都市伏見区御香宮門前町 御香宮神社(ごこうのみやじんじゃ)文久3年(1863)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4516 初版第一刷,大阪教育図書株式会社,大阪市.

絵馬堂が意外におもしろい!(3) - 御香宮神社 -
http://flattravel.blog.fc2.com/blog-entry-16.html

キーワード:円3個,長方形

大円の一部を欠き取る長方形内に,天円と地円を入れる。長方形の長辺と短辺がそれぞれ 4寸,3寸で,大円の径が 13寸のとき,天円と地円の径はいかほどか。

長方形の左下角の座標を (x, y) とおく。
大円の半径,中心座標を r0, (0, 0)
天円の半径,中心座標を r1, (x1, y1) # (x + r1, y + 3 - r1)
地円の半径,中心座標を r2, (x2, y2) # (x + 4 - r2, y + r2)
とおき,以下の方程式を解く。
注: x, x1, x2 は ::negative である。

using SymPy

@syms r0::positive, x::negative, y::positive,
   r1::positive, x1::negative, y1::positive,
   r2::positive, x2::negative, y2::positive;
r0 = 13//2
(x1, y1) = (x+r1, y + 3 - r1)
(x2, y2) = (x + 4 - r2, y + r2)
eq1 = x1^2 + y1^2 - (r0 + r1)^2
eq2 = x2^2 + y2^2 - (r0 - r2)^2
eq3 = x^2 + y^2 - r0^2
eq4 = (x + 4)^2 + (y + 3)^2 - r0^2

res = solve([eq1, eq2, eq3, eq4], (r1, r2, x, y))

   2-element Vector{NTuple{4, Sym}}:
    (4/5, 6/5, -28/5, 33/10)
    (36, 6/5, -28/5, 33/10)

最初の解が妥当なものである。すなわち天円の径は 8/5 = 1.6 つまり 1寸6分,地円の径は 12/5 = 2.4 つまり 2寸4分である。

   r1 = 0.8000000;  r2 = 1.2000000;  x = -5.6000000;  y = 3.3000000
   天円径 = 1.6000000,  地円径 = 2.4000000

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 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")
   (r1, r2, x, y) = res[1]
   r0 = 13//2
   (x1, y1) = (x+r1, y + 3 - r1)
   (x2, y2) = (x + 4 - r2, y + r2)
   @printf("r1 = %.7f;  r2 = %.7f;  x = %.7f;  y = %.7f\n", r1, r2, x, y)
   @printf("天円径 = %.7f,  地円径 = %.7f\n", 2r1, 2r2)
   plot([x, x+4, x+4, x, x], [y, y, y+3, y+3, y], color=:blue, lw=0.5)
   circle(0, 0, r0)
   circle(x1, y1, r1, :green)
   circle(x2, y2, r2, :gold)
   if more == true
       point(x1, y1, "(x1,y1)", :green, :center, :top)
       point(x1, y1, "天円", :green, :center, :bottom)
       point(x2, y2, "(x2,y2)", :gold, :center, :top)
       point(x2, y2, "地円", :gold, :center, :bottom)
       point(x, y, "(x,y)", :blue)
       point(x+4, y+3, "(x+4,y+3)", :blue, :right, :bottom)
       point(0, r0, " r0", :red)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false, xlims=(-6, -0.5), ylims=(3.1, 6.35))
   end
end;

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

算額(その229)

2023年05月13日 | Julia

算額(その229)

京都府京都市上京区 北野天満宮 明治12年(1879)
http://www.wasan.jp/kyoto/kitano2.html

北野天満宮の算額について,永井信一
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/kitano.pdf

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

二等辺三角形(2つの等辺(上斜)の長さは 1014寸,底辺(下斜)の長さは 1428寸)の中に甲円が 2 個,乙円,丙円が 1 個ずつ入っている。それぞれの径を求めよ。

算額(その110)https://blog.goo.ne.jp/r-de-r/e/f9651e4698b7baef75c28a0647b16424 の愛媛県松山市の伊佐爾波神社の算額では「三角」とあり正三角形であるが,ここでは「圭」なので,二等辺三角形である。

算額(その2031)は甲円の直径のみ求めよとなっている。

甲円の半径,中心座標を r1, (x1, r1)
乙円の半径,中心座標を r2, (0, r3 + r2)
丙円の半径,中心座標を r3, (0, r3)
とおき,以下の方程式を解く。
なお,ピタゴラスの定理から,二等辺三角形の高さは sqrt(1014^2 - 714^2) = 720 である。

using SymPy

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))^2  # 二乗距離を返す!!
end;

@syms r1::positive, r2::positive, r3::positive, x1::positive;
eq1 = x1^2 + (2r3 + r2 - r1)^2 - (r1 + r2)^2
eq2 = x1^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = distance(0, 720, 714, 0, 0, 2r3 + r2) - r2^2  # 上斜への二乗距離
eq4 = distance(0, 720, 714, 0, x1, r1) - r1^2  # 上斜への二乗距離

res = solve([eq1, eq2, eq3, eq4], (r1, r2, r3, x1))

   3-element Vector{NTuple{4, Sym}}:
    (357/2, 15232/75, 2856/25, 1428/5)
    (1447992/1315 + 205632*sqrt(61)/1315, 723996/1315 - 84966*sqrt(61)/1315, -402696/1315 + 102816*sqrt(61)/1315, 17136*sqrt(61)/263 + 308448/263)
    (111384/55 + 34272*sqrt(14)/55, -414596/275 + 113288*sqrt(14)/275, 476/11 + 952*sqrt(14)/11, 2856*sqrt(14)/11 + 17136/11)

2 番めの解が妥当なものである。すなわち甲円,乙円,丙円の直径はそれぞれ 357 寸, 30464/75 寸, 5712/25 寸 である。

   r1 = 178.5000000;  r2 = 203.0933333;  r3 = 114.2400000;  x1 = 285.6000000
   甲円径 = 357.0000000,  乙円径 = 406.1866667;  丙円径 = 228.4800000

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 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")
   (r1, r2, r3, x1) = res[1]
   @printf("r1 = %.7f,  r2 = %.7f;  r3 = %.7f;  x1 = %.7f\n", r1, r2, r3, x1)
   @printf("甲円径 = %.7f,  乙円径 = %.7f;  丙円径 = %.7f\n", 2r1, 2r2, 2r3)
   plot([714, 0, -714, 714], [0, 720, 0, 0], color=:blue, lw=0.5)
   circle(0, 2r3 + r2, r2, :green)
   circle(0, r3, r3, :gold)
   circle(x1, r1, r1)
   circle(-x1, r1, r1)
   if more == true
       point(0, 720, " 720", :blue, :left, :bottom)
       point(714, 0, "714", :blue, :left, :bottom)
       point(0, 2r3 + r2, " 2r3+r2", :red)
       point(0, r3, " r3", :red)
       point(x1, r1, "(x1,r1)", :red)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その228)

2023年05月13日 | Julia

算額(その228)

冲方丁の「天地明察」中で,金王八幡神社に奉納されたという算額。実際には渋川春海が算額を奉納したという記録はないそうだ。

鉤が 9寸,股が 12寸の鉤股弦に等円を2個図のように入れる。円の径はいくつか。

等円の半径を r とする。
左上の等円の中心座標を (r, y)
右下の等円の中心座標を (x, r)
とおき,以下の方程式を解く。

using SymPy

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))^2  # 二乗距離を返す!!
end;

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

eq1 = (x - r)^2 + (y - r)^2 - 4r^2
eq2 = distance(0, 9, 12, 0, r, y) - r^2  # 弦への二乗距離
eq3 = distance(0, 9, 12, 0, x, r) - r^2  # 弦への二乗距離

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

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (15/7, 39/7, 33/7)
    (90/17, 234/17, 198/17)

2 番めの解が妥当なものである。すなわち等円の半径は 15/7,直径は 30/7 である。

   r = 2.1428571,  x = 5.5714286;  y = 4.7142857

「天地明察」では,「勾股を掛け合わせ,これを二倍し,勾股弦の総和で割り,弦を掛け,勾股の和で割る」術が述べられている。SymPy では以下のようになる。

using SymPy
@syms 勾, 股
(勾, 股) = (9, 12)
弦 = Int(sqrt(勾^2 + 股^2))
2鉤*股//(鉤 + 股 + 弦)*弦//(鉤 + 股)

30/2

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 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")
   (r, x, y) = res[1]
   @printf("r = %.7f,  x = %.7f;  y = %.7f\n", r, x, y)
   plot([0, 12, 0, 0], [0, 0, 9, 0], color=:blue, lw=0.5)
   circle(r, y, r)
   circle(x, r, r)
   if more == true
       point(0, 9, " 9", :blue, :left, :bottom)
       point(12, 0, " 12", :green, :left, :bottom)
       point(r, y, "(r,y)", :red)
       point(x, r, "(x,r)", :red)
       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アクセスランキング にほんブログ村