算額（その428）

算額（その428）

長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

中村信弥(1999)：算額への招待
http://www.wasan.jp/syotai/syotai.html

外円の中に長方形が内接している。一辺の長さが 41.86 寸の正方形が 6個，大円が 2 個，小円が 4 個入っている。小円の直径を求めよ。

正方形の一辺の長さを a
大円の半径と中心座標を r1, (r0 - r1)
小円の半径と中心座標を r2, (r0 - 2r1 + r2, y)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r0::positive, r1::positive,
     r2::positive, y::positive;
@syms a, r0, r1, r2, y
eq1 = (sqrt(Sym(2))a/2 + a)^2 + (r0 - sqrt(Sym(2))a + a)^2 - r0^2
eq2 = (r0 - 2r1)^2 + (r0 - sqrt(Sym(2))a)^2 - r0^2
eq3 = (r0 - 2r1 + r2)^2 + y^2 - (r0 - r2)^2
eq4 = (r1 - r2)^2 + y^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3, eq4], (r0, r1, r2, y));

2 組の解が得られるが，2 番目のものが適解である。

names = ["r0", "r1", "r2", "y"]
for (i, name) in enumerate(names)
   println("\n", name)
   println(simplify(res[2][i]))
   println(res[1][i](a => 4186).evalf())
end

   r0
   a*(5 + 7*sqrt(2))/4
   15592.3214511641

   r1
   (-23922439916777731675282139353831061069516954103552970349841905*a^2 + 16915719487681281681493317896073504500369711065425523483418778*sqrt(2)*a^2 - 6986054008647705185868918429789392*sqrt(5)*sqrt(1144573653449326586061314076037948929840593626875394274 - 809335791921480250432937095033206217367486276799679753*sqrt(2))*sqrt(a^4) + 4939886163250256057660601952272630*sqrt(10)*sqrt(1144573653449326586061314076037948929840593626875394274 - 809335791921480250432937095033206217367486276799679753*sqrt(2))*sqrt(a^4))/(8*a*(4882633868649679478319412976906635182914530485393454527631129 - 3452543518573294933348514588454588356006400466465731641999825*sqrt(2)))
   1681.32809749726

   r2
   a*(205 + 151*sqrt(2))/1168
   1500.02961796760

   y
   -sqrt(-6848627247131292427665721780*sqrt(2)*a^2 + 9685421536530988262785792422*a^2 - 4*sqrt(10)*sqrt(1144573653449326586061314076037948929840593626875394274 - 809335791921480250432937095033206217367486276799679753*sqrt(2))*sqrt(a^4))/(8*sqrt(8442762866313367764895228377 - 5969934874720155329318365400*sqrt(2)))
   3176.18761647798

小円の半径は a*(205 + 151*√2)/1168 なので，直径は a*(205 + 151*√2)/584 である。

a = 41.86 を代入すると，小円の直径は 30 寸あまりあり。

41.86*(205 + 151*√2)/584

   30.00059235935206

res2 = Float64.([res[2][i](a => 4186).evalf()/100 for i in 1:4])

   4-element Vector{Float64}:
    155.92321451164108
     16.81328097497264
     15.00029617967603
     31.761876164779796

r0 = 155.923;  r1 = 16.8133;  r2 = 15.0003;  y = 31.7619
小円の直径 = 30.0006

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, r2, y) = res2
   @printf("r0 = %g;  r1 = %g;  r2 = %g;  y = %g\n", r0, r1, r2, y)
   @printf("小円の直径 = %g\n", 2r2)
   a = 41.86
   (x0 , y0) = (r0 - 2r1, r0 - √2a)
   plot([x0, x0, -x0, -x0, x0], [-y0, y0, y0, -y0, -y0], color=:orange, lw=0.5)
   plot!([0, √2a/2, 0, -√2a/2, 0], [r0 - √2a, r0 - √2a + √2a/2, r0, r0 - √2a + √2a/2, r0 - √2a], color=:red, lw=0.5)
   rect(√2a/2, r0 - √2a, √2a/2 + a, r0 - √2a + a, :red)
   rect(-√2a/2, r0 - √2a, -√2a/2 - a, r0 - √2a + a, :red)
   plot!([0, √2a/2, 0, -√2a/2, 0], -[r0 - √2a, r0 - √2a + √2a/2, r0, r0 - √2a + √2a/2, r0 - √2a], color=:red, lw=0.5)
   rect(√2a/2, -r0 + √2a, √2a/2 + a, -r0 + √2a - a, :red)
   rect(-√2a/2, -r0 + √2a, -√2a/2 - a, -r0 + √2a - a, :red)
   circle(0, 0, r0, :black)
   circle(r0 - r1, 0, r1, :blue)
   circle(-r0 + r1, 0, r1, :blue)
   circle4(r0 - 2r1 + r2, y, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r0, 0, " r0", :black, :left, :bottom, delta=delta)
       point(r0 - r1, 0, "r0-r1", :blue, :center, delta=-delta/2)
       point(r0 - 2r1 + r2, y, "(r0-2r1+r2,y) ", :green, :right, :vcenter)
       point(0, r0 - √2a, "r0-√2a ", :orange, :right, :top, delta=-delta/2)
       point(a/√2, r0-√2a, "(a/√2,r0-√2a)", :red, :center, delta=-3delta)
       point(a/√2 + a, r0-√2a, "(a/√2+a,r0-√2a)", :red, :center, delta=-delta/2)
       point(a/√2 + a, r0 - √2a + a, "(a/√2+a,r0-√2a+a)", :red, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その427）

算額（その427）

長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

中村信弥(1999)：算額への招待
http://www.wasan.jp/syotai/syotai.html

大円は鈎股弦（直角三角形）の，鈎と弦に接しており，中円，小円は弦と大円に接しかつそれぞれも互いに接している。鈎，股がそれぞれ　3000 寸，4000 寸，小円の直径が 456 寸のとき，大円の直径を求めよ。

鈎，股，弦の長さをそれぞれ鈎，股，弦
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (x2, x2)
小円の半径と中心座標を r3, (x31, y31), (x32, y32)
として，連立方程式を解く。

solve() の能力的に一度では解けないので，まず最初に大円と中円の半径を求める。このためには，図を反時計回りに回転させ，弦（下図の紫の線）が水平になるように回転させて方程式を解くとよい。eq13 はr1 と鈎股弦の面積の関係式である（上図を参照）。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     r1::positive, r2::positive,
     r3::positive, x3::positive;

eq11 = x3^2 + (r1 - 2r2 + r3)^2 - (r1 - r3)^2
eq12 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2
eq13 = r1*(鈎 + 股) + (r1 - 2r2)*弦 - 鈎*股
res1 = solve([eq11, eq12, eq13], (r1, r2, x3))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    ((-2*r3*弦^2 - 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 - sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), -sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) + 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))
    ((-2*r3*弦^2 - 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 - sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) + 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))
    ((-2*r3*弦^2 + 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 + sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), -sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) - 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))
    ((-2*r3*弦^2 + 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 + sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) - 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))

4 組の解が得られるが，そのうちの最初の 2 組が適解である。1 番目と 2 番目は，小円の位置（x3）が対称なもので実質的に同じである。

println("r1 = ", res1[2][1](鈎 => 3000, 股 => 4000, 弦 => 5000, r3 => 456//2))
println("r2 = ", res1[2][2](鈎 => 3000, 股 => 4000, 弦 => 5000, r3 => 456//2))
println("x3 = ", res1[2][3](鈎 => 3000, 股 => 4000, 弦 => 5000, r3 => 456//2))

   r1 = 1250
   r2 = 300
   x3 = 120*sqrt(19)

算額の解としては大円の半径が 1250，直径は 2500 でここまででよい。

次いで，図を描くために，r1, r2 が既知として，小円の中心座標を求める。

using SymPy
@syms r1::positive, r2::positive, x2::positive, y2::positive,
     r3::positive, x31::positive, y31::positive,
     x32::positive, y32::positive;

(r1, r2, r3) = (1250, 300, 456//2)
x2 = r1 + (r1 - r2)*3//5
y2 = r1 + (r1 - r2)*4//5
eq1 = (x2 - x31)^2 + (y2 - y31)^2 - (r2 + r3)^2
eq2 = (x2 - x32)^2 + (y2 - y32)^2 - (r2 + r3)^2
eq5 = (x31 - r1)^2 + (y31 - r1)^2 - (r1 - r3)^2
eq6 = (x32 - r1)^2 + (y32 - r1)^2 - (r1 - r3)^2
res2 = solve([eq1, eq2, eq5, eq6], (x31, y31, x32, y32))

   4-element Vector{NTuple{4, Sym}}:
    (8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5, 8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5)
    (8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5, 96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19))
    (96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19), 8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5)
    (96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19), 96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19))

4 組の解が得られるが，2 番目と 3 番目のものが適解である（小円の位置関係が逆になるだけ）。

(8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5, 96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19)) |> println
(96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19), 8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5) |> println

   (1358.3457014200953, 2266.2407239349286, 2195.2542985799046, 1638.5592760650716)
   (2195.2542985799046, 1638.5592760650716, 1358.3457014200953, 2266.2407239349286)

2 つの解をまとめると，以下のようになる。

   鈎 = 3000;  股 = 4000;  r3 = 228
   r1 = 1250;  r2 = 300;  x31 = 2195.25;  y31 = 1638.56;  x32 = 1358.35;  y32 = 2266.24
   大円の直径 = 2500

using Plots

function draw0(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, x3) = (1250, 300, 456//2, 120*sqrt(19))
   plot()
   circle(0, 0, r1)
   circle(0, r1 - r2, r2, :blue)
   circle(x3, r1 - 2r2 + r3, r3, :green)
   circle(-x3, r1 - 2r2 + r3, r3, :green)
   x = sqrt(r1^2 - (r1 - 2r2)^2)
   segment(-x, r1 - 2r2, x, r1 - 2r2, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r1-r2, " r1-r2", :blue)
       point(0, r1-2r2, " r1-2r2", :magenta, delta=-delta)
       point(x3, r1 - 2r2 + r3, "(x3,r1-2r2+r3)", :green, :center, :top, delta=-delta)
       point(r1, 0, "r1 ", :red, :right, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, r3) = (3000, 4000, 456//2)
   (r1, r2) = (1250, 300)
   (x31, y31, x32, y32) = res2[3]
   x2 = r1 + (r1 - r2)*3/5
   y2 = r1 + (r1 - r2)*4/5
   @printf("鈎 = %g;  股 = %g;  r3 = %g\n", 鈎, 股, r3)
   @printf("r1 = %g;  r2 = %g;  x31 = %g;  y31 = %g;  x32 = %g;  y32 = %g\n", r1, r2, x31, y31, x32, y32)
   @printf("大円の直径 = %.8g\n", 2r1)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:magenta, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, y2, r2, :blue)
   circle(x31, y31, r3, :green)
   circle(x32, y32, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1, r1, " 大円:r1,(r1,r1)", :red, :left, :vcenter)
       point(r1 + (r1 - r2)*3/5, r1 + (r1 - r2)*4/5, " 中円:r2,(x2,y2)", :blue, :left, :vcenter)
       point(x31, y31, " 小円:r3,(x31,y31)", :green, :left, :vcenter)
       point(x32, y32, " 小円:r3,(x32,y32)", :green, :left, :bottom, delta=delta)
       point(股, 0, "股", :black, :left, :bottom, delta=delta)
       point(0, 鈎, " 鈎", :black, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;