算額（その635）

算額（その635）

算額（その360）と同じ図形

埼玉県加須市騎西 玉敷神社 大正4年(1915)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

東京都中央区日本橋 福徳神社（芽吹稲荷）
https://mebuki.jp/2125/
玉敷神社に奉納された算額中の問題を復元したもの。

甲円 2 個と 2 本の直線と 2 個の乙円がある。甲円の直径を 3 とするとき，乙円の直径を求めよ。

算額（その628）を単純化したもの。
三角形の底辺の長さを 2a とする。
甲円の半径と中心座標を r1, (0, r1), (0, 3r1)
乙円の半径と中心座標を r2, (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms r1::positive, r2::positive, x2::positive, y2::positive, a::positive

sinθ = r1/3r1
#eq1 = (4r1 - (2r1 + r3))sinθ - r3
eq2 = a/sqrt((4r1)^2 + a^2) - sinθ
eq3 = r1*sinθ - (r1 - 2r2)
eq4 = x2^2 + (y2 - 3r1)^2 - (r1 - r2)^2
eq5 = (y2 - 3r1)/x2 *(-4r1)/a + 1;

res = solve([eq2, eq3, eq4, eq5], (a, r2, x2, y2))

   1-element Vector{NTuple{4, Sym}}:
    (sqrt(2)*r1, r1/3, 4*sqrt(2)*r1/9, 29*r1/9)

乙円の直径は甲円の直径の 1/3 倍である。
甲円の直径が 3 寸のとき，乙円の直径は 1 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 3//2
   (a, r2, x2, y2) = r1 .* (√2, 1/3, 4√2/9, 29/9)
   @printf("乙円の直径 = %g\n", 2r2)
   @printf("a = %g;  r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", a, r1, r2, x2, y2)
   plot([a, 0, -a, 0], [0, 4r1, 0, 0], color=:orange, lw=0.5)
   circle(0, r1, r1)
   circle(0, 3r1, r1)
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(0, r1, " 甲円:r1,(0,r1)", :red, :left, :vcenter)
       point(0, 3r1, " 3r1", :red, :left, :vcenter)
       point(x2, y2, "乙円:r2\n(x2,y2)", :blue, :center, :top, delta=-delta/2)
       point(a, 0, "a", :black, :left, :bottom, delta=delta)
       point(0, 4r1, " 4r1", :black, :left, :bottom)
   end
end;

算額（その634）

算額（その634）

埼玉県加須市馬内 秋葉神社

東京都中央区日本橋 福徳神社（芽吹稲荷）
https://mebuki.jp/2125/

埼玉県加須市馬内の秋葉神社に奉納された算額中の問題を復元したもの。

大円の中に上円，中円，下円が円に接する 2 本の直線に挟まれている。上円，下円の直径を 16，64 とするとき，外円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
上円の半径と中心座標を r1, (r1, y1)
中円の半径と中心座標を r2, (0, y2)
下円の半径と中心座標を r3, (r3, y3)
2直線の交点は y 軸上にあり，その座標を (0, a) とする。
描かれる図形の正しさを保証する座標として，直線と外円との交点座標を (x, sqrt(R^2 - x^2)) とする。
以下の連立方程式を解く。
直線と y 軸のなす角を 2θ，(0, a) と下円の中心座標 (r3, y3) がなす角は θ で，tan(2θ) と tan(θ) の関係式が連立方程式 eq6 である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms a::positive, R::positive, r1::positive, y1::positive,
     r2::positive, y2::positive, r3::positive, y3::negative,
     x::positive
eq1 = r1^2 + y1^2 - (R - r1)^2
eq2 = r1^2 + (y1 - y2)^2 - (r1 + r2)^2
eq3 = r3^2 + y3^2 - (R - r3)^2
eq4 = r3^2 + (y3 - y2)^2 - (r2 + r3)^2
eq5 = r1/(a - y1) - r3/(a - y3)
tanθ = r1/(a - y1)
tan2θ = r2/sqrt((a - y2)^2 - r2^2)
eq6 = tan2θ - 2tanθ/(1 - tanθ^2)  # tan の倍角公式
eq7 = x/(a + sqrt(R^2 - x^2)) - tan2θ;  # 外円との交点

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 v, r.f_converged
end;

function H(u)
   (a, R, y1, r2, y2, y3, x) = u
   return [
       r1^2 + y1^2 - (R - r1)^2,  # eq1
       r1^2 - (r1 + r2)^2 + (y1 - y2)^2,  # eq2
       r3^2 + y3^2 - (R - r3)^2,  # eq3
       r3^2 - (r2 + r3)^2 + (-y2 + y3)^2,  # eq4
       r1/(a - y1) - r3/(a - y3),  # eq5
       -2*r1/((a - y1)*(-r1^2/(a - y1)^2 + 1)) + r2/sqrt(-r2^2 + (a - y2)^2),  # eq6
       -r2/sqrt(-r2^2 + (a - y2)^2) + x/(a + sqrt(R^2 - x^2)),  # eq7
   ]
end;

(r1, r3) = (16, 64).//2
iniv = BigFloat[100, 75, 70, 35, 25, -30, 70]
res = nls(H, ini=iniv)

   (BigFloat[101.0160878276326834706245414713648871161063810216537006427555673980074763401048, 76.49999999999999999999999999999999999999999999999999999999762230488260311415535, 68.03124282269139907205326262357227091492878721866473716757506133442061889385097, 33.99999999999999999999999999999999999999999999999999999999738290330058695649773, 26.80018656651479357383916406383150066345679496492853282359138444046472151565239, -30.9232921921324541236605739198055776886039941903021532579865373591499745691432, 69.04228072652930803351703141806364055524282066404273049867895237707748909712568], true)

外円の直径 = 153;  a = 101.0160878;  R = 76.5;  y1 = 68.0312;  r2 = 34;  y2 = 26.8002;  y3 = -30.9233

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3) = (16, 64) .// 2
   (a, R, y1, r2, y2, y3) = (230, 160, 140, 75, 50, -70) .// 2
   (a, R, y1, r2, y2, y3, x) = res[1]
   @printf("外円の直径 = %g;  a = %.10g;  R = %g;  y1 = %g;  r2 = %g;  y2 = %g;  y3 = %g\n", 2R, a, R, y1, r2, y2, y3)
   plot()
   circle(0, 0, R)
   circle(r1, y1, r1, :blue)
   circle(-r1, y1, r1, :blue)
   circle(0, y2, r2, :green)
   circle(r3, y3, r3, :magenta)
   circle(-r3, y3, r3, :magenta)
   segment(0, a, x, -sqrt(R^2 - x^2), :orange)
   segment(0, a, -x, -sqrt(R^2 - x^2), :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:gray80, lw=0.5)
       hline!([0], color=:gray80, lw=0.5)
       point(0, a, " a", :black, :left, :vcenter)
       point(x, -sqrt(R^2-x^2), "(x,-√(R^2-x^2)) ", :black, :right, delta=-delta)
       point(r1, y1, " 上円:r1,(r1,y1)", :blue, :left, :vcenter)
       point(0, y2, " 中円:r2\n (0,y2)", :green, :left, :vcenter)
       point(r3, y3, "下円:r3,(r3,y3)", :magenta, :center, :bottom, delta=delta/2)
   end
end;

算額（その633）

算額（その633）

一〇五 加須市騎西町中種足 雷神社 大正元年(1912)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

東京都中央区日本橋 福徳神社（芽吹稲荷）
https://mebuki.jp/2125/

正方形の中に対角線と小さい正方形，正三角形と大小の円がある。小円の直径を 1 とするとき，大円の直径を求めよ。

正方形の一辺の長さを a とする。x 軸にある正三角形の頂点の座標を (x, 0) とする。
大円，小円の半径を r1, r2 として，まず x を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms l::positive, a::positive, x::positive, r1::positive, r2::positive
eq1 = (x - a/2)^2 +(a/2)^2 - 2(a - x)^2
res = solve(eq1, x);
res |> println
res[1](a => 1).evalf(),res[2](a => 1).evalf()

   Sym[a*(3 - sqrt(3))/2, a*(sqrt(3) + 3)/2]

   (0.633974596215561, 2.36602540378444)

解が 2 個求められるが，最初のものが適解である。
x = a*(3 - √3)/2

x = res[1]
x |> println

   a*(3 - sqrt(3))/2

正三角形の一辺の長さ l は √2(a - x) である。

l = sqrt(Sym(2))*(a - x)
l |> println

   sqrt(2)*(-a*(3 - sqrt(3))/2 + a)

大円，小円は（二等辺）直角三角形に内接しているので，「鈎股弦中の円の直径」公式により求められる。

r1 = (2(a/2) - (a/2)*sqrt(Sym(2)))/2
r2 = (2(a - x) - l)/2
r1 |> println
r2 |> println

   -sqrt(2)*a/4 + a/2
   -a*(3 - sqrt(3))/2 + a - sqrt(2)*(-a*(3 - sqrt(3))/2 + a)/2

r2 を 1 としたときの r1 は，(-sqrt(2)*a/4 + a/2)/(-a*(3 - sqrt(3))/2 + a - sqrt(2)*(-a*(3 - sqrt(3))/2 + a)/2) のようになるが，これは簡約化すると 1/2 + sqrt(3)/2 = (1 + √3)/2 = 1.36602540378444 になる。

ratio = r1/r2
ratio |> println
apart(ratio) |> simplify |> println
ratio.evalf() |> println

   (-sqrt(2)*a/4 + a/2)/(-a*(3 - sqrt(3))/2 + a - sqrt(2)*(-a*(3 - sqrt(3))/2 + a)/2)
   1/2 + sqrt(3)/2
   1.36602540378444

小円の直径を 1 とするとき，大円の直径は約 1.366 である。

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

a = 1;  x = 0.633975;  l = 0.517638;  r1 = 0.146447;  r2 = 0.107206;  r1/r2 = 1.36603

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1
   x = a*(3 - √3)/2
   l = √2*(a - x)
   r1 = (2(a/2) - (a/2)√2)/2
   r2 = (2(a - x) - l)/2
   @printf("a = %.10g;  x = %g;  l = %g;  r1 = %g;  r2 = %g;  r1/r2 = %g\n", a, x, l, r1, r2, r1/r2)
   plot()
   rect(0, 0, a, a, :black)
   rect(0, a/2, a/2, a, :blue)
   segment(0, 0, a, a, :black)
   plot!([a/2, x, a, a/2], [a/2, 0, a - x, a/2], color=:orange, lw=0.5)
   circle(r1, a/2 - r1, r1,)
   circle(a/2 + r1, a - r1, r1,)
   circle(a - r2, r2, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       # vline!([0], color=:gray80, lw=0.5)
       # hline!([0], color=:gray80, lw=0.5)
       point(0, a/2, " a/2", :black, :left, :bottom, delta=delta/2)
       point(0, a, " a", :black, :left, :bottom, delta=delta/2)
       point(a/2, a/2, "(a/2,a/2) ", :black, :right, :bottom, delta=delta/2)
       point(x, 0, "x ", :black, :right, :bottom, delta=delta/2)
       point(a, 0, "a ", :black, :right, :bottom, delta=delta/2)
       point(a/2, a, "(a/2,a)", :black, :center, :bottom, delta=delta/2)
       point(a, a - x, "(a,a-x)  ", :black, :right, :vcenter)
       point(r1, a/2 - r1, "大円:r1\n(r1,a/2-r1)", :red, :center, delta=-delta/2)
       point(a/2 + r1, a - r1, "大円:r1\n(a/2+r1,a-r1)", :red, :center, delta=-delta/2)
       point(a - r2, r2, "小円:r2\n(a-r2,r2)", :green, :center, delta=-delta/2)
   end
end;

算額（その629）

算額（その629）

和算図形問題あれこれ
令和4年5月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html

正方形内に正六角形を入れる。正方形の一辺の長さが 1219 寸のとき，正六角形の一辺の長さはいかほどか。

正方形の一辺の長さを a とする。b，c を図のように定義し，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms l::positive, a::positive, b::positive, c::positive

eq0 = 2(a/2 - c)^2 - l^2
eq1 = (a/2)^2 + (b - a/2)^2 - l^2
eq2 = c^2 + (b - c)^2 - l^2
res = solve([eq0, eq1, eq2], (l, b, c))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (a*sqrt(2 - sqrt(3)), -sqrt(3)*a/2 + 3*a/2, a*(1 - sqrt(3)/2))
    (a*sqrt(sqrt(3) + 2), sqrt(3)*a/2 + 3*a/2, a*(sqrt(3)/2 + 1))

最初のものが適解である。

正六角形の一辺の長さは a*sqrt(2 - sqrt(3)) = a√(2 - √3) である。

正方形の一辺の長さ a が 1219 寸のとき，正六角形の一辺の長さは a√(2 - √3) = 631.0008319599457 寸である。

a = 1219
a√(2 - √3)

   631.0008319599457

その他のパラメータは以下の通り。

a = 1219;  b = 772.815;  c = 163.315

function draw(more=false)
   pyplot(size=(300, 300), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1219
   (l, b, c) = (a*sqrt(2 - sqrt(3)), -sqrt(3)*a/2 + 3*a/2, a*(1 - sqrt(3)/2))
   @printf("正六角形の一辺の長さ = %g;  a = %g;  b = %g;  c = %g\n", l, a, b, c)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:cyan, lw=0.5, seriestype=:shape, alpha=0.5)
   plot!([c, b, a, a - c, a - b, 0, c], [c, 0, a - b, a - c, a, b, c], color=:orange, lw=0.5, seriestype=:shape, alpha=0.8)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(a, 0, "a ", :black, :right, :bottom, delta=delta/2)
       point(c, c, "(c,c)", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :vcenter)
   end
end;