算額（その2121）

算額（その2121）

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

直角三角形の中に数個の等円を容れる。等円の直径と赤積が与えられたとき，等円の個数を求めよ。

注：図がカラーではないし，赤積が何なのかも明示されていないので，「赤積　= 直角三角形の面積から等円の面積（n 個分）を引いたもの」とする。しかし，そのようにしても，与える赤積の精度により解が求まらないこともあるかもしれない。そこで，まず，等円の個数が n 個のときに，鈎，股を求める一般解を求める方法から始める。

直角三角形の直角を挟む二辺の短い方を「鈎」，長い方を「股」
等円の個数を n
等円の半径と中心座標を r, (r, r), (r, 3r), ..., (x, r)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms n, r::positive,  x::positive,  鈎::positive, 股::positive;
eq1 = 4r^2 + (x - r)^2 - ((n - 2)*2r)^2
eq2 = 2r/(x - r) - 鈎/股
eq3 = dist2(0, 鈎, 股, 0, x, r, r)
res = solve([eq1, eq2, eq3], (x, 鈎, 股))[3]  # 3 of 4

   (r*(6*n^2 + 10*n*sqrt(n^2 - 4*n + 3) - 19*n - 25*sqrt(n^2 - 4*n + 3) + 4)/(5*n + 3*sqrt(n^2 - 4*n + 3) - 14), 2*r*(3*n^2 + 5*n*sqrt(n^2 - 4*n + 3) - 12*n - 14*sqrt(n^2 - 4*n + 3) + 9)/(n^2 + 3*n*sqrt(n^2 - 4*n + 3) - 4*n - 9*sqrt(n^2 - 4*n + 3) + 3), r*(n + 3*sqrt(n^2 - 4*n + 3) - 1))

鈎，股は以下の通りである。

res[2] |> println
res[3] |> println

   2*r*(3*n^2 + 5*n*sqrt(n^2 - 4*n + 3) - 12*n - 14*sqrt(n^2 - 4*n + 3) + 9)/(n^2 + 3*n*sqrt(n^2 - 4*n + 3) - 4*n - 9*sqrt(n^2 - 4*n + 3) + 3)
   r*(n + 3*sqrt(n^2 - 4*n + 3) - 1)

n = 4, 5, ... と変化させて，鈎，股，赤積（と後に導く n の推定値）を書き出す。

r = 1/2
for n = 4:20 # n = 4:449 にすれば全ての結果が表示される
   鈎 = 2*r*(3*n^2 + 5*n*sqrt(n^2 - 4*n + 3) - 12*n - 14*sqrt(n^2 - 4*n + 3) + 9)/(n^2 + 3*n*sqrt(n^2 - 4*n + 3) - 4*n - 9*sqrt(n^2 - 4*n + 3) + 3)
   股 = r*(n + 3*sqrt(n^2 - 4*n + 3) - 1)
   赤積 = 鈎*股/2 - n*pi*r^2
   推定されたn = round(赤積*0.8223 + 2.522, digits=0)
   @printf("n = %2d;  鈎 = %g;  股 = %g;  赤積 = %g;  推定値 n = %d\n", n, 鈎, 股, 赤積, 推定されたn)
   round(赤積*0.8223 + 2.522, digits=0) != n &&        println(n)
end

   n =  4;  鈎 = 2.36603;  股 = 4.09808;  赤積 = 1.70648;  推定値 n = 4
   n =  5;  鈎 = 2.20711;  股 = 6.24264;  赤積 = 2.9621;  推定値 n = 5
   n =  6;  鈎 = 2.1455;  股 = 8.30948;  赤積 = 4.20159;  推定値 n = 6
   n =  7;  鈎 = 2.11237;  股 = 10.3485;  赤積 = 5.43212;  推定値 n = 7
   n =  8;  鈎 = 2.09161;  股 = 12.3741;  赤積 = 6.65772;  推定値 n = 8
   n =  9;  鈎 = 2.07735;  股 = 14.3923;  赤積 = 7.88035;  推定値 n = 9
   n = 10;  鈎 = 2.06695;  股 = 16.4059;  赤積 = 9.10106;  推定値 n = 10
   n = 11;  鈎 = 2.05902;  股 = 18.4164;  赤積 = 10.3205;  推定値 n = 11
   n = 12;  鈎 = 2.05277;  股 = 20.4248;  赤積 = 11.539;  推定値 n = 12
   n = 13;  鈎 = 2.04772;  股 = 22.4317;  赤積 = 12.7567;  推定値 n = 13
   n = 14;  鈎 = 2.04356;  股 = 24.4374;  赤積 = 13.974;  推定値 n = 14
   n = 15;  鈎 = 2.04006;  股 = 26.4422;  赤積 = 15.1909;  推定値 n = 15
   n = 16;  鈎 = 2.03709;  股 = 28.4464;  赤積 = 16.4075;  推定値 n = 16
   n = 17;  鈎 = 2.03452;  股 = 30.4499;  赤積 = 17.6238;  推定値 n = 17
   n = 18;  鈎 = 2.03229;  股 = 32.4531;  赤積 = 18.8399;  推定値 n = 18
   n = 19;  鈎 = 2.03033;  股 = 34.4558;  赤積 = 20.0558;  推定値 n = 19
   n = 20;  鈎 = 2.02859;  股 = 36.4583;  赤積 = 21.2716;  推定値 n = 20

以上より，「等円の直径が 1 で，赤積 < 542.3553863336211 であれば，赤積を 0.8223 倍して 2.522 を加え，小数点以下を四捨五入すると等円の個数が得られる」
なお，赤積の二乗項まで使えば，予測可能領域は広がる。

n = 449
鈎 = 2*r*(3*n^2 + 5*n*sqrt(n^2 - 4*n + 3) - 12*n - 14*sqrt(n^2 - 4*n + 3) + 9)/(n^2 + 3*n*sqrt(n^2 - 4*n + 3) - 4*n - 9*sqrt(n^2 - 4*n + 3) + 3)
股 = r*(n + 3*sqrt(n^2 - 4*n + 3) - 1)
赤積 = 鈎*股/2 - n*pi*r^2
#@printf("n = %2d;  鈎 = %g;  股 = %g;  赤積 = %g\n", n, 鈎, 股, 赤積)
println(赤積, ";  ", round(赤積*0.8223 + 2.522, digits=0))

   542.3553863336211;  449.0

以下の描画プログラムは，精度の良くない赤積の値を与えると，それに最も近い赤積を与える n を求める。

r = 1/2, 赤積 = 9 を与えると「与えられた赤積 = 9;  等円の直径 = 2;  正確な赤積 = 6.82593;  n = 4」を表示する。

もし，等円の半径が 1/2 以外のときは，面積は半径の二乗に比例するので，換算して計算する。
r = 4, n = 10 のときは赤積は r = 1/2 のときの 8^2 倍なので，8^2 * 9.10106 = 582 を与えると，「与えられた赤積 = 582;  等円の直径 = 8;  正確な赤積 = 582.468;  n = 10」を表示する。

function draw(r0, 赤積0, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   f = 2r0
   赤積 = 赤積0 / f^2
   n = round(赤積*0.8223 + 2.522, digits=0)
   (x, 鈎, 股) = (r*(6*n^2 + 10*n*sqrt(n^2 - 4*n + 3) - 19*n - 25*sqrt(n^2 - 4*n + 3) + 4)/(5*n + 3*sqrt(n^2 - 4*n + 3) - 14), 2*r*(3*n^2 + 5*n*sqrt(n^2 - 4*n + 3) - 12*n - 14*sqrt(n^2 - 4*n + 3) + 9)/(n^2 + 3*n*sqrt(n^2 - 4*n + 3) - 4*n - 9*sqrt(n^2 - 4*n + 3) + 3), r*(n + 3*sqrt(n^2 - 4*n + 3) - 1))
   赤積 = 鈎*股/2 - n*pi*r^2
   string = @sprintf("与えられた赤積 = %g; 等円の直径 = %g\n正確な赤積 = %g;  n = %d", 赤積0, 2r0, 赤積*f^2, n)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
   circle(r, r, r)
   for i = 0:n - 2
       circle(r + (x - r)/(n - 2)*i, 3r - 2r/(n - 2)*i, r)
   end
   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, 鈎, string, :black, :left, :bottom, delta=delta, mark=false)
   end
end;

draw(4, 582, true)