算額（その458）

算額（その458）

長野県長野市 久保寺観音堂 享和3年(1803)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

円に正三角形が内接している。
正三角形の一辺と外積（円と正三角形の間の面積）の和を A，矢（弦の中点と弧の中点を結ぶ線分の長さ）と円周の和を B とするとき，円の面積を求めよ。

最初に，簡単な解を求める方法について書く。条件 A は不要で，かつ，条件 B はとてつもなく簡単なものである。

図において，「矢」は直径の 1/4 である。
したがって，円の半径を r とすれば，「矢と円周の和を B とする」は 「r/2 + 2r*π = B」 である。
半径は，この方程式を r について解けば得られる。

include("julia-source.txt")

using SymPy
@syms r, B
solve(r/2 + 2r*PI - B, r)[1] |> println

   2*B/(1 + 4*pi)

円の半径が 2B/(1 + 4π) なので，面積は π*(2B/(1 + 4π))^2 である。

B = 8.139822368615503
r = 2B/(1 + 4π)
(r, π*(2B/(1 + 4π))^2)

   (1.2, 4.523893421169302)

---

まともにやるととてつもなく面倒なことになる。

円の半径を r，正三角形の一辺の長さを 2a とする。
円，正三角形の面積をそれぞれ s1, s2 とおく。

include("julia-source.txt")

using SymPy

@syms r::positive, a::positive, s1::positive, s2::positive, 矢::positive, A::positive, B::positive
s1 = PI*r^2
a = r*sqrt(Sym(3))/2
s2 = a^2*sqrt(Sym(3))
矢 = r - r/2;

A, B を定義する。

eq1 = 2a + (s1 - s2) - A
eq2 = 矢 + PI*2r - B;

差を取って，r について解くのが常套手段ではあるが，eq2 を見るとこれを r について解けばよいことがわかる。

solve(eq2, r) |> println

   Sym[2*B/(1 + 4*pi)]

ということで，冒頭の簡単な方法へ戻る。

なお，eq1 は不要であるし，差を取ったのでは条件 A が残り，面倒な式になる。無理にやると解は二通り得られるが，どちらが適解かは単純にはいえない。以下に示すように，どちらを取るかは r  に依存するのである。これでは自分の尻尾を飲み込もうとしている蛇と同じである。

eq3 = eq1 - eq2
res = solve(eq3, r);

res[1] |> println

   (-sqrt(-12*sqrt(3)*A + 16*pi*A - 16*pi*B + 12*sqrt(3)*B - 16*sqrt(3)*pi - 4*sqrt(3) + 13 + 8*pi + 16*pi^2) - 4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3))

res[2] |> println

   (sqrt(-12*sqrt(3)*A + 16*pi*A - 16*pi*B + 12*sqrt(3)*B - 16*sqrt(3)*pi - 4*sqrt(3) + 13 + 8*pi + 16*pi^2) - 4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3))

円の半径を与えて，A, B を計算し，それを eq3 の解 res により円の半径を再計算する関数 rev を書いてみる。
r > (-4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3) の場合には最初のもの，そうでない場合には2番めのものが適解である。

function rev(r)  # 半径を与えて検算する
   s1 = pi*r^2
   a = r*sqrt(3)/2
   s2 = a^2*sqrt(3)
   矢 = r - r/2
   A0 = 2a + (s1 - s2)
   B0 = 矢 + pi*2r
   println(A0, ", ", B0)
   term1 = sqrt(-12*sqrt(3)*A0 + 16*pi*A0 - 16*pi*B0 + 12*sqrt(3)*B0 - 16*sqrt(3)*pi - 4*sqrt(3) + 13 + 8*pi + 16*pi^2)
   term2 = - 4*pi - 1 + 2*sqrt(3)
   ans1 = (-term1 + term2)/(-4*pi + 3*sqrt(3))
   ans2 = ( term1 + term2)/(-4*pi + 3*sqrt(3))
   return (r > (-4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3)) ? ans1 : ans2)
end;

半径が 1.2 の場合

rev(1.2)

   4.731739518077568, 8.139822368615503

   1.1999999999999988

半径が 3.45 の場合

rev(3.45)

   27.906580792648718, 23.401989309769576

   3.4499999999999993

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1.2
   a = r*√3/2
   b = -r/2
   c = 3r*(2 - √3)/2
   d = r*(11 - 6*√3)/2
   plot([a, 0, -a, a], [-r/2, r, -r/2, -r/2], color=:blue, lw=0.5)
   circle(0, 0, r)
   segment(0, -r/2, 0, -r, :green, lw=2)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r, 0, "r ", :red, :right, :bottom, delta=delta/2)
       point(a, -r/2, "(a,-r/2) ", :blue, :right, :bottom, delta=delta/2)
       point(0, -r/2, " -r/2", :blue, :left, :bottom, delta=delta/2)
       point(0, -3r/4, " 矢", :green, mark=false)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その457）

算額（その457）

長野県長野市 久保寺観音堂 享和3年(1803)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

円に正三角形，正三角形に正方形が内接している。円の直径，正三角形の辺，正方形の辺のそれぞれの 3 乗の和を A，円の直径，正方形の辺の 3 乗の差を B とするとき，円の直径を A, B で表わせ。

include("julia-source.txt")

using SymPy

円の半径，正三角形の辺，正方形の辺をそれぞれ r, 2a, 2c とおく。
正方形の下辺と上辺がy軸と交差するy座標を b, d とする。
a = r√3/2, b = -r/2 から，c, d を求める。

@syms a, b, c, d, r
a = r*√Sym(3)/2
b = -r/2
eq1 = d - b - 2c
eq2 = 2c/(a - c) - sqrt(Sym(3))
solve([eq1, eq2], (c, d))

   Dict{Any, Any} with 2 entries:
     d => -3*sqrt(3)*r + 11*r/2
     c => -3*sqrt(3)*r/2 + 3*r

r, a, c と A, B の関係式を記述する。

@syms r1::positive, r2::positive,
     r3::positive, x3::positive,
     A::positive, B::positive;
a = r*√Sym(3)/2
c = 3r*(2 - √Sym(3))/2
eq1 = (2r)^3 + (2a)^3 + (2c)^3 - A
eq2 = (2r)^3 - (2c)^3 - B;

eq1 + eq2 としてA, B, r の関係式を得る。

eq3 = eq1 + eq2
eq3 |> println

   -A - B + 3*sqrt(3)*r^3 + 16*r^3

eq3 を r について解く（r を A, B で表す式を求める）。

solve(eq3, r)[1] |> factor |> println

   (A + B)^(1/3)/(3*sqrt(3) + 16)^(1/3)

(中身)^(1/3) の形なので，3 乗根の中身を簡約化する。

@syms x, A::positive, B::positive;
apart((A + B)/(3sqrt(Sym(3)) + 16), x) |> simplify |>factor |> println

   -(-16 + 3*sqrt(3))*(A + B)/229

求める関係式は (16 - 3√3)*(A + B)/229 の 3 乗根である。
求まるのは半径なので，2倍して直径にする関数 f を定義する。

f(A, B) = 2(A + B)^(1/3)/(3*sqrt(3) + 16)^(1/3);

検算してみる。

r = 1; a = r*√3/2; c = 3r*(2 - √3)/2 のとき，
A = 13.715575357311328; B = 7.480577065395304
である。

r = 1
a = r*√3/2
c = 3r*(2 - √3)/2
A = (2r)^3 + (2a)^3 + (2c)^3
B = (2r)^3 - (2c)^3
(A, B) |> println

   (13.715575357311328, 7.480577065395304)

計算された A, B から f() によって円の直径を求める。
確かに 2.0 になる。

f(13.715575357311328, 7.480577065395304)

   2.0

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1
   a = r*√3/2
   b = -r/2
   c = 3r*(2 - √3)/2
   d = r*(11 - 6*√3)/2
   plot([a, 0, -a, a], [-r/2, r, -r/2, -r/2], color=:blue, lw=0.5)
   plot!([c, c, -c, -c, c], [-r/2, d, d, -r/2, -r/2], color=:green, lw=0.5)
   circle(0, 0, r)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r, 0, "r ", :red, :right, :bottom, delta=delta/2)
       point(0, b, " b", :blue, :left, :bottom, delta=delta/2)
       point(0, d, " d", :green, :left, :bottom, delta=delta/2)
       point(c, 0, "c ", :green, :right, :bottom, delta=delta/2)
       point(a, 0, "a ", :blue, :right, :bottom, delta=delta/2)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;