算額（その1380）

算額（その1380）

一七 大里郡岡部村岡 稲荷社 文化13年(1816)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：直角三角形，中鈎，正方形
#Julia, #SymPy, #算額, #和算

問：只云大積四寸小積一寸問方積

非常に短く書かれているが，もう少し分かりやすく書き直す。
直角三角形の中に正方形を容れる。直角三角形の直角を挟む二辺の短い方を鈎，長い方を股，斜辺を弦と呼ぶ。直角の頂点から弦へおろした垂線を中鈎と呼ぶ。中鈎と弦の交点座標により区分される弦の短い方を短弦，長い方を長と弦呼ぶ。短弦と中鈎で構成される直角三角形の面積を小積，長弦と中鈎で構成される直角三角形の面積を大積と呼ぶ。大積が 4 平方寸(注1），小積が 1 平方寸のとき，正方形の一辺の長さはいかほどか。

注1：当時は，単位の次元の概念がなかった。「寸」×「寸」は「平方寸」であるが，「寸」と表していた。
注2：大積，小積の定義は明記されていないが，図に「中鈎」が描かれているので，短弦，長弦を底辺，中鈎を高さとする直角三角形の面積と解釈した。

正方形の一辺の長さを a
鈎，股，弦，短弦，長弦，中鈎に対する変数をそのまま，「鈎」，「股」，「弦」，「短弦」，「長弦」，「中鈎」
大積，小積を K1, K2
とおき，以下の連立方程式を解く。

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

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive, 短弦::positive, 長弦::positive,
     中鈎::positive, x0::positive, y0::positive, a::positive, K1::positive, K2::positive
eq1 = 短弦/中鈎 - 鈎/股
eq2 = 中鈎/長弦 - 鈎/股
eq3 = 鈎^2 + 股^2 - 弦^2
eq4 = 短弦 + 長弦 - 弦
eq5 = 長弦*中鈎/2 - K1  # 大積
eq6 = 短弦*中鈎/2 - K2  # 小積
eq7 = (鈎 - a)/a - a/(股 - a)  # 正方形により区分される 2 つの直角三角形は相似であること
eq8 = sqrt(x0^2 + y0^2) - 中鈎
eq9 = (鈎 - y0)/x0 - y0/(股 - x0)
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9],
   (a, 鈎, 股, 弦, 短弦, 長弦, 中鈎, x0, y0))[1]

   (K1^(1/4)*K2^(1/4)*(sqrt(K1) - sqrt(K2))*sqrt(2*K1 + 2*K2)/(K1 - K2), K2^(1/4)*sqrt(2*K1 + 2*K2)/K1^(1/4), K1^(1/4)*sqrt(2*K1 + 2*K2)/K2^(1/4), sqrt(2)*(K1 + K2)/(K1^(1/4)*K2^(1/4)), sqrt(2)*K2^(3/4)/K1^(1/4), sqrt(2)*K1^(3/4)/K2^(1/4), sqrt(2)*K1^(1/4)*K2^(1/4), sqrt(2)*K1^(1/4)*K2^(3/4)/sqrt(K1 + K2), sqrt(2)*K1^(3/4)*K2^(1/4)/sqrt(K1 + K2))

for i = 1:9
   println("$i: $(res[i](K1 => 4, K2 => 1).evalf())")
end

   1: 1.49071198499986
   2: 2.23606797749979
   3: 4.47213595499958
   4: 5.00000000000000
   5: 1.00000000000000
   6: 4.00000000000000
   7: 2.00000000000000
   8: 0.894427190999916
   9: 1.78885438199983

正方形の一辺の長さ

res[1] |> println
res[1](K1 => 4, K2 => 1) |> println
res[1](K1 => 4, K2 => 1).evalf() |> println

   K1^(1/4)*K2^(1/4)*(sqrt(K1) - sqrt(K2))*sqrt(2*K1 + 2*K2)/(K1 - K2)
   2*sqrt(5)/3
   1.49071198499986

正方形の面積

S = res[1]^2
S |>  println
S(K1 => 4, K2 => 1) |> println
S(K1 => 4, K2 => 1).evalf() |> println

   sqrt(K1)*sqrt(K2)*(sqrt(K1) - sqrt(K2))^2*(2*K1 + 2*K2)/(K1 - K2)^2
   20/9
   2.22222222222222

短弦 = 1, 中鈎 = 2, 小積 = 1
長弦 = 4, 中鈎 = 2, 大積 = 4
正方形の一辺の長さ = 1.49071, 面積 = 2.22222

全てのパラメータは以下のとおりである。
K1 = 4;  K2 = 1;  a = 1.49071;  鈎 = 2.23607;  股 = 4.47214;  弦 = 5;  短弦 = 1;  長弦 = 4;  中鈎 = 2;  x0 = 0.894427;  y0 = 1.78885

function draw(K1, K2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, 鈎, 股, 弦, 短弦, 長弦, 中鈎, x0, y0) = (K1^(1/4)*K2^(1/4)*(sqrt(K1) - sqrt(K2))*sqrt(2*K1 + 2*K2)/(K1 - K2), K2^(1/4)*sqrt(2*K1 + 2*K2)/K1^(1/4), K1^(1/4)*sqrt(2*K1 + 2*K2)/K2^(1/4), sqrt(2)*(K1 + K2)/(K1^(1/4)*K2^(1/4)), sqrt(2)*K2^(3/4)/K1^(1/4), sqrt(2)*K1^(3/4)/K2^(1/4), sqrt(2)*K1^(1/4)*K2^(1/4), sqrt(2)*K1^(1/4)*K2^(3/4)/sqrt(K1 + K2), sqrt(2)*K1^(3/4)*K2^(1/4)/sqrt(K1 + K2))
   @printf("短弦 = %g, 中鈎 = %g, 小積 = %g\n", 短弦, 中鈎, 短弦*中鈎/2)
   @printf("長弦 = %g, 中鈎 = %g, 大積 = %g\n", 長弦, 中鈎, 長弦*中鈎/2)
   @printf("正方形の一辺の長さ = %g, 面積 = %g\n", a, a^2)
   @printf("K1 = %g;  K2 = %g;  a = %g;  鈎 = %g;  股 = %g;  弦 = %g;  短弦 = %g;  長弦 = %g;  中鈎 = %g;  x0 = %g;  y0 = %g\n", K1, K2, a, 鈎, 股, 弦, 短弦, 長弦, 中鈎, x0, y0)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   rect(0, 0, a, a, :red)
   segment(0, 0, x0, y0, :green)
   if more == true
       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(a, a, "(a,a)", :red, :left, :bottom, delta=delta/2)
       point(x0, y0, "(x0,y0)", :green, :left, :bottom, delta=delta/2)
       point(股, 0, "股", :blue, :left, :bottom, delta=delta)
       point(0, 鈎, " 鈎", :blue, :left, :bottom, delta=delta)
       point(x0/2, y0/2, "中鈎", :green, mark=false)
   end
end;

draw(4, 1, true)

 

算額（その1379）

算額（その1379）

七十三 群馬県安中市下後閑 威徳神社 嘉永3年(1850)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円5個，楕円，正方形
#Julia, #SymPy, #算額, #和算

第七問
「群馬の算額」によれば，「十問あるが文章は不明」とのことであるが，以下のようなものでもあろう。

正三角形の中に大円 1 個，小円 4 個，斜線 3 本を容れる。大円は正三角形の三辺と 3 箇所で接し，正三角形の隅にある小円は正三角形の二辺と大円に接する。また，斜線はそれぞれ正三角形の隅にある小円 2 個と中央にある小円に接する。
正三角形の一辺の長さが与えられたとき，小円の直径はいかほどか。
または，小円の直径が与えられたとき，正三角形の一辺の長さはいかほどか。

正三角形の一辺の長さを 2a
大円の半径と中心座標を r1, (0, 0); r1 = a/√3
小円の半径と中心座標を r2, (0, r1 + r2), (x2, r2); r2 = (a - x2)/√3
とおき，以下の方程式を解き，x2 を求める。

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

using SymPy
@syms a, r1, r2, x2
r2 = (a - x2)/√Sym(3)
eq1 = x2^2 + (r2 - a/√Sym(3))^2 - (a/√Sym(3) + r2)^2
solve(eq1, x2)[2] |> println

   2*a/3

x2 = 2a/3 であり，小円の半径は r2 = (a - x2)/√3 = √3a/9 である。

正三角形の一辺の長さが 1 のとき，小円の直径は √3/9 = 0.19245008972987523 である。
小円の直径が 1 のとき，正三角形の一辺の長さは 5.196152422706632 である。

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = a/√3
   x2 = 2a/3
   r2 = √3a/9  # (a - x2)/√3
   @printf("正三角形の一辺の長さが %g のとき，小円の直径は %.15g である。\n", 2a, 2r2)
   plot([a, 0, -a, a], [-a/√3, 2a/√3, -a/√3, -a/√3], color=:blue, lw=0.5)
   circle(0, 0, r1)
   circle(0, 0, r2, :green)
   rotate(0, r1 + r2, r2, :green)
   segment((2a/√3 - r1 - r2)/√3, r1 + r2, -a + 2(2a/√3 - r1 - r2)/√3, -a/√3, :orange)
   segment(-(2a/√3 - r1 - r2)/√3, r1 + r2, a - 2(2a/√3 - r1 - r2)/√3, -a/√3, :orange)
   segment(-a + 2r2/√3, 2r2 - a/√3, a - 2r2/√3, 2r2 - a/√3, :orange)
   if more == true
       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 + r2, " r1+r2", :green, :left, :vcenter)
       point(0, 2a/√3, " 2a/√3", :green, :left, :vcenter)
       point(a, -a/√3, "(a,-a/√3)", :green, :right, delta=-delta/2)
       point(x2, r2 - a/√3, "(x2,r2-a/√3)", :green, :center, delta=-delta/2, deltax=2delta)
       ylims!(-a/√3 - 5delta, 2a/√3 + 5delta)
   end
end;

draw(1/2, true)

 

算額（その1378）

算額（その1378）

七十三 群馬県安中市下後閑 威徳神社 嘉永3年(1850)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円5個，楕円，正方形
#Julia, #SymPy, #算額, #和算

第六問
「群馬の算額」によれば，「十問あるが文章は不明」とのことであるが，以下のようなものでもあろう。

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

または，等円の直径が与えられたとき，正方形の一辺の長さを求めよ。

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

しかし，SymPy の能力的に解析解を求めることができないので，a = 1/2 とおき，数値解を求める（得られる図形は相似なので，一般性を損なわない） 。

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

using SymPy
@syms a, b, r, x0, y0
b = r
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0*(a - r - y0) + (a - r - x0)*(a^2*y0)
eq3 = (x0 - (a - r))^2 + (y0 - (a - r))^2 - r^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)
   (r, x0, y0) = u
   return [
       -1 + y0^2/r^2 + x0^2/a^2,  # eq1
       a^2*y0*(a - r - x0) - r^2*x0*(a - r - y0),  # eq2
       -r^2 + (-a + r + x0)^2 + (-a + r + y0)^2,  # eq3
   ]
end;

a = 1/2
iniv = BigFloat[0.18, 0.28, 0.15]
res = nls(H, ini=iniv)

   ([0.1785016026524881, 0.27969014973090045, 0.14796196722191884], true)

正方形の一辺の長さが 1 のとき，等円の直径は 0.357003205304976 である。
（等円の直径が 1 のとき，正方形の一辺の長さは 2.8010952986983217 である。）

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

draw(true)

 

算額（その1377）

算額（その1377）

七十二 群馬県富岡市一ノ宮 貫前神社 嘉永2年(1849)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円6個，長方形
#Julia, #SymPy, #算額, #和算

長方形の中に6個の等円を容れる。長方形の短辺が 5 寸のとき，長辺はいかほどか。

長方形の長辺と短辺を 2a, 2b とする。
等円 A を中心とする（補助円を含む）周りの6円の中心は正六角形を構成する。
等円の半径と中心座標を r, A:(r*cos(30°), -r*sin(30°))), B:(r*cos(30°), b - r), C:(a - r, r -  b)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms a::positive, b::positive, r::positive
eq1 = 2r + (r - r*sind(Sym(30))) - b
eq2 = 3r*cosd(Sym(30)) + r - a
res = solve([eq1, eq2], (a, r));

長辺は短辺の (2 + 3√3)/5 倍である。
短辺が 5 寸のとき，長辺は 2 + 3√3 = 7.19615242270663 寸である。

res[a] |> factor |> println
2*res[a](b => 5/2).evalf() |> println

   b*(2 + 3*sqrt(3))/5
   7.19615242270663

ちなみに，等円の直径は短辺の 2/5 倍である。
短辺が 5 寸のとき，等円の直径は 2 寸である。

2res[r](b => 5/2).evalf() |> println

   2.00000000000000

function draw(b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r) = (b*(2 + 3√3)/5, 2b/5)
   @printf("長方形の短辺が %g のとき，長辺は %g，等円の直径は %g である。\n", 2b, 2a, 2r)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:blue, lw=0.5)
   circle(r - a, b - r, r)
   circle(a - r, r - b, r)
   circle(r*cosd(30), -r*sind(30), r)
   circle(-r*cosd(30), r*sind(30), r)
   circle(a - r, r*sind(30), r, :gray70)
   circle(r*cosd(30), b - r, r)
   circle(r*cosd(30), -b, r, :gray70)
   circle(-r*cosd(30), r - b, r)
   if more == true
       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*cosd(30), -r*sind(30), " A")
       point(-r*cosd(30), r*sind(30), " D")
       point(r*cosd(30), b - r, " B")
       point(a - r, r*sind(30), " D'")
       point(a - r, r - b, " C")
       point(-r*cosd(30), r - b, " C'")
       point(r*cosd(30), -b, " B'", :green, :left, :bottom)
       point(a, b, "(a,b)", :blue, :right, :bottom)
       point(0, 0, " O", :black, :left)
   end
end;

draw(5/2, true)

算額（その1376）

算額（その1376）

七十二 群馬県富岡市一ノ宮 貫前神社 嘉永2年(1849)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：楕円2個，正方形3個，直線上
#Julia, #SymPy, #算額, #和算

直線上に合同な 2 個の楕円が載っており，その隙間に 3 個の正方形を容れる。楕円の長径が 2 寸のとき，短径はいかほどか。

縦長の楕円であるが慣例により 短半径を a，長半径を b とする（a < b）。
中心座標は (a, b)
正方形の一辺の長さを c
とおき，以下の連立方程式を解く。

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

using SymPy
@syms a::positive, b::positive, c::positive
eq1 = (c/2 - a)^2/a^2 + (2c - b)^2/b^2 - 1
eq2 = (c - a)^2/a^2 + (c - b)^2/b^2 - 1
res = solve([eq1, eq2], (a, c))[1]  #1 of 4

   (b/2, b/5)

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

function draw(b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, c) = (b/2, b/5)
   @printf("楕円の長径が %g のとき，短径は %g，正方形の一辺の長さは %g である。\n", 2b, 2a, c)
   plot()
   ellipse(a, b, a, b, color=:red)
   ellipse(-a, b, a, b, color=:red)
   rect(-c/2, c, c/2, 2c, :blue)
   rect(0, 0, c, c, :blue)
   rect(0, 0, -c, c, :blue)
   if more == true
       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(a, b, "楕円:a,b,(a,b)", :red, :center, delta=-delta)
       point(c/2, 2c, " (c/2,2c)", :blue, :left, :vcenter)
       point(c, c, " (c,c)", :blue, :left, :vcenter)
   end
end;

draw(2/2, true)

 

算額（その1375）

算額（その1375）

七十二 群馬県富岡市一ノ宮 貫前神社 嘉永2年(1849)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円1個，外円，正五角形
#Julia, #SymPy, #算額, #和算

外円の中に 3 個の合同な正五角形を入れる。外円の直径が 10 寸のとき，正五角形の一辺の長さはいかほどか。

問題文では，「円の中に 3 個の正五角形を容れる」と書いているが，3 個ではなく位置と方向を明示すれば「1 個を容れる」だけで十分である。もちろん図のように「10 個を容れる」としてもよいが，それではヒントを与えすぎると考えたのだろうか。

ここでは，青で描いた上部の正五角形について考える。赤で描いた外円の半径と中心座標を R, (0, 0) とする。正五角形が内接する灰色で描いた円の半径を AB = BD = r とする。
与えられるのは OA = R，求めるのは x = CD である。
∠CBD = 36°，∠COD = 18° である。
BC = r*cos(36°)
CD = r*sin(36°) より，r = CD/sin(36°) = x/sin(36°)
OC*tan(18°) = CD
OA = AB + BC + OC
R = r + r*cos(36°) + x/tan(18°)
R = x/sin(36°)*(1 + cos(36°)) + x/tan(18°)
これを x について解いて
x = 2R*sqrt(35 - 15√5)/(√10sqrt(5 - 2√5) + 5√2*sqrt(5 - 2√5) + 2√5*sqrt(5 - √5))
を得る。

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

using SymPy
@syms R, r, x
eq = x/sind(Sym(36))*(1 + cosd(Sym(36))) + x/tand(Sym(18)) - R;
ans_x = solve(eq, x)[1]
ans_x |> println

   2*R*sqrt(35 - 15*sqrt(5))/(sqrt(10)*sqrt(5 - 2*sqrt(5)) + 5*sqrt(2)*sqrt(5 - 2*sqrt(5)) + 2*sqrt(5)*sqrt(5 - sqrt(5)))

SymPy では変数が残っている式は簡約化されにくいという特徴（欠点）を持っているので，R の係数部分（下記）だけを簡約化する。

ans_x/R |> println

   2*sqrt(35 - 15*sqrt(5))/(sqrt(10)*sqrt(5 - 2*sqrt(5)) + 5*sqrt(2)*sqrt(5 - 2*sqrt(5)) + 2*sqrt(5)*sqrt(5 - sqrt(5)))

@syms d
apart(ans_x/R, d) |> simplify |> println

   sqrt(25 - 10*sqrt(5))/10

長い式は sqrt(25 - 10√5)/10 に簡約化される。
術では sqrt(0.25 - sqrt(0.05)) としているが 1/10 を外側の平方根の中に入れると同じである。

正五角形の一辺の長さ 2x は 外円の半径が R のとき　2R*sqrt(25 - 10√5)/10 である。
外円の直径が 10 寸のとき，正五角形の一辺の長さは sqrt(25 - 10√5) = 1.6245984811645313 寸である。

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   # x = 2R*sqrt(35 - 15√5)/(√10sqrt(5 - 2√5) + 5√2*sqrt(5 - 2√5) + 2√5*sqrt(5 - √5))
   x = R*sqrt(25 - 10√5)/10
   r = x/sind(36)
   θ = 18:72:378
   xa = @. r*cosd(θ)
   ya = @. r*sind(θ) + (R - r)
   plot([-r*sind(36), 0, r*sind(36), 0], [R - r*(1 + cosd(36)), 0, R - r*(1 + cosd(36)), R - r], color=:gray70, lw=0.5)
   circle(0, 0, R, :red)
   circle(0, R - r, r, :gray70)
   for i = 1:10
       plot!(xa, ya, color=:green, lw=1)
       (xa, ya) = transform(xa, ya, deg=36)
   end
   plot!(xa, ya, color=:blue, lw=1)
   if more == true
       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, R, " A", :blue, :left, :bottom, delta=delta/2)
       point(0, R - r, " B", :blue, :left, :bottom, delta=delta/2)
       point(0, R - r*(1 + cosd(36)), " C", :blue, :left, :bottom, delta=delta/2)
       point(r*sind(36), R - r*(1 + cosd(36)), " D", :blue, :left, :bottom, delta=delta/2)
       point(0, 0, " O", :blue, :left, :bottom, delta=delta/2)
   end
end;

draw(10/2, true)

算額（その1374）

算額（その1374）

神壁算法 天明5年乙巳 秋九月 第二術 關流四傅藤田権平定資六弟子 東都青山 西村治郎右衛門政央
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
キーワード：球5個，外球，3次元，ソディ・ゴセットの n+1 球定理
#Julia, #SymPy, #算額, #和算

大球の中に，甲球，乙球，丙球，丁球の 4 球を容れる。甲球，乙球，丙球，丁球の直径が 3 寸，2 寸，1.2寸， 1 寸のとき，大球の直径はいかほどか。

注：甲球，乙球，丙球，丁球は互いに外接し，それぞれは大球に内接している。

この問題を解くには，2 次元平面上で 4 個の円の半径 rᵢ から計算される曲率 kᵢ = 1/rᵢ の関係についてのデカルトの円定理を3次元に拡張したものを使う。

一般的には，「ソディ・ゴセットの n+1 球定理」https://ikuro-kotaro.sakura.ne.jp/koramu/8023_p2.htm というもので，n 次元空間では，互いに接する n + 2 個の球の間に

　　　n*(k₁² + k₂² + k₃² + ... + kₙ₊₂²) = (k₁ + k₂ + k₃ + ... ± kₙ₊₂)²

が成り立つ。

n = 2 の場合がデカルトの円定理である。

本問は n = 3 の場合であり，以下の式が成り立つ。
互いに外接する 4 個の球のとき，k₅² の符号がマイナスの場合は 4 個の球は 5 番目の球に内接し，プラスの場合は 5 番目の球に外接する。

　　　3(k₁² + k₂² + k₃² + k₄² + kₙ₊₂²) = (k₁ + k₂ + k₃ + k₄ ± k₅)²

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

4 個の球が 5 番目の球に内接する場合 --- 　5 番目の球の半径

@syms r1, r2, r3, r4, r5, k1, k2, k3, k4, k5
eq01 = 3(k1^2 + k2^2 + k3^2 + k4^2 + k5^2) - (k1 + k2 + k3 + k4 - k5)^2;
K5 = solve(eq01, k5)[2]  # 2 of 2
K5 |> println

   -k1/2 - k2/2 - k3/2 - k4/2 + sqrt(3)*sqrt(-k1^2 + 2*k1*k2 + 2*k1*k3 + 2*k1*k4 - k2^2 + 2*k2*k3 + 2*k2*k4 - k3^2 + 2*k3*k4 - k4^2)/2

1/K5(k1 => 1/r1, k2 => 1/r2, k3 => 1/r3, k4 => 1/r4)(r1 => 3//2, r2 => 1, r3 => 12//20, r4 => 1//2) |> println

   3

4 個の球が 5 番目の球に外接する場合 --- 5 番目の球の半径

eq02 = 3(k1^2 + k2^2 + k3^2 + k4^2 + k5^2) - (k1 + k2 + k3 + k4 + k5)^2;
K6 = solve(eq02, k5)[2]
K6 |> println

   k1/2 + k2/2 + k3/2 + k4/2 + sqrt(3)*sqrt(-k1^2 + 2*k1*k2 + 2*k1*k3 + 2*k1*k4 - k2^2 + 2*k2*k3 + 2*k2*k4 - k3^2 + 2*k3*k4 - k4^2)/2

1/K6(k1 => 1/r1, k2 => 1/r2, k3 => 1/r3, k4 => 1/r4)(r1 => 3//2, r2 => 1, r3 => 12//20, r4 => 1//2) |> println

   3/17

以下の関数は，4個の球の半径 r1, r2, r3, r4 を与えると，4 球が内接する球と4 球が外接する球の両方の半径を返す。

function radius_of_fifth_sphere(r1, r2, r3, r4)
   (k1, k2, k3, k4) = 1 ./ (r1, r2, r3, r4)
   #return 1/(-k1/2 - k2/2 - k3/2 - k4/2 + sqrt(3)*sqrt(-k1^2 + 2*k1*(k2 + k3 + k4) - k2^2 + 2*k2*(k3 + k4) - k3^2 + 2*k3*k4 - k4^2)/2),
   #       1/( k1/2 + k2/2 + k3/2 + k4/2 + sqrt(3)*sqrt(-k1^2 + 2*k1*(k2 + k3 + k4) - k2^2 + 2*k2*(k3 + k4) - k3^2 + 2*k3*k4 - k4^2)/2)
   (a, b) = (k1/2 + k2/2 + k3/2 + k4/2, sqrt(3)*sqrt(-k1^2 + 2*k1*(k2 + k3 + k4) - k2^2 + 2*k2*(k3 + k4) - k3^2 + 2*k3*k4 - k4^2)/2)
   return 1/(b - a), 1/(b + a)
end;

問では，甲球，乙球，丙球，丁球の直径が 3 寸，2 寸，1.2寸，1 寸のときなので，4 球が内接する大球の半径は 3 寸と計算される。

radius_of_fifth_sphere(3/2, 2/2, 1.2/2, 1/2)

   (3.0000000000000027, 0.17647058823529413)

算額の問は，大球の直径だけを問うているので，「大球の直径は 6 寸」で完了である。

以下では，図を描くために必要なパラメータを求める。

座標軸は任意に設定できるので，大球の中心を原点とし，甲球の中心は x 軸上にあり，乙球は x−y 平面上にあるように設定する。後に明らかになるが，乙球の中心は y 軸上にある。
大球の半径と中心座標を R, (0, 0, 0)
甲球の半径と中心座標を r1, (R - r1, 0, 0)
乙球の半径と中心座標を r2, (x2, y2, 0)
丙球の半径と中心座標を r3, (x3, y3, z3)
丁球の半径と中心座標を r4, (x4, y4, z4)
とおき，以下の連立方程式を（逐次的に）解く。

using SymPy
@syms R::positive, r1::positive,
     r2::positive, x2::positive, y2::positive,
     r3::positive, x3::positive, y3::positive, z3::positive,
     r4::positive, x4::positive, y4::positive, z4::positive;

eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = x3^2 + y3^2 + z3^2 - (R - r3)^2
eq3 = x4^2 + y4^2 + z4^2 - (R - r4)^2
eq4 = (R - r1)^2 + y2^2 - (r1 + r2)^2
eq5 = (R - r1 - x3)^2 + y3^2 + z3^2 - (r1 + r3)^2
eq6 = (R - r1 - x4)^2 + y4^2 + z4^2 - (r1 + r4)^2
eq7 = (x2 - x3)^2 + (y2 - y3)^2 + z3^2 - (r2 + r3)^2
eq8 = (x2 - x4)^2 + (y2 - y4)^2 + z4^2 - (r2 + r4)^2
eq9 = (x3 - x4)^2 + (y3 - y4)^2 + (z3 - z4)^2 - (r3 + r4)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (R, x2, y2, x3, y3, z3, x4, y4, z4))

1.　R を求める

「ソディ・ゴセットの n+1 球定理」から，R = 3 なので，甲球の半径と中心座標は 1.5, (1.5, 0, 0) である。

2.　x2, y2 を求める

eq1, eq4 を解いて，x2, y2 を求める。
乙球の半径と中心座標は 1, (0, 2, 0) である。

res = solve([eq1, eq4], (x2, y2))[1]

   (sqrt(2)*sqrt(R^2 - R*r1 - R*r2 - r1*r2), sqrt(R + r2)*sqrt(-R + 2*r1 + r2))

x2
res[1](R => 3, r1 => 3//2, r2 => 2//2) |> println
y2
res[2](R => 3, r1 => 3//2, r2 => 2//2) |> println

   0
   2

3.　x3, y3, z3 を求める

eq2, eq5, eq7 を解いて，x3, y3, z3 を求める。
丙球の半径と中心座標は 0.6, (6/5, 9/5, 3√3/5) である。

res2 = solve([eq2, eq5, eq7], (x3, y3, z3))[1]  # 1 of 2

   ((-R^2 + R*r1 + R*r3 + r1*r3)/(-R + r1), (-R^3 + R^2*r1 + 2*R^2*r3 + 2*R^2*x2 - 2*R*r1*r3 - 2*R*r1*x2 + R*r2^2 + 2*R*r2*r3 - 2*R*r3*x2 - R*x2^2 - R*y2^2 - r1*r2^2 - 2*r1*r2*r3 - 2*r1*r3*x2 + r1*x2^2 + r1*y2^2)/(2*y2*(-R + r1)), -sqrt(-R^6/4 + R^5*r1/2 + R^5*r3 + R^5*x2 - R^4*r1^2/4 - 2*R^4*r1*r3 - 2*R^4*r1*x2 + R^4*r2^2/2 + R^4*r2*r3 - R^4*r3^2 - 3*R^4*r3*x2 - 3*R^4*x2^2/2 - R^4*y2^2/2 + R^3*r1^2*r3 + R^3*r1^2*x2 - R^3*r1*r2^2 - 2*R^3*r1*r2*r3 + 2*R^3*r1*r3^2 + 4*R^3*r1*r3*x2 + 3*R^3*r1*x2^2 + R^3*r1*y2^2 - R^3*r2^2*r3 - R^3*r2^2*x2 - 2*R^3*r2*r3^2 - 2*R^3*r2*r3*x2 + 2*R^3*r3^2*x2 + 3*R^3*r3*x2^2 + R^3*r3*y2^2 + R^3*x2^3 + R^3*x2*y2^2 + R^2*r1^2*r2^2/2 + R^2*r1^2*r2*r3 - R^2*r1^2*r3^2 - R^2*r1^2*r3*x2 - 3*R^2*r1^2*x2^2/2 - R^2*r1^2*y2^2/2 + 2*R^2*r1*r2^2*r3 + 2*R^2*r1*r2^2*x2 + 4*R^2*r1*r2*r3^2 + 4*R^2*r1*r2*r3*x2 - 2*R^2*r1*r3*x2^2 + 2*R^2*r1*r3*y2^2 - 2*R^2*r1*x2^3 - 2*R^2*r1*x2*y2^2 - R^2*r2^4/4 - R^2*r2^3*r3 - R^2*r2^2*r3^2 + R^2*r2^2*r3*x2 + R^2*r2^2*x2^2/2 + R^2*r2^2*y2^2/2 + 2*R^2*r2*r3^2*x2 + R^2*r2*r3*x2^2 + R^2*r2*r3*y2^2 - R^2*r3^2*x2^2 - R^2*r3*x2^3 - R^2*r3*x2*y2^2 - R^2*x2^4/4 - R^2*x2^2*y2^2/2 - R^2*y2^4/4 - R*r1^2*r2^2*r3 - R*r1^2*r2^2*x2 - 2*R*r1^2*r2*r3^2 - 2*R*r1^2*r2*r3*x2 - 2*R*r1^2*r3^2*x2 - R*r1^2*r3*x2^2 - 3*R*r1^2*r3*y2^2 + R*r1^2*x2^3 + R*r1^2*x2*y2^2 + R*r1*r2^4/2 + 2*R*r1*r2^3*r3 + 2*R*r1*r2^2*r3^2 - R*r1*r2^2*x2^2 - R*r1*r2^2*y2^2 - 2*R*r1*r2*r3*x2^2 - 2*R*r1*r2*r3*y2^2 - 2*R*r1*r3^2*x2^2 - 4*R*r1*r3^2*y2^2 + R*r1*x2^4/2 + R*r1*x2^2*y2^2 + R*r1*y2^4/2 - r1^2*r2^4/4 - r1^2*r2^3*r3 - r1^2*r2^2*r3^2 - r1^2*r2^2*r3*x2 + r1^2*r2^2*x2^2/2 + r1^2*r2^2*y2^2/2 - 2*r1^2*r2*r3^2*x2 + r1^2*r2*r3*x2^2 + r1^2*r2*r3*y2^2 - r1^2*r3^2*x2^2 + r1^2*r3*x2^3 + r1^2*r3*x2*y2^2 - r1^2*x2^4/4 - r1^2*x2^2*y2^2/2 - r1^2*y2^4/4)/(y2*(-R + r1)))

x3
res2[1](R => 3, r1 => 3//2, r2 => 2//2, r3 => 12//20) |> println
y3
res2[2](R => 3, r1 => 3//2, r2 => 2//2, r3 => 12//20, x2 => 0, y2 => 2) |> println
z3
res2[3](R => 3, r1 => 3//2, r2 => 2//2, r3 => 12//20, x2 => 0, y2 => 2) |> println

   6/5
   9/5
   3*sqrt(3)/5

4.　x4,y4,z4 を求める

eq3, eq6, eq8 を解いて，x4, y4, z4 を求める。
丁球の半径と中心座標は 0.5, (3/2, 2, 0) である。

res3 = solve([eq3, eq6, eq8], (x4, y4, z4))[1]  # 1 of 2

   ((-R^2 + R*r1 + R*r4 + r1*r4)/(-R + r1), (-R^3 + R^2*r1 + 2*R^2*r4 + 2*R^2*x2 - 2*R*r1*r4 - 2*R*r1*x2 + R*r2^2 + 2*R*r2*r4 - 2*R*r4*x2 - R*x2^2 - R*y2^2 - r1*r2^2 - 2*r1*r2*r4 - 2*r1*r4*x2 + r1*x2^2 + r1*y2^2)/(2*y2*(-R + r1)), -sqrt(-R^6/4 + R^5*r1/2 + R^5*r4 + R^5*x2 - R^4*r1^2/4 - 2*R^4*r1*r4 - 2*R^4*r1*x2 + R^4*r2^2/2 + R^4*r2*r4 - R^4*r4^2 - 3*R^4*r4*x2 - 3*R^4*x2^2/2 - R^4*y2^2/2 + R^3*r1^2*r4 + R^3*r1^2*x2 - R^3*r1*r2^2 - 2*R^3*r1*r2*r4 + 2*R^3*r1*r4^2 + 4*R^3*r1*r4*x2 + 3*R^3*r1*x2^2 + R^3*r1*y2^2 - R^3*r2^2*r4 - R^3*r2^2*x2 - 2*R^3*r2*r4^2 - 2*R^3*r2*r4*x2 + 2*R^3*r4^2*x2 + 3*R^3*r4*x2^2 + R^3*r4*y2^2 + R^3*x2^3 + R^3*x2*y2^2 + R^2*r1^2*r2^2/2 + R^2*r1^2*r2*r4 - R^2*r1^2*r4^2 - R^2*r1^2*r4*x2 - 3*R^2*r1^2*x2^2/2 - R^2*r1^2*y2^2/2 + 2*R^2*r1*r2^2*r4 + 2*R^2*r1*r2^2*x2 + 4*R^2*r1*r2*r4^2 + 4*R^2*r1*r2*r4*x2 - 2*R^2*r1*r4*x2^2 + 2*R^2*r1*r4*y2^2 - 2*R^2*r1*x2^3 - 2*R^2*r1*x2*y2^2 - R^2*r2^4/4 - R^2*r2^3*r4 - R^2*r2^2*r4^2 + R^2*r2^2*r4*x2 + R^2*r2^2*x2^2/2 + R^2*r2^2*y2^2/2 + 2*R^2*r2*r4^2*x2 + R^2*r2*r4*x2^2 + R^2*r2*r4*y2^2 - R^2*r4^2*x2^2 - R^2*r4*x2^3 - R^2*r4*x2*y2^2 - R^2*x2^4/4 - R^2*x2^2*y2^2/2 - R^2*y2^4/4 - R*r1^2*r2^2*r4 - R*r1^2*r2^2*x2 - 2*R*r1^2*r2*r4^2 - 2*R*r1^2*r2*r4*x2 - 2*R*r1^2*r4^2*x2 - R*r1^2*r4*x2^2 - 3*R*r1^2*r4*y2^2 + R*r1^2*x2^3 + R*r1^2*x2*y2^2 + R*r1*r2^4/2 + 2*R*r1*r2^3*r4 + 2*R*r1*r2^2*r4^2 - R*r1*r2^2*x2^2 - R*r1*r2^2*y2^2 - 2*R*r1*r2*r4*x2^2 - 2*R*r1*r2*r4*y2^2 - 2*R*r1*r4^2*x2^2 - 4*R*r1*r4^2*y2^2 + R*r1*x2^4/2 + R*r1*x2^2*y2^2 + R*r1*y2^4/2 - r1^2*r2^4/4 - r1^2*r2^3*r4 - r1^2*r2^2*r4^2 - r1^2*r2^2*r4*x2 + r1^2*r2^2*x2^2/2 + r1^2*r2^2*y2^2/2 - 2*r1^2*r2*r4^2*x2 + r1^2*r2*r4*x2^2 + r1^2*r2*r4*y2^2 - r1^2*r4^2*x2^2 + r1^2*r4*x2^3 + r1^2*r4*x2*y2^2 - r1^2*x2^4/4 - r1^2*x2^2*y2^2/2 - r1^2*y2^4/4)/(y2*(-R + r1)))

x4
res3[1](R => 3, r1 => 3//2, r2 => 2//2, r3 => 12//20, r4 => 1//2) |> println
y4
res3[2](R => 3, r1 => 3//2, r2 => 2//2, r3 => 12//20, r4 => 1//2, x2 => 0, y2 => 2) |> println
z4
res3[3](R => 3, r1 => 3//2, r2 => 2//2, r3 => 12//20, r4 => 1//2, x2 => 0, y2 => 2) |> println

   3/2
   2
   0

以下は，図を描くためのプログラムである。3次元の図は GeoGebra によって描いた。

function radius_of_fifth_sphere(r1, r2, r3, r4)
   (k1, k2, k3, k4) = 1 ./ (r1, r2, r3, r4)
   (a, b) = (k1/2 + k2/2 + k3/2 + k4/2, sqrt(3)*sqrt(-k1^2 + 2*k1*(k2 + k3 + k4) - k2^2 + 2*k2*(k3 + k4) - k3^2 + 2*k3*k4 - k4^2)/2)
   return 1/(b - a), 1/(b + a)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, r4) = (3/2, 2/2, 1.2/2, 1/2)
   (R, junk) = radius_of_fifth_sphere(r1, r2, r3, r4)
   (x2, y2) = (0, 2)
   (x3, y3, z3) = (6/5, 9/5, 3√3/5)
   (x4, y4, z4) = (3/2, 2, 0)
   plot()
   circle(0, 0, R, :brown)
   circle(R - r1, 0, r1, :green)
   circle(x2, y2, r2, :orange)
   circle(x3, y3, r3, :blue)
   circle(x4, y4, r4, :magenta)
   if more == true
       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, -1.57, "z 軸の正の方向から原点方向を見た透視図\n甲球，乙球，丁球の中心は x-y 平面にあり,\nx-y 平面による切断面は外接し合う円である。\n丙球の中心はそれらの上方にあり,\n甲球，乙球，丁球と外接している。", :black, :center, mark=false)
       point(R - r1, 0, "甲球:r1,(R-r1,0,0)", :green, :center, delta=-delta/2)
       point(0, 2, "乙球:r2,(0,2,0)", :orange, :right, delta=-delta/2, deltax=7delta)
       point(6/5, 9/5, "丙球:r3,(6/5,9/5,3√3/5)", :blue, :left, delta=-delta/2, deltax=-6delta)
       point(3/2, 2, "丁球:r4,(3/2,2,0)", :magenta, :left, :bottom, delta=delta/2, deltax=-6delta)
   end
end;

draw(true)

 

算額（その1373）

算額（その1373）

四十 群馬県前橋市山王町 日枝神社 文政7年(1824)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：楕円，等脚台形，最大面積
#Julia, #SymPy, #算額, #和算

斜辺と上底が等しい（等斜と呼ぶ）等脚台形の中に楕円を容れる。等斜が 10 寸のとき，等脚台形の面積が最大になるときの楕円の長径はいかほどか。

等脚台形の，上底，下底，高さを 2j, 2k, h
楕円の長半径，短半径，中心座標を a, b, (0, b); b = h/2
楕円と等脚台形の斜辺の接点座標を (x0, y0)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms j::positive, k::positive, h::positive,
     S::positive, a::positive, b::positive,
     x0::positive, y0::positive;
b = h/2
eq1 = (k - j)^2 + h^2 - (2j)^2
eq2 = x0^2/a^2 + (y0 - b)^2/b^2 - 1
eq3 = -b^2*x0/(a^2*(y0 - b)) + h/(k - j)
eq4 = y0/(k - x0) - h/(k - j);

k, a, x0, y0 は h, j を含む式になる。

res = solve([eq1, eq2, eq3, eq4], (k, a, x0, y0))[2]  # 1 of 2

   (j + sqrt(-h^2 + 4*j^2), sqrt(j)*sqrt(j + sqrt(-h^2 + 4*j^2)), 2*j*(h^2 - 2*j^2 + j*sqrt(-h^2 + 4*j^2))/h^2, (h^2 - 2*j^2 + j*sqrt(-h^2 + 4*j^2))/h)

等脚台形の面積は h*(4*j + 2*sqrt(-h^2 + 4*j^2))/2 で，等斜(=2j)が与えられれば，h だけの関数 h*(4j + 2sqrt(4*j^2 - h^2))/2　になる。

S = (2j + 2res[1])*h/2
S |> println

   h*(4*j + 2*sqrt(-h^2 + 4*j^2))/2

等斜 = 10（j = 5）のとき，面積は下図のようになり，h が 9 前後のときに最大値をとる。
注：高さは 10 より大きくなれない。高さが 10 は極限状態であり，そのとき等脚台形は正方形になる。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(S(j => 10/2), xlims=(0, 10), xlabel="h", ylabel="S")

面積 S が最大になるときの高さ h を求めるには，面積を表す式を h について微分し，その式が 0 になるときの h を求めればよい。

diff(S, h)

   -h^2/sqrt(-h^2 + 4*j^2) + 2*j + sqrt(-h^2 + 4*j^2)

面積 S が最大になるときの h は √3j である。等斜が 10 寸のときは，h = 8.66025403784439 である。

h0 = solve(diff(S, h))[1]
h0[h] |> println

   sqrt(3)*j

h = √3j のときの面積は 3*√3*j^2/2 である。

S(h => sqrt(Sym(3))*j) |> println

   3*sqrt(3)*j^2

等斜が 10 寸のとき，面積は 3√3*5^2= 75√3/2 = 129.9038105676658 である。

要求されているのは楕円の長径である。
h = √3j のとき，長径は 2√2*j = √2*等斜= 14.1421356237310 である。

2res[2](h => j*√Sym(3)) |> simplify |> println

   2*sqrt(2)*j

なお，h = √3j のとき，下底は 20 である。

2res[1](h => √Sym(3)j)(j => 10//2) |> println

   20

面積は確かに (10 + 20) * 8.66025403784439 / 2 = 129.90381056766586 になる。

function draw(等斜, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   j = 等斜/2
   h = √3j
   (k, a, x0, y0) = (sqrt(100 - h^2) + 5, sqrt(5*sqrt(100 - h^2) + 25), 10*(h^2 + 5*sqrt(100 - h^2) - 50)/h^2, (h^2 + 5*sqrt(100 - h^2) - 50)/h)
   S = h*(sqrt(100 - h^2) + 10)/2
   b = h/2
   @printf("等斜（上底）が %g，下底が %g のとき，楕円の長径は %g である。\n", 等斜, 2k, 2a)
   plot([k, j, -j, -k, k], [0, h, h, 0, 0], color=:red, lw=0.5)
   ellipse(0, b, a, b, color=:blue)
   if more == true
       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(j, h, "(j,h)", :red, :right, :bottom, delta=delta)
       point(k, 0, "k", :red, :left, :bottom, delta=delta, deltax=delta/2)
       point(0, b, " b", :blue, :left, :vcenter)
       point(0, 2b, " 2b", :blue, :left, :bottom, delta=delta)
       point(x0, y0, "(x0,y0)", :blue, :left, :vcenter, deltax=2delta)
   end
end;

draw(10, true)

算額（その1372）

算額（その1372）

七十二 群馬県富岡市一ノ宮 貫前神社 嘉永2年(1849)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：内接円，外円，個数
#Julia, #SymPy, #算額, #和算

大円の中に互いに外接し合い大円に内接する等円を描く。大円と等円の直径が与えられたとき，等円の個数を求める術を述べよ。

大円と等円の半径を R, r，等円の個数を n とすれば，以下の関係式が成り立つ。

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

using SymPy
@syms R::positive, r::positive, n::positive
eq1 = asind(r/(R - r)) - 180/n
eq1 |> println

   180*asin(r/(R - r))/pi - 180/n

1.　等円の半径を求める

大円の半径と等円の個数を指定して等円の半径を求めるには，この関係式を等円の半径について解けばよい。

res = solve(eq1, r)[1]
res |> println

   R*sin(pi/n)/(sin(pi/n) + 1)

たとえば，外円の半径が 1，等円の個数が 8 のとき，等円の半径は以下により求めることができる（直径を与えれば直径が求まる）。

res(R => 1/2, n => 8).evalf() |> println
res(R => 1, n => 8).evalf() |> println

   0.138384326957078
   0.276768653914155

sin を計算するのはコンピュータにとっては何の苦もないことであるが，例えば電卓で計算するときには，近似式を使えばよい。
近似式 sin(x) ≈ x - x^3/6 + x^5/120 は小数点以下 5，6 桁まで正しい。

sin0(x) = x - x^3/6 + x^5/120

x = pi/n
(R*sin0(x)/(sin0(x) + 1))(R => 1/2, n => 8).evalf() |> println

   0.138384401530346

2.　等円の個数を求める

大円と等円の半径（直径）を指定して等円の個数を求めるには，前述の関係式を等円の個数について解く（必要ならば四捨五入する）。

res2 = solve(eq1, n)[1]
res2 |> println

   pi/asin(r/(R - r))

たとえば，外円，等円の半径が 0.5，0.138384326957078 のとき，等円の個数は 8 個である。

res2(R => 0.5, r => 0.138384326957078).evalf() |> println

   7.99999999999997

asin(arcsin) を計算するのはコンピュータにとっては何の苦もないことであるが，例えば電卓で計算するときには，近似式を使えばよい。
近似式 arcsin(x) ≈ x + x^3/6 + 3x^5/40 は小数点以下 4，5 桁まで正しい。

arcsin0(x) = x + x^3/6 + 3x^5/40

R = 1/2
r = 0.138384326957078
x = r/(R - r)
asin(x) |>  println
arcsin0(x) |> println

   0.39269908169872586
   0.392639425546951

等円の個数は，8.001215489793262 もしくは 7.997159214528561 になり，四捨五入すれば正しい答えになる。

pi/arcsin0(x) |> println
3.14/arcsin0(x) |> println

   8.001215489793262
   7.997159214528561

なお，「術」では，円周率*(大円の直径 - 等円の直径)/等円の直径 としている。これは上記近似式の最初の項だけを使用している。つまり arcsin(x) ≈ x としている。

大円の直径 = 1
等円の直径 = 0.276768653914155
円周率 = 3.14
円周率*(大円の直径 - 等円の直径)/等円の直径 |> println
円周率/(等円の直径/(大円の直径 - 等円の直径)) |> println

   8.205215419423654
   8.205215419423654

正確な asin を使おうが近似式を使おうが大円の直径，等円の直径が正確でなければ，計算結果は整数値からかなり外れることになるが，多くの場合は四捨五入すれば答えが得られるであろう。

function draw(R, r, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   n = pi/asin(r/(R - r))
   n0 = 3.14/(r/(R - r))
   @printf("外円の直径が %g，等円の直径が %g のとき，等円の個数は %d である(近似値は %d)。\n", 2R, 2r, n, n0)
   plot()
   circle(0, 0, R)
   rotate(R - r, 0, r, angle=360/n, :blue)
   if more == true
       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)
   end
end;

R = 1/2
n = 13
r = R*sin(pi/n)/(sin(pi/n) + 1)
draw(R, r, true)

 

算額（その1371）

算額（その1371）

六十二 群馬県高崎市石原町 清水寺 天保10年(1837)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：三角形，二等辺三角形3個
#Julia, #SymPy, #算額, #和算

三角形の中に斜線を 2 本引き，3 個の二等辺三角形を作る。それぞれの二等辺三角形の斜辺（×をつけた）の長さは小斜に等しい。中斜，小斜が与えられたとき，大斜を求める術を述べよ。

3 個の二等辺三角形の頂点の座標を，(x1, y1), (x2, y2), (x3, 0), (x4, 0)
中斜，小斜の長さを 「中斜」，「小斜」
とおき，以下の連立方程式を解く。
大斜 = x4 である。

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

using SymPy
@syms x1::positive, y1::positive,
     x2::positive, y2::positive,
     x3::positive, x4::positive,
     小斜::positive, 中斜::positive;
eq1 = x2^2 + y2^2 - 小斜^2
eq2 = (x3 - x2)^2+y2^2 - 小斜^2
eq3 = (x1 - x3)^2 + y1^2 - 小斜^2
eq3 = (x4 - x1)^2 + y1^2 - 小斜^2
eq4 = 小斜 + sqrt((x1 - x2)^2 + (y1 - y2)^2) - 中斜
eq5 = (x3 + x4)/2 - x1
eq6 = y2/x2 - y1/x1
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (x1, y1, x2, y2, x3, x4))[2]  # 2 of 2

   (中斜*(中斜^3 - 3*中斜^2*小斜 + 3*中斜*小斜^2 - 小斜^3)*sqrt(中斜^3 - 中斜^2*小斜 - 中斜*小斜^2 + 小斜^3)/(2*sqrt(小斜)*(中斜^4 - 4*中斜^3*小斜 + 6*中斜^2*小斜^2 - 4*中斜*小斜^3 + 小斜^4)), 中斜*sqrt(-中斜 + 3*小斜)/(2*sqrt(小斜)), sqrt(小斜)*(中斜^3 - 3*中斜^2*小斜 + 3*中斜*小斜^2 - 小斜^3)*sqrt(中斜^3 - 中斜^2*小斜 - 中斜*小斜^2 + 小斜^3)/(2*(中斜^4 - 4*中斜^3*小斜 + 6*中斜^2*小斜^2 - 4*中斜*小斜^3 + 小斜^4)), sqrt(小斜)*sqrt(-中斜 + 3*小斜)/2, sqrt(小斜)*(中斜^3 - 3*中斜^2*小斜 + 3*中斜*小斜^2 - 小斜^3)*sqrt(中斜^3 - 中斜^2*小斜 - 中斜*小斜^2 + 小斜^3)/(中斜^4 - 4*中斜^3*小斜 + 6*中斜^2*小斜^2 - 4*中斜*小斜^3 + 小斜^4), sqrt(中斜^3 - 中斜^2*小斜 - 中斜*小斜^2 + 小斜^3)/sqrt(小斜))

res[6] |> factor |> println

   sqrt(中斜 + 小斜)*Abs(中斜 - 小斜)/sqrt(小斜)

大斜（x4）は，「中斜と小斜の和を小斜で割ったもの平方根に中斜と小斜の差を掛ける」と得られる。

たとえば，中斜 = 1.5，小斜 = 0.6 のとき，大斜は 1.5459002853702086 である。

中斜 = 1.5
小斜 = 0.65
sqrt((中斜 + 小斜)/小斜)*(中斜 - 小斜)

   1.5459002853702086

function draw(中斜, 小斜, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, y1, x2, y2, x3, x4) = (中斜*(中斜^3 - 3*中斜^2*小斜 + 3*中斜*小斜^2 - 小斜^3)*sqrt(中斜^3 - 中斜^2*小斜 - 中斜*小斜^2 + 小斜^3)/(2*sqrt(小斜)*(中斜^4 - 4*中斜^3*小斜 + 6*中斜^2*小斜^2 - 4*中斜*小斜^3 + 小斜^4)), 中斜*sqrt(-中斜 + 3*小斜)/(2*sqrt(小斜)), sqrt(小斜)*(中斜^3 - 3*中斜^2*小斜 + 3*中斜*小斜^2 - 小斜^3)*sqrt(中斜^3 - 中斜^2*小斜 - 中斜*小斜^2 + 小斜^3)/(2*(中斜^4 - 4*中斜^3*小斜 + 6*中斜^2*小斜^2 - 4*中斜*小斜^3 + 小斜^4)), sqrt(小斜)*sqrt(-中斜 + 3*小斜)/2, sqrt(小斜)*(中斜^3 - 3*中斜^2*小斜 + 3*中斜*小斜^2 - 小斜^3)*sqrt(中斜^3 - 中斜^2*小斜 - 中斜*小斜^2 + 小斜^3)/(中斜^4 - 4*中斜^3*小斜 + 6*中斜^2*小斜^2 - 4*中斜*小斜^3 + 小斜^4), sqrt(中斜^3 - 中斜^2*小斜 - 中斜*小斜^2 + 小斜^3)/sqrt(小斜))
   @printf("中斜が %g，小斜が %g のとき，大斜は %g である。\n", 中斜, 小斜, x4)
   @printf("x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g;  x3 = %g;  x4 = %g\n", x1, y1, x2, y2, x3, x4)

   plot([0, x4, x1, 0], [0, 0, y1, 0], color=:magenta, lw=0.5)
   segment(x2, y2, x3, 0)
   segment(x3, 0, x1, y1)
   if more == true
       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(x2, y2, "(x2,y2) ", :black, :right, :bottom, delta=delta)
       point(x1, y1, "(x1,y1) ", :black, :right, :bottom, delta=delta)
       point(0, 0, "0", :black, :center, delta=-3delta)
       point(x3, 0, "x3", :black, :center, delta=-3delta)
       point(x4, 0, "x4", :black, :center, delta=-3delta)
       dimension_line(x1 + 2.5delta, y1, x4 + 2.5delta, 0, "小斜", :green, :left, deltax=2delta)
       dimension_line(-delta/2, 2delta, x1 - delta/2, y1 + 2delta, "中斜", :blue, :center, delta=6delta)
       dimension_line(0, -2delta, x4, -2delta, "大斜", :red, :center, delta=-4delta)
       point(x2/2, y2/2, "×", :black, :center, :vcenter, mark=false)
       point((x2+x3)/2, y2/2, "×", :black, :center, :vcenter, mark=false)
       point((x1+x3)/2, y1/2, "×", :black, :center, :vcenter, mark=false)
       point((x1+x4)/2, y1/2, "×", :black, :center, :vcenter, mark=false)
       segment(0, 0, x2, y2, :black)
       plot!(xlims=(-5delta, x4 + 5delta), ylims=(-10delta, y1 + 8delta))
   end
end;

draw(1.5, 0.65, true)

 

算額（その1370）

算額（その1370）

六十 群馬県佐波郡玉村町飯塚 光琳寺 天保10年(1837)
六十一 群馬県高崎市石原町 清水寺 天保10年(1837)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：累円
#Julia, #SymPy, #算額, #和算

等脚台形の中に半円，全円，次円，三円，四円，五円を容れる。上底（大頭）が 90 寸のとき，各円の直径はいかほどか。

上底，下底を 2b, 2a
半円の半径と中心座標を r0, (0, r0); r0 = 2r1 = b/2
全円の半径と中心座標を r1, (0, r1)
次円の半径と中心座標を r2, (x21, r2), (x22, r0 - r2)
三円の半径と中心座標を r3, (x3, r3)
四円の半径と中心座標を r4, (x4, r4)
五円の半径と中心座標を r5, (x5, r5)
とおき，以下の連立方程式を解く。
四円のパラメータまでは一度に求めることができるが，SymPy の性能上，五円のパラメータも一緒に求めることができないので，分けて解く。

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

using SymPy
@syms a::positive, r0::positive, r1::positive,
     r2::positive, x21::positive, x22::positive,
     r3::positive, x3::positive,
     r4::positive, x4::positive,
     r5::positive, x5::positive;
r0 = 2r1
eq1 = x21^2 + (r0 - r2)^2 - (r0 + r2)^2
eq2 = x22^2 + r2^2 - (r0 - r2)^2
eq3 = x22^2 + (r0 - r2 - r1)^2 - (r1 + r2)^2
eq4 = dist2(r0, r0, a, 0, x21, r2, r2)
eq5 = x3^2 + (r0 - r3)^2 - (r0 + r3)^2
eq6 = (x21 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq7 = x4^2 + (r0 - r4)^2 - (r0 + r4)^2
eq8 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (r2, x21, x22, r3, x3, r4, x4, a))[4]  # 4 of 6

   (r1/2, 2*r1, sqrt(2)*r1, 2*r1/9, 4*r1/3, r1/8, r1, r1*(sqrt(2)/2 + 2))

@syms r5::positive, x5::positive
(r2, x21, x22, r3, x3, r4, x4, a) = (r1/2, 2*r1, sqrt(Sym(2))*r1, 2*r1/9, 4*r1/3, r1/8, r1, r1*(sqrt(Sym(2))/2 + 2))
eq9 = x5^2 + (r0 - r5)^2 - (r0 + r5)^2
eq10 = (x4 - x5)^2 + (r4 - r5)^2 - (r4 + r5)^2
#res2 = solve([eq7, eq8], (r4, x4,))
res2 = solve([eq9, eq10], (r5, x5))[1]  # 1 of 2

   (2*r1/25, 4*r1/5)

次円（二円），三円，四円，五円の半径は 全円の半径の 1/2, 2/9, 1/8, 2/25 倍である。

n 番目の円の半径は，全円の半径の 2/n^2 である。

上底が 90 のとき，次円の直径 = 22.5，三円の直径 = 10，四円の直径 = 5.625，五円の直径 = 3.6 である。

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

   r2 = 11.25;  x21 = 45;  x22 = 31.8198;  r3 = 5;  x3 = 30;  r4 = 2.8125;  x4 = 22.5;  r5 = 1.8;  x5 = 18;  a = 60.9099

function draw(b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = b/2
   r1 = r0/2
   (r2, x21, x22, r3, x3, r4, x4, a) = (r1/2, 2*r1, sqrt(2)*r1, 2*r1/9, 4*r1/3, r1/8, r1, r1*(sqrt(2)/2 + 2))
   (r5, x5) = (2*r1/25, 4*r1/5)
   @printf("上底が %g のとき，次円の直径 = %g，三円の直径 = %g，四円の直径 = %g，五円の直径 = %g である。\n", b, 2r2, 2r3, 2r4, 2r5)
   @printf("r2 = %g;  x21 = %g;  x22 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  x4 = %g;  r5 = %g;  x5 = %g;  a = %g\n", r2, x21, x22, r3, x3, r4, x4, r5, x5, a)
   plot([a, r0, -r0, -a, a], [0, r0, r0, 0, 0], color=:magenta, lw=0.5)
   circle(0, r0, r0, beginangle=180, endangle=360)
   circle(0, r1, r1, :blue)
   circle2(x21, r2, r2, :green)
   circle2(x22, r0 - r2, r2, :green)
   circle2(x3, r3, r3, :orange)
   circle2(x4, r4, r4, :brown)
   circle2(x5, r5, r5, :tomato)
   if more == true
       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(b/2, b/2, "(b/2,b/2)", :magenta, :right, :bottom, delta=delta)
       point(a, 0, "a", :magenta, :center, delta=-4delta)
       point(0, r1, "全円:r1,(0,r1)", :blue, :center, delta=-delta)
       point(x21, r2, "次円:r2\n(x21,r2)", :green, :center, :bottom, delta=delta)
       point(x22, r0 - r2, "次円:r2\n(x22,r0-r2)", :green, :center, :bottom, delta=delta)
       point(x3, r3, "三円:r3,(x3,r3)", :orange, :left, delta=-18delta)
       point(x4, r4, "四円:r4,(x4,r4)", :brown, :left, delta=-22delta)
       point(x5, r5, "五円:r5,(x5,r5)", :tomato, :left, delta=-28delta)
   end
end;

draw(90, true)

算額（その1369）

算額（その1369）

五十九 群馬県佐波郡玉村町飯塚 光琳寺 天保8年(1837)
六十一 群馬県高崎市石原町 清水寺 天保10年(1839)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円22個，部分円（円欠），面積
#Julia, #SymPy, #算額, #和算

水平な弦で欠き取られた円の一部の中に等円を 22 個容れる。等円の直径が 1 寸のとき，黒積（図で黒塗りした面積の2倍）はいかほどか。

等円の半径を r，部分円の半径を R とする。

⊿ABC は一辺の長さが 8r，高さ（OD）が 2√3r の正三角形である。
A を中心として B, C を通る円（灰色で示す）は O を中心とする元の円と相似で，半径は AC = 8r である。
したがって，A, C は O を中心とする円の直径の上にある。

黒積は，半径 R の面積の 1/4 から，台形（OABD）の面積と等円の面積の (2 + 120/360 +  240/360) 倍を差し引いたものの 2 倍である。

等円の面積の 120/360 は A の等円から差し引くべき面積，240/360 は B の等円から差し引くべき面積である。

r = 1/2
R = 8r
S = 2*(pi*R^2/4 - (r + 5r)*4√3r/2 - pi*r^2*(2 + 120/360 +  240/360))

   10.028047402920391

黒積は 10.028047402920391 である。

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

function draw(r, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 8r
   plot()
   circlef(0, 0, R, :black, beginangle=0, endangle=180)
   circlef(0, 0, R, :white, beginangle=90, endangle=270)
   plot!([5r, r, -r, -5r, 5r], [0, 4√3r, 4√3r, 0, 0], seriestype=:shape, color=:white, lw=0.5)
   circle(0, 0, R, :blue)
   [circle2f(i*r, 0, r, :white) for i in (1,3,5,7)]
   [circle2f(i*r, 2√3r, r, :white) for i in (1,3)]
   circle2f(r, 4√3r, r, :white)
   [circle2f(i*r, √3r, r, :white) for i in (2, 4)]
   circle2f(2r, 3√3r, r, :white)
   circlef(0, 3√3r, r, :white)
   circlef(0, √3r, r, :white)

   [circle2(i*r, 0, r) for i in (1,3,5,7)]
   [circle2(i*r, 2√3r, r) for i in (1,3)]
   circle2(r, 4√3r, r)
   [circle2(i*r, √3r, r) for i in (2, 4)]
   circle2(2r, 3√3r, r)
   circle(0, 3√3r, r)
   circle(0, √3r, r)
   circle(5r, 0, R, :gray70)
   plot!([5r, r, 0, 0, 5r], [0, 4√3r, 4√3r, 0, 0], color=:green, lw=0.5)
   if more == true
       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", :blue, :left, :bottom, delta=delta/2, deltax=delta/2)
       point(5r, 0, "A", :red, :center, delta=-delta/2)
       point(r, 4√3r, "B", :red, :center, :bottom, delta=delta/2)
       point(-3r, 0, "C", :red, :center, delta=-delta/2)
       point(0, 0, "O", :red, :left, delta=-delta/2, deltax=delta)
       point(0, 4√3r, "D", :red, :right, :bottom, delta=delta/2, deltax=-delta)
   end
end;

draw(1/2, true)

算額（その1368）

算額（その1368）

五十九 群馬県佐波郡玉村町飯塚 光琳寺 天保8年(1837)
六十一 群馬県高崎市石原町 清水寺 天保10年(1839)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：直角三角形，長方形，最大面積
#Julia, #SymPy, #算額, #和算

直角三角形の中に長方形を容れる。斜辺（弦）が 64 寸，長方形の長辺が 27 寸とする。長方形の面積が最大になるのは底辺（股）はいかほどか。

基本的に，算額（その816）と同じである。

鈎，股を A, B，鈎，股上にある点を a（平）, b（長） とおくと，長方形の面積 S は S = a*b である。
以下の連立方程式を解く。

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

using SymPy

@syms A, B, a, b, S

b = 27
eq1 = (A - a)/b - A/B
eq2 = a/(B - b) - A/B
eq3 = a*b - S
eq1 = A^2 + B^2 - 64^2
res = solve([eq1, eq2, eq3], (a, A, S))[2]

   (sqrt(4096 - B^2)*(B - 27)/B, sqrt(4096 - B^2), 27*sqrt(-(B - 64)*(B + 64))*(B - 27)/B)

面積 S は，B の関数である。

S = res[3]
S |> println

   27*sqrt(-(B - 64)*(B + 64))*(B - 27)/B

面積は B が 50 前後で最大値を取る。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(S, xlims=(27, 64), xlabel="B", ylabel="S(B)")

面積が最大になるのは上図の曲線の接線の傾きが 0 になるときである。数値的に解くには，S(B) を B で微分し，S'(B) = 0 になるときの B の値を求める。

diff(S, B) |> println

   27*sqrt(-(B - 64)*(B + 64))*(B - 27)/((B - 64)*(B + 64)) + 27*sqrt(-(B - 64)*(B + 64))/B - 27*sqrt(-(B - 64)*(B + 64))*(B - 27)/B^2

solve(diff(S, B), B)[1] |> println

   48

面積が最大になるのは，B が 48 のときである。

S(B => 48) |> println
S(B => 48).evalf() |> println

   189*sqrt(7)
   500.046997791208

B が 48 のとき，面積は最大値 500.046997791208 になる。

function draw(more)
    pyplot(size=(300, 300), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   S = 500.046997791208
   B = 48
   A = sqrt(4096 - B^2)
   a = sqrt(4096 - B^2)*(B - 27)/B
   b = S/a
   plot([0, B, 0, 0], [0, 0, A, 0], showaxis=false)
   rect(0, 0, b, a, :red)
   if more == true
       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)
       #plot!(showaxis=false, xlims=(-0.1, 9.5))
       point(0, A, " A（鈎）", :blue, :left, :bottom, delta=delta/2)
       point(B, 0, " B（股）", :blue, :left, :bottom, delta=delta/2)
       point(0, a, " a（平）", :red, :left, :bottom, delta=delta/2)
       point(b, 0, " b（長）", :red, :left, :bottom, delta=delta/2)
   end
end;

draw(true)

算額（その1367）

算額（その1367）

四十一 群馬県高崎市下小鳥町 幸宮神社 文政7年(1824)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：直角三角形，鈎股弦
#Julia, #SymPy, #算額, #和算

直角三角形の面積が 756 歩，股寸と弦寸の和が 1 丈 4 尺 7 寸のとき，釣，股，弦はいかほどか。

直角三角形の三辺を「鈎」，「股」，「弦」
直角三角形の面積を K1
股と弦の和を K2
とおき，以下の連立方程式を解く。

using SymPy
@syms 鈎, 股, 弦, K1, K2
eq1 = 鈎^2 + 股^2 - 弦^2  # こちらで与えると解けない
eq1 = sqrt(鈎^2 + 股^2) - 弦
eq2 = 鈎*股/2 - K1
eq3 = 股 + 弦 - K2;
res = solve([eq1, eq2, eq3], (鈎, 股, 弦))[3]  # 3 of 3

   (6^(2/3)*(12*2^(1/3)*K2^2 + 18^(1/3)*(1 + sqrt(3)*I)^2*(18*K1*K2 + sqrt(3)*sqrt(K2^2*(108*K1^2 - K2^4)))^(2/3))/(36*(1 + sqrt(3)*I)*(18*K1*K2 + sqrt(3)*sqrt(K2^2*(108*K1^2 - K2^4)))^(1/3)), 12*K1*(1 + sqrt(3)*I)*(108*K1*K2 + 6*sqrt(3)*sqrt(K2^2*(108*K1^2 - K2^4)))^(1/3)/(12*2^(1/3)*K2^2 + 18^(1/3)*(1 + sqrt(3)*I)^2*(18*K1*K2 + sqrt(3)*sqrt(K2^2*(108*K1^2 - K2^4)))^(2/3)), (12*K1*(1 + sqrt(3)*I)*(108*K1*K2 + 6*sqrt(3)*sqrt(K2^2*(108*K1^2 - K2^4)))^(1/3) + K2*(-12*2^(1/3)*K2^2 - 18^(1/3)*(1 + sqrt(3)*I)^2*(18*K1*K2 + sqrt(3)*sqrt(K2^2*(108*K1^2 - K2^4)))^(2/3)))/(-12*2^(1/3)*K2^2 - 18^(1/3)*(1 + sqrt(3)*I)^2*(18*K1*K2 + sqrt(3)*sqrt(K2^2*(108*K1^2 - K2^4)))^(2/3)))

解は 3 組得られるが，3 番目のものが適解である。虚数解として得られるものがあるが，虚部は実質的に 0 である。

鈎，股，弦はそれぞれ，21 寸，72 寸，75 寸である。

res[1](K1 => 756, K2 => 147).evalf() |> println
res[2](K1 => 756, K2 => 147).evalf() |> println
res[3](K1 => 756, K2 => 147).evalf() |> println

   21.0000000000000
   72.0 - 8.67361737988404e-19*I
   75.0000000000000

 

算額（その1366）

算額（その1366）

三十六 群馬県多野郡新町 稲荷神社 文政3年(1820)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円8個，等脚台形
#Julia, #SymPy, #算額, #和算

等脚台形の中に，大円 5 個，小円 3 個を容れる。台形の上底が与えられたとき，小円の直径はいかほどか。

注：算額（その1364）の設問の不備を修正し，枝葉を取り除いたものである。

等脚台形の下底，上底，高さを 2a, 2b, h
大円の半径と中心座標を r1, (2r1, r1), (r1, h - r1)
小円の半径と中心座標を r2, (0, 2r1 + r2)
甲円の半径と中心座標を r3, (x3, r3)
乙円の半径と中心座標を r4, (0, h - r4)
とおき，以下の連立方程式を解く。
方程式は算額（その1364）と同じで，与える条件と求めるパラメータが違うだけである。
甲円と乙円のパラメータも求めておく（図には描かない）。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, h::positive,
     r1::positive, r2::positive,
     r3::positive, x3::positive, r4::positive;
eq1 = r1^2 + (h - 2r1)^2 - 4r1^2
eq2 = r1/(a - 2r1) - 1/√Sym(3)
eq3 = r1/(b - r1) - √Sym(3);
# res1 = solve([eq1, eq2, eq3], (r1, a, h));

eq4 = r1^2 + (h - r1 - (2r1 + r2))^2 - (r1 + r2)^2
eq5 = r1^2 + (r1 - r4)^2 - (r1 + r4)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, r4, a, h))

eq6 = r3/(a - x3) - 1/√Sym(3)
eq7 = x3 - 2r1 - 2sqrt(r1*r3)
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7],
   (r1, r2, r3, x3, r4, a, h))[1]  # 2 of 2

   (b*(3 - sqrt(3))/2, -989*b*sqrt(151316*sqrt(3) + 262087) - 3*b/2 + sqrt(3)*b/2 + 571*sqrt(3)*b*sqrt(151316*sqrt(3) + 262087), b*(3 - sqrt(3))/6, 2*b, b*(3 - sqrt(3))/8, b*(sqrt(3) + 3)/2, b*(-2967*sqrt(151316*sqrt(3) + 262087) - 2*sqrt(3) + 6 + 1713*sqrt(453948*sqrt(3) + 786261))/2)

パラメータによっては二重根号を外して簡約化できる。

r1
res[1] |> println

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

r2
res[2] |> sympy.sqrtdenest |> simplify |> println

   b*(-5 + 3*sqrt(3))/2

r3
res[3] |> sympy.sqrtdenest |> simplify |> println

   b*(3 - sqrt(3))/6

x3
res[4] |> sympy.sqrtdenest |> simplify |> println

   2*b

r4
res[5] |> println

   b*(3 - sqrt(3))/8

a
res[6] |> println

   b*(sqrt(3) + 3)/2

h
res[7] |> sympy.sqrtdenest |> simplify |> println

   b*(sqrt(3) + 3)/2

小円の半径は b*(3 - sqrt(3))/6 なので，上底が 14 寸のとき小円の直径は 14*(3 - sqrt(3))/6 = 2.95854811567262 寸である。
このとき，下底は 32*(3 - sqrt(3))/3 = 13.524791385931977 である。

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

   a = 16;  r1 = 4.28719;  r2 = 0.66323;  r3 = 1.42906;  x3 = 13.5248;  r4 = 1.0718;  b = 6.7624;  h = 16

---

なお，小円の直径を知るだけならば，以下のような手順で小円の直径を（紙と鉛筆だけで）求めることができる。

上底 2b と大円の半径 r1 の関係式は以下であり，大円の半径は r1 = b*(3 - √3)/2 である。

@syms b, r1
eq01 = (1 + 1/√Sym(3))r1 - b
eq01 |> println

   -b + r1*(sqrt(3)/3 + 1)

r1 = solve(eq01, r1)[1] |> simplify
r1 |> println

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

一辺の長さが 2r1 の正三角形の頂点を中心とする3個の大円が作る中央の隙間にピッタリハマる円の半径は以下のようになる（図中の灰色で描いた小さな直角三角形を参照）。

r1*2/√Sym(3) - r1 |> simplify |> println

    b*(-5 + 3*sqrt(3))/2

上底が 14 寸のとき，小円の直径は 1.3730669589464242 寸である。

b = 14/2
2*b*(-5 + 3*sqrt(3))/2

   1.3730669589464242

function draw(b, more)
   pyplot(size=(800, 800), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, x3, r4, a, h) = (b*(3 - sqrt(3))/2, -989*b*sqrt(151316*sqrt(3) + 262087) - 3*b/2 + sqrt(3)*b/2 + 571*sqrt(3)*b*sqrt(151316*sqrt(3) + 262087), b*(3 - sqrt(3))/6, 2*b, b*(3 - sqrt(3))/8, b*(sqrt(3) + 3)/2, b*(-2967*sqrt(151316*sqrt(3) + 262087) - 2*sqrt(3) + 6 + 1713*sqrt(453948*sqrt(3) + 786261))/2)
   r1 = b*(3 - √3)/2
   r2 = b*(3√3 - 5)/2
   r3 = b*(3 - √3)/6
   x3 = 2b
   r4 = b*(3 - √3)/8
   a = b*(√3 + 3)/2
   h = b*(√3 + 3)/2
   @printf("上底が %g のとき，小円の直径は %g である（下底は %g）。\n", 2b, 2r2, 2a)
   @printf("b = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  a = %g;  h = %g\n", b, r1, r2, r3, x3, r4, a, h)
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:green, lw=0.5)
   circle(0, r1, r1)
   circle2(2r1, r1, r1)
   circle2(r1, h - r1, r1)

   circle(0, 2r1 + r2, r2, :blue)
   circle2(r1, h - 2r1 - r2, r2, :blue)

   #=
   circle(0, h - r4, r4, :green)
   circle2(r1, r4, r4, :green)

   circle2(x3, r3, r3, :orange)
   l = h - (2r1 + r2) - r4
   circle2(r1 + l*cosd(30), h - (2r1 + r2) + l*sind(30), r4, :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(2r1, r1, "大円:r1,(2r1,r1)", :red, :center, delta=-delta)
       point(0, r1, "大円:r1,(0,r1)", :red, :center, delta=-delta)
       point(r1, (1 + √3)r1, "大円:r1,(r1,(1+√3)r1)", :red, :center, delta=-delta)
       point(0, 2r1 + r2, "小円:r2,(0,2r1+r2)", :blue, :left, :vcenter, deltax=4delta)
       #=
       point(0, h - r4, "乙円:(0,h-r4)", :green, :left, :vcenter, deltax=7delta)
       point(x3, r3, "甲円:(x3,r3)", :orange, :right, :vcenter, deltax=-9delta)
       =#
       point(b, h, "(b,h)", :brown, :right, :bottom, delta=delta/2)
       point(a, 0, "a", :brown, :left, :bottom, delta=delta/2)
       plot!([r1, 2r1, r1, r1], [r1, r1, h - 2r1 - r2, r1], color=:gray70, lw=0.5)
   end  
end;

draw(14/2, true)