算額（その1092）

算額（その1092）

五十一 岩手県一関市西風 西風日山神社 明治18年(1885)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html
キーワード：円2個，楕円1個，正三角形

正三角形の中に楕円 1 個と 2 個の等円を容れる。等円の直径が与えられるとき，楕円の長径を求めよ。

等円の半径と中心座標を r, (0, r), (0, 3r)
楕円の長半径と短半径を a, b; b = r
正三角形の一辺の長さを 2l
楕円と正三角形の辺の接点座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms l::positive, r::positive, a::positive, x0::positive, y0::positive
b = r
eq1 = r/(√Sym(3)l- 3r) - 1//2
eq2 = x0^2/a^2 + (y0 - r)^2/b^2 - 1
eq3 = -b^2*x0/(a^2*(y0 - r)) + √Sym(3)
eq4 = y0/x0 - 1/√Sym(3)
res = solve([eq1, eq2, eq3, eq4], (l, a, x0, y0))[1]

   (5*sqrt(3)*r/3, sqrt(5)*r, 5*sqrt(3)*r/4, 5*r/4)

楕円の長半径 a は，等円の半径 r の √5 倍である。
等円の直径が 1 寸のとき，楕円の長径は √5 = 2.23606797749979 寸である。

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

   r1 = 0.5;  l = 1.44338;  a = 1.11803;  x0 = 1.08253;  y0 = 0.625

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (l, a, x0, y0) = (5*sqrt(3)*r/3, sqrt(5)*r, 5*sqrt(3)*r/4, 5*r/4)
   @printf("等円の直径が %g のとき，楕円の長径は %g である。\n", 2r, 2a)
   @printf("r1 = %g;  l = %g;  a = %g;  x0 = %g;  y0 = %g\n", r, l, a, x0, y0)
   plot([l, 0, -l, l], [0, √3l, 0, 0], color=:blue, lw=0.5)
   ellipse(0, r, a, r, color= :green)
   circle(0, r, r)
   circle(0, 3r, r)
   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)", :green, :left, :vcenter)
       point(0, r, "等円:r,(0,r)", :red, :center, delta=-delta/2)
       point(0, 3r, "等円:r,(0,3r)", :red, :center, delta=-delta/2)
       point(a, r, "(a,r) ", :green, :right, :vcenter)
       point(0, 2r, "2r = 2b", :green, :center, delta=-delta/2)
       point(l, 0, " l", :blue, :left, :bottom, delta=delta/2)
       point(0, √3l, " √3l", :blue, :left, :vcenter)
   end
end;

算額（その1091）

算額（その1091）

四十三 岩手県一関市真滝 熊野白山滝神社 弘化3年(1846)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html
キーワード：円1個，四分円2個，正方形

正方形の中に四分円を 2 個，円を 1 個容れる。円の直径が与えられたときに黒積（灰色の部分の 2 倍）を求めよ。

正方形の一辺の長さを R，円の半径を r とおき，以下の方程式を解く。

include("julia-source.txt")

using SymPy

@syms R::positive, r::positive
eq1 = (R/2)^2 + r^2 - (R - r)^2
res = solve(eq1, R)[1]
res[1] |> println

   8*r/3

正方形の一辺の長さは円の半径の 8/3 倍である。

黒積の半分（図の灰色の部分）は，
扇形 OAR の面積（半径　R の円の面積の 1/6 = πR^2/6）から，
⊿AOC の面積（√3R/2 * R/2 / 2 = √3R^2/8）と，
円の面積の半分（πr^2/2）を引いたものの，
2倍である。

r = 1/2 のとき，R = 4/3 となり，黒積は 0.30648601314366886 である。

r = 1/2
R = 4/3
2(pi*R^2/6 - R^2*√3/8 - pi*r^2/2)

   0.30648601314366886

術は，円径^2*1.5593 とあるが，おかしい?

(2r)^2*1.5593

   1.5593

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   R = 8r/3
   黒積 = 2(pi*R^2/6 - R^2*√3/8 - pi*r^2/2)
   @printf("円の直径が %g(正方形の一辺の長さが %g) のとき，黒積は %g である。\n", 2r, R, 黒積)
   θ = 0:0.1:60
   x = @. R*cosd(θ)
   y = @. R*sind(θ)
   append!(x, [R/2, R])
   append!(y, [0, 0])
   plot([0, R, R, 0, 0], [0, 0, R, R, 0], color=:blue, lw=0.5)
   plot!(x, y, seriestype=:shape, fillcolor=:gray80)
   circlef(R/2, r, r, :white, beginangle=-90, endangle=90)
   circle(0, 0, R, beginangle=0, endangle=90)
   circle(R, 0, R, beginangle=90, endangle=180)
   circle(R/2, r, r, :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)
       segment(R/2, 0, R/2, √3R/2)
       segment(0, 0, R/2, √3R/2)
       point(R/2, √3R/2, "A", :black, :center, :bottom, delta=delta/2)
       point(R, 0, " R", :black, :left, :bottom, delta=delta/2)
       point(0, 0, "O ", :black, :right, :bottom, delta=delta/2)
       circle(0, 0, 0.05R,  :black, beginangle=0, endangle=60)
       circle(0, 0, 0.055R,  :black, beginangle=0, endangle=60)
       point(0.13R, 0.01R, "60°", :black, :right, :bottom, delta=delta/2, mark=false)
       point(R/2, r, " B", :black, :left, :vcenter)
       point(R/2, 0, " C", :black, :left, :bottom, delta=delta)
       plot!(xlims=(-3delta, R + 3delta))
   end
end;

 

モンテカルロ法で，解析解のチェック

 一辺の長さが R の正方形の中に，半径 R の四分円 2 個と，半径 r の円を描く。

四分円と円の隙間（緑の点が散らばっている領域）の面積をモンテカルロ法により求める。

1000000 個の点を打って，領域内に入った点の割合から面積を推定すると，0.172231*R^2 となった。解析解は 0.17239838239331384*R^2 だったので，解析解を導いた式は間違っていなさそうだ。

using Random, Distributions, StatsBase
R = 1
r = 3R/8
n = 1000000
x = rand(Uniform(0.5, 1.0), n)
y = rand(Uniform(0, 1.0), n)
z = @. (x^2 + y^2 < R^2) && ((x - R/2)^2 + (y - r)^2 > r^2)
2*(R^2/2 * mean(z))

   0.172231

算額（その1090）

算額（その1090）

四十三 岩手県一関市真滝 熊野白山滝神社 弘化3年(1846)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html
キーワード：円3個，二等辺三角形，正方形

正方形の中に二等辺三角形と等円 3 個を容れる。等円の直径が与えられたとき，正方形の一辺の長さを求めよ。

正方形の一辺の長さを a，
正方形の辺の上にある二等辺三角形の頂点座標を (b, 0), (a - b, 0)
等円の半径と中心座標を r, (r, r), (a - r, r), (a - r, a - r)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms a::positive, b::positive, r::positive
eq1 = a + b - sqrt(a^2 + b^2) - 2r
eq2 = 2(a - b) - sqrt(2(a - b)^2) - 2r
res = solve([eq1, eq2], (a, b))[1]

   (2*r*(-sqrt(2) + sqrt(26 - 16*sqrt(2)) + 2*sqrt(13 - 8*sqrt(2)))/(-2 - sqrt(2) + sqrt(2)*sqrt(13 - 8*sqrt(2)) + 2*sqrt(13 - 8*sqrt(2))), r*(-sqrt(2) + sqrt(26 - 16*sqrt(2)) + 2 + 2*sqrt(13 - 8*sqrt(2)))/2)

t = sqrt(13 - 8√2) とおいて，
正方形の一辺の長さ a は r*(6 + √2 + (2 + √2)t)/2
斜線と正方形の一辺の交点座標 b は r*(2 - √2 + (2 + √2)t)/2
である。

等円の半径が 1/2 寸のとき，正方形の一辺の長さは 2.9619546674751045 寸である。

r = 1/2
t = sqrt(13 - 8√2)
r*(6 + √2 + (2 + √2)t)/2

   2.9619546674751045

「術」も正しい。

等円径 = 1
位 = (√2 + 6)
(sqrt(4位 -10) + 位)*等円径/4

   2.961954667475105

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   t = sqrt(13 - 8√2)
   a = r*(6 + √2 + (2 + √2)t)/2
   b = r*(2 - √2 + (2 + √2)t)/2
   @printf("等円の直径が %g のとき，正方形の一辺の長さは %g である。\n", 2r, a)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
   plot!([0, b, a, 0], [a, 0, a - b, a], color=:green, lw=0.5)
   circle(r, r, r)
   circle(a - r, r, r)
   circle(a - r, a - r, r)
   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, r, "等円:r,(r,r)", :red, :center, delta=-delta/2)
       point(a - r, r, "等円:r,(a-r,r)", :red, :center, delta=-delta/2)
       point(a - r, a - r, "等円:r,(a-rr,a-rr)", :red, :center, delta=-delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, a, " a", :blue, :left, :bottom, delta=delta/2)
       point(b, 0, "b ", :green, :right, :bottom, delta=delta/2)
       point(a, a - b, "(a,a-b)  ", :green, :right, :vcenter)
   end
end;

算額（その1089）

算額（その1089）

四十三 岩手県一関市真滝 熊野白山滝神社 弘化3年(1846)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：楕円2個（楕円3個），正六角形（正三角形）

正六角形の中に等楕円（同じ大きさの楕円）を 2 個容れる。楕円の長径，短径が与えられたとき正六角形の一辺の長さを求めよ。

注：多分もとの算額の図でもそうなのであろうが，山村の図では楕円の書き方がまずく，図形を正確に表現できていない。正しくは下図のように，楕円は互いに接し，それぞれは正六角形の 2 辺にも接している。
算額では楕円を 2 個としているが正六角形の中心を点対称として，3 個の楕円を考えるとわかりやすい。しかも，第 3 の楕円は長径，短径が水平，垂直なので，計算しやすい。以下では，第 3 の楕円について解く。

原点を 第 3 の楕円の中心に置く。
楕円の長半径，短半径を a, b
正六角形の一辺の長さを l
楕円同士の接点（符号は違うが正六角形の一辺との接点と同じ）の座標を (x0, y0)
とおき，以下の連立方程式を解く。
作図する際には，(x0, y0) は不要である。

include("julia-source.txt")

using SymPy
@syms a::positive, b::positive, l::positive,
     x0::positive, y0::positive
@syms a, b, l, x0, y0
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + 1/√Sym(3)
eq3 = √Sym(3)y0 - √Sym(3)l/2 + x0
res = solve([eq1, eq2, eq3], (l, x0, y0))[2];

l，x0，y0 は以下のように，少し簡約化できる。

(2√Sym(3)sqrt(a^2 + 3b^2)/3, a^2/sqrt(a^2 + 3b^2), √Sym(3)b^2/sqrt(a^2 + 3b^2));

たとえば，長半径 が 4 寸，短半径が 2 寸のとき，正六角形の一辺の長さは 6.110100926607786 寸である。

a = 4
b = 2
2√3sqrt(a^2 + 3b^2)/3

   6.110100926607786

術は「邪術」である。これでは答えは出ない。
たとえば，長径と短径がともに 1 である（これは，円である）ときに以下のように 0.8611513361757016 になる。

長径 = 1
短径 = 1
位 = 長径^2 + 短径^2
A = (sqrt(位^2 + 2長径^2*短径^2) + 位)/6
sqrt(A)

   0.8611513361757016

図を描いてみると正六角形の一辺の長さは 0.8611513361757016 ではないことが明らかである。
円の半径が 1/2 なので，正三角形の一辺の長さは 2 になる。更に，正六角形の一辺の長さは，正三角形の一辺の長さの半分の 2/√3 倍で 1.1547005383792515 になる。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (4, 2)
   (l, x0, y0) = (2*sqrt(3)*sqrt(a^2 + 3*b^2)/3, a^2/sqrt(a^2 + 3*b^2), sqrt(3)*b^2/sqrt(a^2 + 3*b^2))
   @printf("l = %.15g;  x0 = %g;  y0 = %g\n", l, x0, y0)
   plot([√3l/2, 0, -√3l/2, √3l/2], [0, 3l/2, 0, 0], color=:blue, lw=0.5)
   plot!(l .* [0, √3/2, √3/2, 0, -√3/2, -√3/2, 0],
       l .* [-1/2, 0, 1, 3/2, 1, 0, -1/2], color=:green, lw=0.5)
   circle(0, l/2, l)
   ellipse(0, 0, a, b, color=:magenta)
   ellipse(l/2*cosd(30), l/2*sind(30) + l/2, a, b, φ=120, color=:magenta)
   ellipse(-l/2*cosd(30), l/2*sind(30) + l/2, a, b, φ=240, color=:magenta)
   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, 0, "", :magenta)
       point(l/2*cosd(30), l/2*sind(30) + l/2, "", :magenta)
       point(-l/2*cosd(30), l/2*sind(30) + l/2, "", :magenta)
       point(0, b, "b", :magenta, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :magenta, :left, :bottom)
   end
end;