裏 RjpWiki

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

算額(その707)

2024年02月19日 | Julia

算額(その707)

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード:円5個,外円,弦,斜線2本
#Julia, #SymPy, #算額, #和算

直径が 2.5 寸の外円内に弦と斜線を入れ,区画された領域に等円を 4 個入れる。等円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
弦と外円の交点座標を (x1, y1); y1 = 2r - R, x1 = sqrt(R^2 - y1^2)
斜線と外円の交点座標を (x2, y2); x2 = sqrt(R^2 - y2^2)
等円の半径と中心座標を r, (0, r - R), (0, R - r), (x, y1 + r)
とおき以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive,
     x1::positive, y1::positive,
     x2::positive, y2::positive
y1 = 2r - R
x1 = sqrt(R^2 - y1^2)
x2 = sqrt(R^2 - y2^2)
eq1 = dist(0, y1, x2, y2, 0, R - r) - r^2
eq2 = dist(0, y1, x2, y2, x, y1 + r) - r^2
eq3 = x^2 + (y1 + r)^2 - (R - r)^2;
# res = solve([eq1, eq2, eq3], (r, x, y2))

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (r, x, y2) = u
   return [
       -r^2 + (2*R - 3*r)^2*(R^2 - y2^2)*(R - 2*r + y2)^2/(R^2 - y2^2 + (R - 2*r + y2)^2)^2 + (2*R - 3*r - (2*R - 3*r)*(R - 2*r + y2)^2/(R^2 - y2^2 + (R - 2*r + y2)^2))^2,  # eq1
       -r^2 + (r - (r*(R - 2*r + y2) + x*sqrt(R^2 - y2^2))*(R - 2*r + y2)/(R^2 - y2^2 + (R - 2*r + y2)^2))^2 + (x - sqrt(R^2 - y2^2)*(r*(R - 2*r + y2) + x*sqrt(R^2 - y2^2))/(R^2 - y2^2 + (R - 2*r + y2)^2))^2,  # eq2
       x^2 + (-R + 3*r)^2 - (R - r)^2,  # eq3
   ]
end;

R = 25//20
iniv = BigFloat[0.5, 0.8, 1.1]
res = nls(H, ini=iniv)

   ([0.4734152343522919, 0.7576940667754868, 1.0587960931595892], true)

外円の直径が 2.5 寸のとき,等円の直径は 2*0.4734152343522919 = 0.9468304687045838 寸である。

その他のパラメータは以下のとおりである。

r = 0.473415;  x = 0.757694;  x1 = 1.21268;  y1 = -0.30317;  x2 = 0.664418;  y2 = 1.0588

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 25//20
   (r, x, y2) = res[1]
   y1 = 2r - R
   x1 = sqrt(R^2 - y1^2)
   x2 = sqrt(R^2 - y2^2)
   @printf("等円の直径 = %g\n", 2r)
   @printf("r = %g;  x = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g\n", r, x, x1, y1, x2, y2)
   plot()
   circle(0, 0 , R, :blue)
   circle(0, r - R, r)
   circle(x, y1 + r, r)
   circle(-x, y1 + r, r)
   circle(0, R - r, r)
   segment(-x1, y1, x1, y1, :green)
   segment(0, y1, x2, y2, :green)
   segment(0, y1, -x2, y2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x, y1 + r, " (x,y1+r)", :red, :left, :vcenter)
       point(0, R - r, " R-r", :red, :left, :vcenter)
       point(0, r - R, " r-R", :red, :left, :vcenter)
       point(0, y1, " y1", :green, :left, delta=-delta/2)
       point(x1, y1, "(x1,y1) ", :green, :right, delta=-delta/2)
       point(x2, y2, "(x2y2)", :green, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その706)

2024年02月19日 | Julia

算額(その706)

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード:円5個,外円,二等辺三角形
#Julia, #SymPy, #算額, #和算

直径 163 寸の外円内に圭(二等辺三角形),大円 2 個,小円 2 個が入っている。小円の直径はいかほどか。
注:後述するが 163 寸ではなく 162 寸である。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, r1 - R), (0, 3r1 - R)
小円の半径と中心座標を r2, (x2, 2r1 - R - r2)
二等辺三角形の底辺と外円の接点を (x1, y1); y1 = 2r1 - R, x1 = sqrt(R^2 - y1^2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive
eq1 = dist(sqrt(R^2 - (2r1 - R)^2), 2r1 - R, 0, R, 0, 3r1 - R) - r1^2
eq2 = x2^2 + (2r1 - R - r2)^2 - (R - r2)^2
eq3 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (r1, r2, x2))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (4*R/9, 20*R/81, 8*sqrt(5)*R/27)

小円の半径は 外円の半径の 20/81 倍である。

術では「外径を 20 倍し 80 で割る」とあるが,正しくは「81 で割る」であろう。
また,問にある「外径163寸」というのも「外径162寸」であろう。
この問に限らず,この算額にはミスが多い。

以上 2 点を修正すると,「外円の直径が 162 寸のとき,小円の直径は 40 寸」である。

その他のパラメータは以下のとおりである。

r1 = 36;  r2 = 20;  x2 = 53.6656

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 162//2  # 「問」のミスを修正
   (r1, r2, x2) = (4*R/9, 20*R/81, 8*sqrt(5)*R/27)
   @printf("小円の直径 = %g\n", 2r2)
   @printf("r1 = %g;  r2 = %g;  x2 = %g\n", r1, r2, x2)
   y1 = 2r1 - R
   x1 = sqrt(R^2 - y1^2)
   plot([x1, 0, -x1, x1], [y1, R, y1, y1], color=:black, lw=0.5)
   circle(0, 0 , R, :blue)
   circle(0, r1 - R, r1)
   circle(0, 3r1 - R, r1)
   circle(x2, 2r1 - R - r2, r2, :blue)
   circle(-x2, 2r1 - R - r2, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r1 - R, " 大円:r1\n (0,r1-R)", :red, :left, :vcenter)
       point(0, 3r1 - R, " 大円:r1\n (0,3r1-R)", :red, :left, :vcenter)
       point(x2, 2r1 - R - r2, "小円:r2\n(x2,2r1-R-r2)", :blue, :center, delta=-delta/2)
       point(x1, y1, "(x1,y1)  ", :black, :right, :bottom, delta=delta/2)
   end
end;

 

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

算額(その705)

2024年02月19日 | Julia

算額(その705)

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード:円4個,外円,二等辺三角形,弦,矢
#Julia, #SymPy, #算額, #和算

外円内に等斜(この場合は弦と矢),圭(二等辺三角形),等円 3 個をいれる。等斜が 1 寸のとき,等円の直径はいかほどか。

等斜の長さを a,圭の底辺の長さを 2b とする。
外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, a - R + r), (x, a - R + r)
として,いかの連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, x::positive, r::positive, R::positive
eq1 = x^2 + (a - R + r)^2 - (R - r)^2
eq2 = (a - R)^2 + (a/2)^2 - R^2
eq3 = dist(0, R, b, a - R, x, a - R + r) - r^2
eq4 = dist(0, R, b, a - R, 0, a - R + r) - r^2
solve([eq1, eq2, eq3, eq4], (r, x, b, R))

   1-element Vector{NTuple{4, Sym}}:
    (a*(3 - sqrt(5))/8, a*sqrt(-2 + sqrt(5))/2, a*(3 - sqrt(5))/(8*sqrt(-2 + sqrt(5))), 5*a/8)

等円の半径は a*(3 - sqrt(5))/8 である。

等斜を 1 寸とすれば,等円の半径は 0.09549150281252627 寸,直径は 0.19098300562505255 寸である。

その他のパラメータは以下のとおりである。

r = 0.0954915;  x = 0.242934;  b = 0.196538;  R = 0.625

a = 1
a*(3 - sqrt(5))/8, 2*a*(3 - sqrt(5))/8

   (0.09549150281252627, 0.19098300562505255)

「答曰 等円径三分八厘二毛」とあるが,これは直径を二倍してしまっているようだ。

「術曰」では,「置五個開平方以滅▢個内餘乗等斜四除之得等径合問」とある。
欠損している文字は「三」で,「三から五の平方根を引き等斜を掛け四で割る」と読める。上述の式の通り a*(3 - sqrt(5))/4 ということで,直径を表している。これで1分9厘0毛9糸なので,それをまた2倍して3分8厘2毛としてしまっている。

a*(3 - sqrt(5))/4, 2*a*(3 - sqrt(5))/4

   (0.19098300562505255, 0.3819660112501051)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1
   (r, x, b, R) = (a*(3 - sqrt(5))/8, a*sqrt(-2 + sqrt(5))/2, a*(3 - sqrt(5))/(8*sqrt(-2 + sqrt(5))), 5*a/8)
   println("等円の直径 = $(2r)")
   @printf("r = %g;  x = %g;  b = %g;  R = %g\n", r, x, b,R)
   plot()
   circle(0, 0 , R, :blue)
   circle(0, a - R + r, r)
   circle(x, a - R + r, r)
   circle(-x, a - R + r, r)
   segment(0, R, b, a - R)
   segment(0, R, -b, a - R)
   segment(-a/2, sqrt(R^2 - a^2/4), a/2, sqrt(R^2 - a^2/4))
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, a - R, " a-R", :black, :left, delta=-delta/2)
       point(b, a - R, "b", :black, :center, delta=-delta/2)
       point(x, a - R + r, " r,(x,a-R+r)", :black, :left, :vcenter)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その704)

2024年02月19日 | Julia

算額(その704)

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード:四分円4個,正方形,面積
#Julia, #SymPy, #算額, #和算

一辺の長さが 3.5 寸の正方形内に 4 つの四分円を描く。黒積(図に色付した部分の面積の 8 倍)を求めよ。

この問題は,算額(その122)と同趣旨のものである。

正方形の一辺の長さを a とする。

図に色付した部分の面積は,正方形の面積の 1/2 から,半径 a の円の面積の 1/12 と底辺の長さ a/2,高さ √3a/2 の三角形の面積を差し引いたものである(ACDO - ABO - BDO)。
黒積はその 8 倍である。

include("julia-source.txt");

a = 3.5
8(a^2/2 - pi*a^2/12 - √3*a^2/8)

   2.1260376029646118

SymPy の integrate() でも求めてみる。

using SymPy
@syms a, x
a = 3.5
8*integrate(a - sqrt(a^2 - x^2), (x, 0, a/2)) |> println

   2.12603760296460

黒積は約 2.126 平方寸である。

術では 2.195 としているが,円周率として 3.13314827844261 という不思議な値を使っていることになる。

@syms 円周率
eq = 8(a^2/2 - 円周率*a^2/12 - √3*a^2/8) - 2.195
solve(eq, 円周率)[1] |> println

   3.13314827844261

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 3.5
   θ = 60:0.01:90
   x = vcat(a .* cosd.(θ), [a/2, a/2])
   y = vcat(a .* sind.(θ), [a, √3a/2])
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circle(0, 0 , a, beginangle=0, endangle=90)
   circle(a, 0 , a, beginangle=90, endangle=180)
   circle(a, a , a, beginangle=180, endangle=270)
   circle(0, a , a, beginangle=270, endangle=360)
   segment(0, a/2, a, a/2)
   segment(a/2, 0, a/2, a)
   segment(0, 0, a/2, √3a/2)
   plot!(x, y, seriestype=:shape, color=:magenta, fillcolor=:magenta, alpha=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, a, " A", :green, :left, :bottom, delta=delta/2)
       point(a/2, √3*a/2, " B", :green, :left, :bottom, delta=1.5delta)
       point(a/2, a, " C", :green, :center, :bottom, delta=delta/2)
       point(a/2, 0, " D", :green, :left, :bottom, delta=delta/2)
       point(0, 0, "  O", :green, :left, :bottom, delta=delta/4)
   end
end;

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

算額(その703)

2024年02月19日 | Julia

算額(その703)

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード:円3個,斜線,弦,矢
#Julia, #SymPy, #算額, #和算

大円 2 個と小円 1 個が直線上に載っている。小円に接する斜線を引き,大円 2 個の弦とする。弦の中心と弧の中心を「矢」というが,直径から矢の長さを引いたものを甲矢,乙矢とする。甲矢,乙矢の長さが与えられたとき,大円の直径を求めよ。

大円と小円が互いに接しまた直線に接していることから,
r1^2 + (r1 - r2)^2 = (r1 + r2) より,r2 = r1/4 である。

大円の中心から斜線までの距離は,甲矢,乙矢の長さから半径を引いたものと一致する。

斜線の x 切片,y 切片を (a, 0), (0, b)
大円の半径と中心座標を r1, (r1, r1), (-r1, r1)
小円の半径と中心座標を r2, (0, r2)
甲矢,乙矢を「甲矢」,「乙矢」
として,以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms 甲矢::positive, 乙矢::positive,
     a::positive, b::positive,
     r1::positive, r2::positive
r2 = r1/4
eq1 = dist(a, 0, 0,  b, r1, r1) - (甲矢 - r1)^2
eq2 = dist(a, 0, 0,  b, -r1, r1) - (乙矢 - r1)^2
eq3 = dist(a, 0, 0,  b, 0, r2) - r2^2;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (r1, a, b) = u
   return [
       -(-r1 + 甲矢)^2 + (-b*(-a*(-a + r1) + b*r1)/(a^2 + b^2) + r1)^2 + (-a + a*(-a*(-a + r1) + b*r1)/(a^2 + b^2) + r1)^2,  # eq1
       -(-r1 + 乙矢)^2 + (-b*(-a*(-a - r1) + b*r1)/(a^2 + b^2) + r1)^2 + (-a + a*(-a*(-a - r1) + b*r1)/(a^2 + b^2) - r1)^2,  # eq2
       -r1^2/16 + (-a + a*(a^2 + b*r1/4)/(a^2 + b^2))^2 + (-b*(a^2 + b*r1/4)/(a^2 + b^2) + r1/4)^2,  # eq3
   ]
end;

(甲矢, 乙矢) = (15, 10)
iniv = BigFloat[9, 12, 4.7]
res = nls(H, ini=iniv)

   ([8.520833333333334, 14.20138888888889, 4.35848252344416], true)

甲矢 15 寸,乙矢 10 寸のとき,大円の直径は 17.0417 寸である。

その他のパラメータは以下のとおりである。

r1 = 8.52083;  r2 = 2.13021;  a = 14.2014;  b = 4.35848

function 矢(x0, y0, r, a, b)
   x1 = (a*y0 + x0 - a*b)/(a^2 + 1)
   y1 = a*x1 + b
   (x2, y2) = (a*y0 - a*(r/sqrt(a^2 + 1) + y0) + x0, r/sqrt(a^2 + 1) + y0)
   return (x1, y1, x2, y2)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, a, b) = res[1]
   甲矢 = r1 + sqrt(dist(a, 0, 0,  b, r1, r1))
   乙矢 = r1 + sqrt(dist(a, 0, 0,  b, -r1, r1))
   r2 = r1/4
   @printf("甲矢 = %g;  乙矢 = %g;  大円の直径 = %g\n", 甲矢, 乙矢, 2r1)
   @printf("r1 = %g;  r2 = %g;  a = %g;  b = %g\n", r1, r2, a, b)
   plot()
   circle(r1, r1, r1)
   circle(-r1, r1, r1)
   circle(0, r2, r2, :blue)
   abline(0, b, -b/a, -2r1, 2r1)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :top, delta=-2delta)
       point(-r1, r1, " (-r1,r1)", :red, :left, :vcenter)
       point(r1, r1, " (r1,r1)", :red, :left, :vcenter)
       point(0, r2, " r2", :blue, :left, :vcenter)
       (x1, y1, x2, y2) = 矢(r1, r1, r1, -b/a, b)
       segment(x1, y1, x2, y2, :green, lw=1)
       point(10, 12, "甲矢",  mark=false)
       (x3, y3, x4, y4) = 矢(-r1, r1, r1, -b/a, b)
       segment(x3, y3, x4, y4, :green, lw=1)
       point(-10.5, 12, "乙矢",  mark=false)
   end
end;

 

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

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

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