算額（その1098）

算額（その1098）

六十六 岩手県花泉町金沢字大柳 金沢八幡宮 明治29年(1896)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円3個，四分円2個，正三角形

正三角形と 2 個の四分円が重なっている領域に等円 3 個を容れる。等円の直径が与えられたときに正三角形の一辺の長さを求めよ。

正三角形の一辺の長さを 2a
四分円の半径と中心座標を R, (a, 0), (-a, 0)
等円の半径と中心座標を r, (a - r, y1), (0, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms a::positive, R::positive, r::positive, y1::positive, y2::positive
eq1 = r^2 + y1^2 - (R - r)^2
eq2 = a^2 + y2^2 - (R + r)^2
eq3 = 2r - (√Sym(3)a -y2)
eq4 = dist2(0, √Sym(3)a, a, 0, a - r, y1, r)
res = solve([eq1, eq2, eq3, eq4], (a, R, y1, y2))[1];

三角形の一辺の長さは等円の半径の (sqrt(3) + sqrt(4*sqrt(2) + 4*sqrt(3) + 4*sqrt(6) + 11)) 倍である。

a = res[1] |> sympy.sqrtdenest |> simplify
a |> println

   r*(sqrt(3) + sqrt(4*sqrt(2) + 4*sqrt(3) + 4*sqrt(6) + 11))/2

R = res[2] |> sympy.sqrtdenest |> simplify
R |> println

   r*(1 + sqrt(2) + sqrt(6))

y1 = res[3] |> sympy.sqrtdenest |> simplify
y1 |> println

   r*(sqrt(3) + 2)

y2 = res[4] |> sympy.sqrtdenest |> simplify
y2 |> println

   r*(-1 + sqrt(12*sqrt(2) + 12*sqrt(3) + 12*sqrt(6) + 33))/2

等円の直径が 1 寸のとき，三角形の一辺の長さは 3.7549272907866866 寸である。

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

   r = 0.5,  a = 1.87746;  R = 2.43185;  y1 = 1.86603;  y2 = 2.25186

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   a = r*(√3 + sqrt(4√2 + 4√3 + 4√6 + 11))/2
   R = r*(1 + √2 + √6)
   y1 = r*(√3 + 2)
   y2 = r*(-1 + sqrt(12√2 + 12√3 + 12√6 + 33))/2
   @printf("等円の直径が %g のとき，正三角形の一辺の長さは %g である。\n", 2r, 2a)
   @printf("r = %g,  a = %g;  R = %g;  y1 = %g;  y2 = %g\n", r, a, R, y1, y2)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:green, lw=0.5)
   circle(0, y2, r)
   circle2(a - r, y1, r)
   circle(a, 0, R, :blue, beginangle=90, endangle=180)
   circle(-a, 0, R, :blue, beginangle=0, endangle=90)
   segment(a, 0, a, R, :blue)
   segment(-a, 0, -a, 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)
       point(0, y2, "等円:r,(0,y2)", :red, :center, delta=-delta/2)
       point(a - r, y1, "等円:r,(a-r,y1)", :red, :center, delta=-delta/2)
       point(0, √3a, " √3a", :green, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
       point(a, R, "(a,R)", :blue, :right, :bottom, delta=delta/2)
   end
end;

算額（その1097）

算額（その1097）

六十六 岩手県花泉町金沢字大柳 金沢八幡宮 明治29年(1896)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円11個（円12個），外円，正五角形，五芒星

算額（その1088）にプラスαした図形である。
https://blog.goo.ne.jp/r-de-r/e/8116c7b02599138383dcfbe11d21ff27

正五角形の対角線を引き，区画された領域の中央に円 1 個，小円 5 個，大円 5 個を容れる。小円の直径与えられたとき，大円の直径を求めよ。

正五角形が内接する円の半径と中心座標を R, (0, 0)
中央の円の半径と中心座標を r1, (0, 0)
小円の半径と中心座標を r2, (0, r1 + r2)
大円の半径と中心座標を r3, (0, r3 - R)
とし，以下のように逐次決定してゆく。

計算において，必要な角度がいくつか出てくる。その三角関数の値は無理数ではあるが，きれいな形で表現される。

include("julia-source.txt")

using SymPy
@syms R
s18 = Sym(18)
s36 = Sym(36);

1.　中央の円の半径

中央の円の半径 OB は五芒星の頂点（正五角形の頂点） A の y 座標値に等しい。
∠AOR は 18° である。
中央の円の半径は R*(√5 - 1)/4 である。

r1 = y = R*sind(s18)
r1 |> println

   R*(-1/4 + sqrt(5)/4)

2.　小円の半径

小円の半径 BD は BC * tan(∠BCD) である。
BC は OB * tan(∠BOC) である。
ここで，∠BCD = ∠BOC = 36°である。
まとめると，小円の半径は R*(7√5 - 15)/4 である。

OB = R*sind(s18)
r2 = OB*tand(s36)*tand(s36)
r2 |> simplify |> println

   R*(-15 + 7*sqrt(5))/4

3.　大円の半径

大円の半径 GF は EF * tan(∠FEG) である。
EF は OE * sin(∠EOF) である。
ここで，∠FEG = 18°，∠EOF = 36°である。
まとめると，大円の半径は R\*sqrt(14 - 6\*sqrt(5))/4 である。
なお，中央の円を含む三角形と大円を含む三角形の相似比から求めるのやり方もある。

EF = R*sind(s36)
r3 = EF*tand(s18)
r3 |> simplify |> println

   R*sqrt(14 - 6*sqrt(5))/4

@syms d
apart(r3/r2, d) |> simplify |> println

   1/2 + 3*sqrt(5)/10

1/2 + 3*sqrt(5)/10

   1.170820393249937

4.　小円の直径が既知のときの大円の直径

大円の直径は小円の直径の r1/r2 = 1/2 + 3*sqrt(5)/10 = 2√5/5 + 1 = 1.170820393249937 倍である。
小円の直径が 1 寸のとき，大円の直径は 1.170820393249937 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 1
   r1 = R*(√5 - 1)/4
   r2 = R*(7√5 - 15)/4
   r3 = R*sqrt(14 - 6√5)/4
   θ = 90:-72:-342
   x = R .* cosd.(θ)
   y = R .* sind.(θ)
   plot(x[1:6], y[1:6], seriestype=:shape, fillcolor=:aliceblue, lw=0.5)
   circle(0, 0, R)
   for i in 1:5
       segment(x[i], y[i], x[i + 1], y[i + 1], :black)
       segment(x[i], y[i], x[i + 2], y[i + 2], :black)
   end
   circlef(0, 0, r1, :red)
   circle(0, 0, r1, :black)
   rotatef(0, r1 + r2, r2, :yellow, angle=72)
   rotate(0, r1 + r2, r2, :black, angle=72)
   rotatef(0, r3 - R*cosd(36), r3, :orange, angle=72)
   rotate(0, r3 - R*cosd(36), r3, :black, angle=72)
   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", :black, :left, :bottom, delta=delta/2)
       point(0, r1, "B", :black, :center, delta=-delta/2)
       point(0, r1 + r2, " D", :black, :left, :vcenter)
       point(0, r1 + r2, "小円", :black, :center, delta=4delta)
       point(0, r3 - R*cosd(36), " G", :black, :left, :vcenter)
       point(0, r3 - R*cosd(36), "大円", :black, :center, delta=4delta)
       point(0, 0, "O", :black, :center, delta=-delta/2)
       point(0, 0, "円", :black, :center, delta=-3delta)
       point(r1*tand(36), r1, "C", :black, :left, :bottom, delta=delta/2)
       point(R*cosd(18), r1, " A", :black, :left, :vcenter)
       point(R*sind(36), -R*cosd(36), "E", :black, :left, delta=-delta/2)
       point(0, -R*cosd(36), "F", :black, :center, delta=-delta/2)
   end
end;

算額（その1096）

算額（その1096）

五十八 岩手県花泉町 花泉天満宮 明治3年(1870)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円5個，二等辺三角形

二等辺三角形の頂点を周上に持つ甲円と，二等辺三角形の中に乙円，丁円，丙円を容れる。乙円と丁円の直径がそれぞれ 15 寸，10 寸のとき，丙円の直径はいかほどか。

注：山村（算額も？）の図では丙円と甲円が離れているように描かれているが，それでは解は不定である。また，「答」にあるように丙円の直径が 16 寸のときには，算額のような図にはならない。算額に似るように描いた下図は甲円，丁円の直径が 10 寸，5 寸のときのものである。

二等辺三角形の底辺の長さを 2a，高さを b
甲円の半径と中心座標を r1, (0, 2r4 + r1)
乙円の半径と中心座標を r2, (0, 2r4 + r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (0, r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, r4::positive
eq1 = b - 2r1 - 2r4
#eq2 = r2/(b - r2 - 2r4) - a/sqrt(a^2 + b^2)
eq3 = dist2(0, b, a, 0, x3, r3, r3)
eq4 = dist2(0, b, a, 0, 0, 2r4 + r2, r2)
eq5 = x3^2 + (r4 - r3)^2 - (r3 + r4)^2 |> expand
eq6 = x3^2 + (2r4 + r1 - r3)^2 - (r1 + r3)^2 |> expand;
#res = solve([eq1, eq3, eq4, eq5, eq6], (a, b, r1, r3, x3))

println(eq1, ",  # eq1")
#println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")
println(eq6, ",  # eq6")

   b - 2*r1 - 2*r4,  # eq1
   b*(a^2*b - 2*a^2*r3 - 2*a*b*x3 + 2*a*r3*x3 - b*r3^2 + b*x3^2),  # eq3
   a^2*b^2 - 2*a^2*b*r2 - 4*a^2*b*r4 + 4*a^2*r2*r4 + 4*a^2*r4^2 - b^2*r2^2,  # eq4
   -4*r3*r4 + x3^2,  # eq5
   -4*r1*r3 + 4*r1*r4 - 4*r3*r4 + 4*r4^2 + x3^2,  # eq6

eq1, eq5 をそれぞれ解いて，r1, r3 を求める。

ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println

   b/2 - r4

ans_r3 = solve(eq5, r3)[1]
ans_r3 |> println

   x3^2/(4*r4)

eq3, eq4, eq6 に r1, r3 を代入し整理する。

eq3 = eq3(r1 => ans_r1, r3 => ans_r3) |> simplify |> (x -> x*16r4^2/b)
eq4 = eq4(r1 => ans_r1, r3 => ans_r3)
eq6 = eq6(r1 => ans_r1, r3 => ans_r3)*2r4 |> simplify;

eq6 を解いて b を求める。

ans_b = solve(eq6, b)[1] |> simplify
ans_b |> println

   2*r4*x3^2/(-4*r4^2 + x3^2)

eq3, eq4 に b を代入し整理する。

eq3 = eq3(b => ans_b) |> simplify |> (x -> x*(4r4^2 - x3^2)) |> expand |> simplify |> (x -> x/(2r4*x3^2))
eq4 = eq4(b => ans_b) |> simplify |> numerator |> (x -> x/4r4^2);

eq3, eq4 を解いて，a, x3 を求める。

res = solve([eq3, eq4], (a, x3))[4]
res[1] |> println

   4*r2*r4^(3/2)*sqrt(r2 + r4)*(3*r2^16 + 44*r2^15*r4 + 300*r2^14*r4^2 + 1260*r2^13*r4^3 + 3640*r2^12*r4^4 + 7644*r2^11*r4^5 + 12012*r2^10*r4^6 + 14300*r2^9*r4^7 + 12870*r2^8*r4^8 + 8580*r2^7*r4^9 + 4004*r2^6*r4^10 + 1092*r2^5*r4^11 - 140*r2^3*r4^13 - 60*r2^2*r4^14 - 12*r2*r4^15 - r4^16)/(3*r2^18 + 44*r2^17*r4 + 297*r2^16*r4^2 + 1216*r2^15*r4^3 + 3340*r2^14*r4^4 + 6384*r2^13*r4^5 + 8372*r2^12*r4^6 + 6656*r2^11*r4^7 + 858*r2^10*r4^8 - 5720*r2^9*r4^9 - 8866*r2^8*r4^10 - 7488*r2^7*r4^11 - 4004*r2^6*r4^12 - 1232*r2^5*r4^13 - 60*r2^4*r4^14 + 128*r2^3*r4^15 + 59*r2^2*r4^16 + 12*r2*r4^17 + r4^18)

res[2] |> println

   4*r4^(3/2)/sqrt(r2 + r4)

a はかなり複雑な分数式であるが，分子，分母を別々に簡約化して再度分数式にすることで，簡約化できる。

num = res[1] |> numerator |> factor
den = res[1] |> denominator |> factor;

ans_a = num/den
ans_a |> println

   4*r2*r4^(3/2)/((r2 - r4)*sqrt(r2 + r4))

二等辺三角形の底辺の長さの半分 a は 4*r2*r4^(3/2)/((r2 - r4)*sqrt(r2 + r4)) により計算できる。
r2 =7.5, r4 = 5 のとき，a = 37.9473319220205 である。

ans_a(r2 => 7.5, r4 => 5).evalf() |> println

   37.9473319220205

x3 は 4*r4^(3/2)/sqrt(r2 + r4) により計算できる。
r2 =7.5, r4 = 5 のとき，x3 = 12.6491106406735 である。

res[2](r2 => 7.5, r4 => 5).evalf() |> println

   12.6491106406735

a, x3 が求まった後，b, r1, r3 は逆にたどって以下のようにして求める。

r2 = 7.5
r4 = 5
a = 4*r2*r4^(3/2)/((r2 - r4)*sqrt(r2 + r4))
x3 = 4*r4^(3/2)/sqrt(r2 + r4)
b = 2*r4*x3^2/(-4*r4^2 + x3^2)
r1 = b/2 - r4
r3 = x3^2/(4*r4)
(a, b, r1, r3, x3) |> println

   (37.94733192202055, 26.666666666666657, 8.333333333333329, 8.000000000000002, 12.649110640673518)

乙円，丁円の直径が 15 寸，10 寸のとき，丙円の直径は 16 寸である。　
以下のような図になる。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r4) = (15, 10)./2
   r2 = 10
   r4 = 5
   a = 4*r2*r4^(3/2)/((r2 - r4)*sqrt(r2 + r4))
   x3 = 4*r4^(3/2)/sqrt(r2 + r4)
   b = 2*r4*x3^2/(-4*r4^2 + x3^2)
   r1 = b/2 - r4
   r3 = x3^2/(4*r4)
   @printf("乙円，丁円の直径が %g, %g のとき，丙円の直径は %g である。\n", 2r2, 2r4, 2r3)
   plot([a, 0, -a, a], [0, b, 0, 0], color=:blue, lw=0.5)
   circle(0, 2r4 + r1, r1)
   circle(0, 2r4 + r2, r2, :orange)
   circle2(x3, r3, r3, :green)
   circle(0, r4, r4, :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, 2r4 + r1, "甲円:r1,(0,2r4+r1)", :red, :center, delta=-delta)
       point(0, 2r4 + r2, "乙円:r2,(0,2r4+r2)", :orange, :center, delta=-delta)
       point(x3, r3, "丙円:r3\n(x3,r3)", :green, :center, delta=-delta)
       point(0, r4, "丁円:r4\n(0,r4)", :magenta, :center, delta=-delta)
       point(a, 0, "a", :blue, :left, :bottom, delta=delta/2)
       point(0, b, "b", :blue, :center, :bottom, delta=delta)
   end
end;

算額（その1095）

算額（その1095）

五十七 岩手県花泉町 花泉天満宮 文政13年(1830)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.http://www.wasan.jp/yamamura/yamamura.html
キーワード：円9個，外円，楕円2個，正方形

正方形の中に内接する大円，2 個の楕円，8 個の等円を容れる。正方形の一辺の長さが 1 寸のとき，楕円の短径はいかほどか。

正方形の一辺の長さを 2a
等円の半径と中心座標を r, (a - r, a - r), (a - r - √2r, a - r - √2r)
楕円の長半径と短半径，中心座標を 0, 0, (a, b)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive, r::positive, x0::positive, y0::positive
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (a - r -　√Sym(2)r - x0)/(a - r -　√Sym(2)r - y0)
eq2 = eq2 |> simplify |> numerator
eq3 = 2(a - r)^2 - (a + r)^2 |> expand
eq4 = (a - r -　√Sym(2)r - x0)^2 + (a - r -　√Sym(2)r - y0)^2 - r^2 |> expand;

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)
   (b, r, x0, y0) = u
   return [
       -1 + y0^2/b^2 + x0^2/a^2,  # eq1
       a^2*y0*(-a + r + sqrt(2)*r + x0) - b^2*x0*(-a + r + sqrt(2)*r + y0),  # eq2
       a^2 - 6*a*r + r^2,  # eq3
       2*a^2 - 4*sqrt(2)*a*r - 4*a*r - 2*a*x0 - 2*a*y0 + 5*r^2 + 4*sqrt(2)*r^2 + 2*r*x0 + 2*sqrt(2)*r*x0 + 2*r*y0 + 2*sqrt(2)*r*y0 + x0^2 + y0^2,  # eq4
   ]
end;

a = 1/2
iniv = BigFloat[0.25, 0.086, 0.267, 0.211]
res = nls(H, ini=iniv)

   ([0.24968850668279452, 0.08578643762690495, 0.2670793193663204, 0.2110827336905282], true)

正方形の一辺の長さが 1 寸のとき，楕円の短径は 0.49937701336558904 寸である。

「答」では 二分七厘五毛になっているが？

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1/2
   (b, r, x0, y0) = res[1]
   @printf("正方形の一辺の長さが %g のとき，楕円の短径は %g である。\n", 2a, 2b)
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:blue, lw=0.5)
   circle(0, 0, a, :magenta)
   circle4(a - r, a - r, r)
   circle4(a - r - √2r, a - r - √2r, r)
   ellipse(0, 0, a, b, color=:green)
   ellipse(0, 0, b, a, color=: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(x0, y0, "(x0,y0)", :green, :left, :bottom, delta=delta, deltax=-2delta)
       point(a - r, a - r, "等円:r\n(a-r,a-r)", :red, :center, :bottom, delta=delta/2)
       point(a - r - √2r, a - r - √2r, "等円:r\n(a-r-√2r,a-r2-√r)", :red, :center, :bottom, delta=delta/2)
       point(a - r - √2r, a - r - √2r)
       point(0, b, "b", :green, :center, :bottom, delta=delta/2)
       point(0, a, "a", :green, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
   end
end;