裏 RjpWiki

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

算額(その719)

2024年02月22日 | Julia

算額(その719)

四二 浦和市西堀(現さいたま市桜区西堀) 氷川神社 嘉永5年(1852)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

埼玉県浦和市西堀 氷川神社 嘉永5年
http://www.wasan.jp/saitama/uhikawa.html

直角三角形内に斜線を引き一方の区画に内接する円を描く。直角を挟む二辺の短い方(鈎)と長い方(股)の辺の長さがそれぞれ 8 寸,15 寸,円の直径が 4 寸のとき,斜線の長さはいかほどか。

鈎,股をそれぞれ「鈎」,「股」
斜線と股の交点座標を (a, 0)
円の半径と中心座標を r, (x, r)
とおき,以下の連立方程式を解く。
斜線の長さ length は後で計算してもよいが,一緒に求めておく。

include("julia-source.txt");

using SymPy
@syms r::positive, x::positive, a::positive,
     length::positive,
     鈎::positive, 股::positive
eq1 = dist(0, 鈎, a, 0, x, r) - r^2
eq2 = dist(0, 鈎, 股, 0, x, r) - r^2
eq3 = sqrt(鈎^2 + a^2) - length
res = solve([eq1, eq2, eq3], (a, x, length));
("a:", res[3][1]) |> println  # a
("x:", res[3][2]) |> println  # x
("length:", res[3][3]) |> println  # length

   ("a:", (-2*r^2*股 - 2*r^2*sqrt(股^2 + 鈎^2) + 2*r*股*鈎 + 2*r*鈎*sqrt(股^2 + 鈎^2) - 股*鈎^2)/(鈎*(2*r - 鈎)))
   ("x:", (-r*sqrt(股^2 + 鈎^2) - 股*(r - 鈎))/鈎)
   ("length:", sqrt(鈎^4*(2*r - 鈎)^2 + (2*r^2*股 + 2*r^2*sqrt(股^2 + 鈎^2) - 2*r*股*鈎 - 2*r*鈎*sqrt(股^2 + 鈎^2) + 股*鈎^2)^2)/(鈎*Abs(2*r - 鈎)))

4 組の解が得られるが,3 番目のものが適解である。
sqrt(股^2 + 鈎^2) を 「弦」とおくと少しきれいになるか。
a: (-2*r^2*股 - 2*r^2*弦 + 2*r*股*鈎 + 2*r*鈎*弦 - 股*鈎^2)/(鈎*(2*r - 鈎))
x: (-r*弦 - 股*(r - 鈎))/鈎
length: sqrt(鈎^4*(2*r - 鈎)^2 + (2*r^2*股 + 2*r^2*弦 - 2*r*股*鈎 - 2*r*鈎*弦 + 股*鈎^2)^2)/(鈎*Abs(2*r - 鈎)))

鈎,股,円の半径が 8,15,4/2 のとき斜線の長さは 10 寸である。
他のパラメータは a = 6,x = 7 である。

(鈎, 股) = (8, 15)
r = 4//2
Abs = abs
((-2*r^2*股 - 2*r^2*sqrt(股^2 + 鈎^2) + 2*r*股*鈎 + 2*r*鈎*sqrt(股^2 + 鈎^2) - 股*鈎^2)/(鈎*(2*r - 鈎)), (-r*sqrt(股^2 + 鈎^2) - 股*(r - 鈎))/鈎, sqrt(鈎^4*(2*r - 鈎)^2 + (2*r^2*股 + 2*r^2*sqrt(股^2 + 鈎^2) - 2*r*股*鈎 - 2*r*鈎*sqrt(股^2 + 鈎^2) + 股*鈎^2)^2)/(鈎*Abs(2*r - 鈎)))

   (6.0, 7.0, 10.0)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, r) = (8, 15, 4/2)
   Abs = abs
   (a, x, length) = ((-2*r^2*股 - 2*r^2*sqrt(股^2 + 鈎^2) + 2*r*股*鈎 + 2*r*鈎*sqrt(股^2 + 鈎^2) - 股*鈎^2)/(鈎*(2*r - 鈎)), (-r*sqrt(股^2 + 鈎^2) - 股*(r - 鈎))/鈎, sqrt(鈎^4*(2*r - 鈎)^2 + (2*r^2*股 + 2*r^2*sqrt(股^2 + 鈎^2) - 2*r*股*鈎 - 2*r*鈎*sqrt(股^2 + 鈎^2) + 股*鈎^2)^2)/(鈎*Abs(2*r - 鈎)))
   @printf("斜線の長さ = %g;  a = %g;  x = %g\n", length, a, x)
   plot([0, 股, 0,  0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   circle(x, r, r)
   segment(0, 鈎, a, 0, :green)
   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(0, 鈎, " 鈎", :blue, :left, :bottom)
       point(股, 0, " 股", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, "a ", :blue, :right, :bottom, delta=delta/2)
       point(x, r, "r,(x,r)", :red, :center, delta=-delta/2)
   end
end;

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

算額(その718)

2024年02月22日 | Julia

算額(その718)

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

長方形内に菱形と大円 2 個,小円 4 個をいれる。小円は長方形の二辺に内接し菱形と大円に外接している。
長方形の短辺が寸五分のとき,小円の直径はいかほどか。

長方形の中心を原点として,長辺と短辺を 2a, 2b とする。
大円の半径と中心座標を r1, (0, r1)
小円の半径と中心座標を r2, (a - r1, b - r2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     r1::positive, r2::positive
r1 = b/2
eq1 = (a - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = dist(0, b, a, 0, a - r2, b - r2) - r2^2
res = solve([eq1, eq2], (r2, a))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (b*(9 - sqrt(17))/16, b*(5 + 3*sqrt(17))/16)

「問」では「平寸五分」と一文字が欠損しているがこれは「二」であることがわかる。
小円の直径は 2(短辺/2*(9 - √17)/16 = 0.7620147459972405 つまり 「七分六厘二毛有奇」

短辺 = 2.5  # 二寸五分 = 2b
2(短辺/2*(9 - √17)/16)

   0.7620147459972405

「術」は,「から十七の平方根を引いたものに,短辺を十六で割ったものを掛ける」とあるが,五は九の誤記である。

("誤", (5 - sqrt(17))*(短辺/16)) |> println
("正", (9 - sqrt(17))*(短辺/16)) |> println

   ("誤", 0.13701474599724053)
   ("正", 0.7620147459972405)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = 2.5/2
   r1 = b/2
   (r2, a) = (b*(9 - sqrt(17))/16, b*(5 + 3*sqrt(17))/16)
   @printf("小円の直径 = %g\n", 2r2)
   @printf("r2 = %g, a = %g\n", r2, a)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:gray70, lw=0.5)
   plot!([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:green, lw=0.5)
   circle(0, r1, r1)
   circle(0, -r1, r1)
   circle4(a - r2, b - r2, r2, :blue)
   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(0, r1, " 大円:r1\n (0,r1)", :red, :left, :vcenter)
       point(a - r2, b - r2, " 小円:r2\n (a-r2,b-r2)", :blue, :center, delta=-delta/2)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その717)

2024年02月22日 | Julia

算額(その717)

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

外円内に正三角形と弦を描き,区画された領域に菱形 1 個,甲円 1 個,乙円 3 個,丙円 2 個を入れる。
菱形の長い方の対角線が 1 寸 5 分のとき,丙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, r1 - R)
乙円の半径と中心座標を r2, (x2, y + r2), (0, y + r2)
丙円の半径と中心座標を r3, (x3, 2r1 - R - r3)
菱形の長い方の対角線の長さを 2x
乙円が接している弦と y 軸との交点座標を (0, y) とする。
R = 4r1, r3 = 3r1/4, x3 = 2sqrt(r1*r3) は容易にわかるので,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, x::positive, y::positive,
     r1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive
@syms d
R = 4r1
r3 = 3r1/4
x3 = 2*sqrt(r1*r3)
eq1 = x3^2 + (r3 - r1)^2 - (r1 + r3)^2
eq2 = x3^2 + (2r1 - R - r3)^2 - (R - r3)^2
eq3 = x2^2 + (y + r2)^2 - (R - r2)^2
eq3 = eq3 |> expand
eq4 = dist(sqrt(Sym(3))/2*R, -R/2, 0, R, x2, y + r2) - r2^2
eq4 = eq4 |> expand
eq5 = dist(sqrt(Sym(3))/2*R, -R/2, 0, R, 0, y + r2) - r2^2
eq5 = apart(eq5, d);
eq6 = ((y + 2r1 - R)/2 - 2r1 + R)/(sqrt(Sym(3))/2*R - x) - sqrt(Sym(3))
eq6 = eq6 |> factor;
eq7 = 16R*r3 - 3R^2;

res = solve([eq3, eq4, eq5, eq6], (r1, r2, x2, y))

   1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (25*sqrt(3)*x/111, 8*sqrt(3)*x/37, 32*x/37, 28*sqrt(3)*x/111)

丙円の直径は,甲円の直径の 3/4 倍であり,甲円の直径は 2(25*sqrt(3)*x/111) である
菱形の長い方の対角線の長さが 2x のとき,丙円の直径は 2(25*sqrt(3)*x/111)*(3/4) である。

2x = 1.5 寸のとき 0.4388642248907628 すなわち,「4分3厘8毛有奇」である。

x = 1.5/2
2(25*sqrt(3)*x/111)*(3/4)

   0.4388642248907628

「術」は,「3 の平方根に菱長と 25 を掛けて 148 で割る」sqrt(3)*1.5*25/148 = 0.4388642248907628 である。

sqrt(3)*1.5*25/148

   0.4388642248907628

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

R = 1.1703, r1 = 0.292576;  r2 = 0.280873;  x2 = 0.648649;  r3 = 0.219432;  x3 = 0.506757

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x = 15//20
   (r1, r2, x2, y) = (25*sqrt(3)*x/111, 8*sqrt(3)*x/37, 32*x/37, 28*sqrt(3)*x/111)
   R = 4r1
   r3 = 3r1/4
   x3 = 2*sqrt(r1*r3)
   @printf("丙円の直径 = %g\n", 2r3)
   @printf("R = %g, r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", R, r1, r2, x2, r3, x3)
   plot([√3/2*R, 0, -√3/2*R, √3/2*R], [-R/2, R, -R/2, -R/2], color=:gray70, lw=0.5)
   circle(0, 0, R)
   circle(0, r1 - R, r1)
   circle(x3, 2r1 - R - r3, r3, :blue)
   circle(-x3, 2r1 - R - r3, r3, :blue)
   circle(0, y + r2, r2, :orange)
   circle(x2, y + r2, r2, :orange)
   circle(-x2, y + r2, r2, :orange)
   plot!([x, 0, -x, 0, x], [(y + 2r1 - R)/2, y, (y + 2r1 - R)/2, 2r1 - R, (y + 2r1 - R)/2], color=:green, lw=0.5)
   xa = sqrt(R^2 - y^2)
   segment(-xa, y, xa, y)
   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(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(0, y, " y", :black, :left, :top, delta=-delta)
       point(0, y + r2, " 乙円:r2\n (0,y+r2)", :black, :left, :vcenter)
       point(x2, y + r2, " 乙円:r2\n (x2,y+r2)", :black, :center, delta=-delta/2)
       point(x, (y + 2r1 - R)/2, "(x,(y+2r1-R)/2) ", :green, :right, :vcenter)
       point(0, 2r1 - R, " 2r1-R", :red, :center, :bottom, delta=delta/2)
       point(0, r1 - R, "甲円:r1\n(0,r1-R)", :red, :center, :bottom, delta=delta/2)
       point(x3, 2r1 - R - r3, "丙円:r1\n(x3,2r1-R-r3)", :black, :center, :bottom, delta=delta/2)
   end
end;

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

算額(その716)

2024年02月22日 | Julia

算額(その716)

埼玉県秩父市大宮 秩父神社 明治20年

山口正義(2015): やまぶき, 第27号
https://yamabukiwasan.sakura.ne.jp/ymbk27.pdf

楕円の中に大円4個,小円3個が入っている。小円の直径が 1 寸のとき,楕円の長径はいかほどか。

楕円の長半径,短半径を a, b,中心座標を (0, 0); b = 3r2
大円の半径と中心座標を r1, (x1, 0), (0, r2); r1 = 2r2
楕円と大円の接点座標を (x0, y0)
小円の半径を r2
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, r2::positive,
     a::positive, b::positive,
     x0::positive, y0::positive
r1 = 2r2
b = 3r2
eq1 = x1^2 + r2^2 - 4r1^2
eq2 = (x0 - x1)^2 + y0^2 - r1^2
eq3 = x0^2/a^2 + y0^2/b^2 - 1
eq4 = -b^2*x0/(a^2*y0) - (x1 - x0)/y0
res = solve([eq1, eq2, eq3, eq4], (a, x1, x0, y0))

   1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (6*r2, sqrt(15)*r2, 4*sqrt(15)*r2/3, sqrt(21)*r2/3)

楕円の長半径は小円の半径の 6 倍である。
小円の半径が 1/2 寸のとき,長半径は 3 寸(長径は 6 寸)である。

r2 = 1//2
(6.0*r2, sqrt(15)*r2, 4*sqrt(15)*r2/3, sqrt(21)*r2/3)

   (3.0, 1.9364916731037085, 2.581988897471611, 0.7637626158259733)

楕円の長径 = 6;  a = 3;  x1 = 1.93649;  x0 = 2.58199;  y0 = 0.763763

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   r1 = 2r2
   b = 3r2
   (a, x1, x0, y0) = (6*r2, sqrt(15)*r2, 4*sqrt(15)*r2/3, sqrt(21)*r2/3)
   @printf("楕円の長径 = %g;  a = %g;  x1 = %g;  x0 = %g;  y0 = %g\n", 2a, a, x1, x0, y0)
   plot()
   circle(x1, 0, r1)
   circle(-x1, 0, r1)
   circle.(0, [0, 2r2, -2r2], r2)
   circle(0, r2, r1, :green)
   circle(0, -r2, r1, :green)
   ellipse(0, 0, a, b, color=:blue)
   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(x0, y0, "(x0,y0)", :blue, :left, :bottom, delta=delta/2)
       point(x1, 0, "x1", :magenta, :center, :bottom, delta=delta/2)
       point(0, r2, " r2", :green, :center, delta=-delta/2)
       point(0, b, " b", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
   end
end;

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

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

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