算額（その685）

算額（その685）

和算問題あれこれ 2 令和5年2月の問題-No.4（『算法求積通考』第6問）
https://gunmawasan.web.fc2.com/k-n-mondai.html

長径 2 寸，短径 1 寸の楕円形の周は何寸か。
小数点以下 4 位を四捨五入し，小数点以下 3 位まで答えよ。
円周率は 3.14159 を使用のこと。

コンピュータを使って正確な値を求めるならば，第二種完全楕円積分を使えばよい。
計算するためのパッケージ（ライブラリ）が用意されているのでそれを使えばよい。あるいは，区分求積法を使ってもよいし，逐次計算により求めてもよい。

しかし，今回のように要求される精度がほどほどに指定されている場合には，いくつかの近似式があるので精度に基づいて選択すればよい。

以下では，楕円のパラメータとして，長半径 a，短半径 b を使う。
長半径，短半径は長径，短径の半分である。

a = 2//2   # 長径 2
b = 1//2;  # 短径 1

まずは，関孝和の近似式

2sqrt(4(a - b)^2 + 3.14159^2*a*b)

    4.872286471072899

第二種完全楕円積分の近似式で，Gauss-Kummer の公式

h = (a - b)/(a + b)
3.14159*(a + b)*(1 + (h/2)^2 + h^4/64 + h^6/256)

    4.844218858908822

同じく近似式でシュリニヴァーサ・アイヤンガー・ラマヌジャンの近似式

3.14159*(3(a + b) - sqrt((3a + b)*(a + 3b)))

    4.844206457106038

名前が似ているが，もう少し精度が高いらしい，シュリニヴァーサ・ラマヌジャンの近似式

3.14159* (a + b)*(1 + 3*((a - b)/(a + b))^2/(10 + sqrt(4 - 3*((a - b)/(a + b))^2)))

    4.844220016323984

残念ながら，関孝和の近似式はちょっと精度が足りないが，他の近似公式では「小数点以下 4 位を四捨五入し，小数点以下 3 位まで答えよ」の要件を満たし，4.844 を与える。

なお，第二種完全楕円積分を用いた正確な値は 4.844224110273837 である。

using Elliptic
m = 1.0 - (b / a)^2             # 楕円の離心率を計算
println(4 * a * ellipke(m)[2])  # 楕円の周長を計算(第二種完全楕円積分）

    4.844224110273837

算額（その684）

算額（その684）

和算問題あれこれ 2 令和5年2月の問題-No.2（茨城県 板橋不動願成寺）
https://gunmawasan.web.fc2.com/k-n-mondai.html

正三角形の中に，大円，中円，小円が入っており，中円と小円の共通接線となる斜線がある。
小円の直径が 9 寸のとき，中円の直径を求めよ。

正三角形の左の頂点を原点とする
正三角形の一辺の長さを a
大円の半径と中心座標を r1, (a/2, r1); r1 = √3*a/6
中円の半径と中心座標を r2, (a/2, 2r1 + r2)
小円の半径と中心座標を r3, (x3, y3)
斜線と正三角形の交点座標を (x1, y1) 
とおき，以下の連立方程式を解く。

(1)　まず，正三角形の一辺の長さ a と，大円，中円の半径 r1, r2 の関連を知る。

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

using SymPy
@syms a::positive, r1::positive, r2::positive

sqrt3 = sqrt(Sym(3))
eq1 = (sqrt3*a/2 - r1)/2 - r1
eq2 = (sqrt3*a/2 - 2r1 - r2)/2 - r2
solve([eq1, eq2], (r1, r2))

   Dict{Any, Any} with 2 entries:
     r1 => sqrt(3)*a/6
     r2 => sqrt(3)*a/18

(2)　次に，小円の半径がわかっているので，小円の中心座標 (x3, y3) と正三角形の一辺の長さ a を求める

@syms a::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, y3::positive,
     x1::positive
y1 = sqrt3*(a - x1)
r1 = sqrt3*a/6
r2 = r1/3
r3 = 9//2

eq3 = (a/2 - x3)^2 + (sqrt3*a/6 - y3)^2 - (r1 + r3)^2
eq5 = y3*a/2 - (2r1 + r2)*x3
eq6 = (a/2 - 2sqrt(r1*r3))*r2 - (a/2 +  2sqrt(r1*r2))*(9//2);

res = solve([eq3, eq5, eq6], (a, x3, y3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (75*sqrt(3), 27*sqrt(3)/2, 63/2)
    (75*sqrt(3), 1287*sqrt(3)/38, 3003/38)

2 組の解が得られるが，最初のものが適解である。

(3)　最後に斜線と正三角形の交点座標 (x1, y1) を求める。

@syms a::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, y3::positive,
     x1::positive
(a, x3, y3) = (75*sqrt3, 27*sqrt3/2, 63//2)
y1 = sqrt3*(a - x1)
(r1, r2, r3) = (sqrt3*a/6, sqrt3*a/18, 9//2)
eq4 = dist(0, 0, x1, y1, a/2, 2r1 + r2) - r2^2
res2 = solve(eq4)
res2 |> println

   Sym[75*sqrt(3)/2, 325*sqrt(3)/7]

2番目のものが適解である。x1 = 325√3/7 = 80.4166446371264

x1 = res2[2].evalf()
x1 |> println

   80.4166446371264

y 座標値は y1 = √3(a - x1) = 85.7142857142857

y1 = √3(a - x1).evalf()
y1 |> println

   85.7142857142857

以上のパラメータをまとめると以下のようになる。

   中円の直径 = 25;  a = 129.904;  r1 = 37.5;  r2 = 12.5;  r3 = 4.5;  x3 = 23.3827;  y3 = 31.5;  x1 = 80.4166;  y1 = 85.7143

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   sqrt3 = √3
   (a, x3, y3) = (75√3, 27√3/2, 63/2)
   (r1, r2, r3) = (√3a/6, √3a/18, 9/2)
   x1 = 325*sqrt(3)/7
   y1 = sqrt3*(a - x1)
   @printf("中円の直径 = %.15g;  a = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  x1 = %g;  y1 = %g\n",
       2r2, a, r1, r2, r3, x3, y3, x1, y1)
   plot([0, a, a/2, 0], [0, 0, sqrt3*a/2, 0], color=:black, lw=0.5)
   circle(a/2, r1, r1)
   circle(a/2, 2r1 + r2, r2, :blue)
   circle(x3, y3, r3, :green)
   segment(0, 0, x1, y1, :orange)
   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/2, r1, "大円:r1,(a/2,r1)", :red, :center, delta=-delta)
       point(a/2, 2r1 + r2, "中円:r2,(a/2,2r1+r2) ", :blue, :right, :vcenter)
       point(x3, y3, " 小円:r3,(x3,y3)", :blue, :left, :top, delta=-2delta)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(a/2, 0, "a/2", :black, :center, :bottom, delta=delta/2)
       point(a/2, √3a/2, "(a/2,√3a/2)", :black, :left, :bottom, delta=delta/2)
       point(x1, y1, " (x1,y1)", :black, :left, :vcenter)
   end
end;

算額（その683）

算額（その683）

和算問題あれこれ 2 令和5年2月の問題-No.1（茨城県 板橋不動願成寺）
https://gunmawasan.web.fc2.com/k-n-mondai.html

直角三角形の中に正方形，長方形，等円 2 個がある。
鈎，股 がそれぞれ 27 寸，36 寸のとき，等円の直径を求めよ。

この問題は，算額（70）
https://blog.goo.ne.jp/r-de-r/e/071461791de5889b7aaceea141d47ced
東京都府中市 大国魂神社
http://www2.ttcn.ne.jp/~nagai/sangaku/sangakumondai1.htm
のものに似ているが，鈎，股に接する等円が正方形，長方形との接し方に違いがある。

式の記述を簡単にするために，問の図を左右を反転させた図で考える。

鈎，股，弦の長さをを「鈎」，「股」，「弦」
長方形の長辺，短辺の長さを「長」，「短」
正方形の一辺の長さを「方」
等円のと中心座標を「径」，(x1, 短 + 径/2)
として，以下の連立方程式を解く。

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, 径::positive,
     x1

(鈎, 股) = (Sym(27), Sym(36))
弦 = sqrt(鈎^2 + 股^2)
eq1 = (1 + 4//3 - 5//3)方 - 径
eq2 = 鈎 - (5//4 + 3//5)*方 - 短
eq3 = (弦 - (3//4 + 1 + 4//3)*方)^2 - 短^2 - (股 - 径 - 長)^2
eq4 = dist(4方/5, 短, 0, 鈎 - 5方/4, 径/2, 径/2) - 径^2/4
eq5 = dist(股, 0, 0, 鈎, x1, 短 + 径/2) - 径^2/4
res = solve([eq1, eq2, eq3, eq4, eq5], (方, 長, 短, 径, x1))

   4-element Vector{NTuple{5, Sym}}:
    (12, 108/5, 24/5, 8, 88/5)
    (12, 108/5, 24/5, 8, 464/15)
    (12, 172/5, 24/5, 8, 88/5)
    (12, 172/5, 24/5, 8, 464/15)

4 組の解が得られるが，最初のものが適解である。

鈎，股 が 27寸, 36寸のとき，等円の直径は 8 寸である。

正方形の一辺の長さは 12寸，長方形の長辺，短辺の長さはそれぞれ 108/5 = 21.6寸，24/5 = 4.8寸である。

大国魂神社の算額と比較すると，鈎股に接する等円と長方形の短辺が接するように長方形の長辺を長くしただけであることがわかる。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股) = (27, 36)
   弦 = sqrt(鈎^2 + 股^2)
   (方, 長, 短, 径, x1) = (12, 108/5, 24/5, 8, 88/5)
   r = 径/2
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(r, r, r)
   circle(x1, 短 + r, r)
   rect(2r + 長, 0, 2r, 短, :blue)
   plot!([4方/5, 0, 3方/5, 7方/5, 4方/5],
       [短, 鈎 -5方/4, 鈎 - 9方/20, 鈎 - 21方/20, 短], color=:green, lw=0.5)
   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, " 股", :black, :left, :bottom, delta=delta)
       point(0, 鈎, " 鈎", :black, :left, :bottom, delta=delta/2)
       point(r, r, "等円:r,(r,r)", :red, :center, delta=-delta)
       point(x1, 短 + r, "等円:r\n(x1,短+r)", :red, :center, delta=-delta)
       point(2r + 長, 0, " 2r+長", :blue, :left, :bottom, delta=delta)
       point(4方/5, 短, "a", :green, :center, delta=-delta)
       point(0, 鈎 -5方/4, "  b", :green, :left, :vcenter)
       point(3方/5, 鈎 - 9方/20, " c", :green, :left, :bottom)
       point(7方/5, 鈎 - 21方/20, " d", :green, :left, :bottom)
   end
end;