算額（その515）

算額（その515）

横川良助直胤: 神壁算法追加，文化4年(1807)
一関市博物館>>和算に挑戦>>平成27年度出題問題>>平成27年度出題問題(2)[中級問題]&解答例
https://www.city.ichinoseki.iwate.jp/museum/wasan/h27/normal.html

上底と下底が 4寸5分，8寸の等脚台形の中に円が内接しており，その円に正方形が内接しており，さらにその正方形に円が内接しており...と繰り返されている。
。4 番目の円に内接する正方形の一辺の長さを求めよ。なお，最後の正方形の辺は座標軸に水平・垂直である。

上底，仮定を 2a, 2b とする。
台形の高さ h は sqrt((a + b)^2 - (b - a)^2) である。
よって，第1円の半径は h/2 である。
第1正方形の対角線の長さが第1円の半径に等しい。
第n円の直径は第n-1円の直径の 1/√2 である。
以下，この繰り返し。
最後の正方形は45度回転して描く。

   a = 2.25;  b = 4;  h = 6;  r1 = 3
   1 番目の内接円の半径 = 3
   2 番目の内接円の半径 = 2.12132
   3 番目の内接円の半径 = 1.5
   4 番目の内接円の半径 = 1.06066
   最後の正方形の一辺の長さは 1.5（1 寸 5 分）

include("julia-source.txt");

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (45, 80) .// 20
   h = sqrt((a + b)^2 - (b - a)^2)
   r1 = h/2
   r = zeros(10)
   @printf("a = %g;  b = %g;  h = %g;  r1 = %g\n", a, b, h, r1)
   plot([b, a, -a, -b, b], [0, h, h, 0, 0], color=:black, lw=0.5)
   r[1] = h/2
   n = 4
   for i = 1:n - 1
       @printf("%d 番目の内接円の半径 = %g\n", i, r[i])
       circle(0, h/2, r[i])
       plot!([r[i], 0, -r[i], 0, r[i]], [h/2, h/2 + r[i], h/2, h/2 - r[i], h/2], color=:blue, lw=0.5)
       r[i + 1] = r[i]/√2
   end
   @printf("%d 番目の内接円の半径 = %g\n", n, r[n])
   circle(0, h/2, r[n])
   t = r[n]*√2/2
   plot!([t, -t, -t, t], [h/2 + t, h/2 + t, h/2 - t, h/2 - t, h/2 + t], color=:green, lw=0.5)
   @printf("最後の正方形の一辺の長さは %g\n", 2t)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(a, h, " (a,h)", :black, :left, :vcenter)
       point(b, 0, "b ", :black, :right, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;

算額（その514）

算額（その514）

横川良助直胤: 神壁算法追加，文化4年(1807)
一関市博物館>>和算に挑戦>>平成27年度出題問題>>平成27年度出題問題(2)[中級問題]&解答例
https://www.city.ichinoseki.iwate.jp/museum/wasan/h27/normal.html

大円と小円が交わっており，上部の小円部分に甲円 1 個，乙円 2 個が入っている。

下図は説明のためのものなので，問の条件のときのものとは見た目が違う。

小円，甲円，乙円の直径がそれぞれ，8 寸，5 寸，3 寸のとき，大円の直径はいかほどか。

大円の直径と中心座標を r0, (0, 0)
小円の直径と中心座標を r1, (0, r0 + 2r2 - r1)
甲円の直径と中心座標を r2, (0, r0 + r2)
乙円の直径と中心座標を r3, (x3, y3)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, r2::positive, r3::positive, x3::positive, y3::positive

(r1, r2, r3) = (8, 3, 2) .// 2
eq1 = x3^2 + y3^2 - (r0 + r3)^2
eq2 = x3^2 + (y3 - (r0+ 2r2 - r1))^2 - (r1 - r3)^2
eq3 = x3^2 + (r0 + r2 - y3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (r0, x3, y3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-r2^2*(r1 - r2 - r3)/(r1*r2 - r1*r3 - r2^2), -2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 - r2 - r3)/(r1 - r2), (r1^2*r2^2 - 3*r1^2*r2*r3 + r1^2*r3^2 - 2*r1*r2^3 + 3*r1*r2^2*r3 + r1*r2*r3^2 + r2^4)/((r1 - r2)*(r1*r2 - r1*r3 - r2^2)))
    (-r2^2*(r1 - r2 - r3)/(r1*r2 - r1*r3 - r2^2), 2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 - r2 - r3)/(r1 - r2), (r1^2*r2^2 - 3*r1^2*r2*r3 + r1^2*r3^2 - 2*r1*r2^3 + 3*r1*r2^2*r3 + r1*r2*r3^2 + r2^4)/((r1 - r2)*(r1*r2 - r1*r3 - r2^2)))

2 組の解が得られるが，2 番目のものが適解である。
小円，甲円，乙円の直径がそれぞれ，8 寸，5 寸，3 寸のとき，大円の直径は 27 寸である。

(r1, r2, r3) = (8, 3, 2) .// 2
(
   -r2^2*(r1 - r2 - r3)/(r1*r2 - r1*r3 - r2^2),
   2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 - r2 - r3)/(r1 - r2),
   (r1^2*r2^2 - 3*r1^2*r2*r3 + r1^2*r3^2 - 2*r1*r2^3 + 3*r1*r2^2*r3 + r1*r2*r3^2 + r2^4)/((r1 - r2)*(r1*r2 - r1*r3 - r2^2))
)

   (27//2, 2.3999999999999995, 143//10)

using Plots
function draw(r1, r2, r3, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, x3, y3) = (-r2^2*(r1 - r2 - r3)/(r1*r2 - r1*r3 - r2^2), 2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 - r2 - r3)/(r1 - r2), (r1^2*r2^2 - 3*r1^2*r2*r3 + r1^2*r3^2 - 2*r1*r2^3 + 3*r1*r2^2*r3 + r1*r2*r3^2 + r2^4)/((r1 - r2)*(r1*r2 - r1*r3 - r2^2)))
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  r0 = %g;  x3 = %g;  y3 = %g\n", r1, r2, r3, r0, x3, y3)
   @printf("大円の直径 = %g\n", 2r0)
   plot()
   circle(0, 0, r0)
   circle(0, r0 + 2r2 - r1, r1, :blue)
   circle(0, r0 + r2, r2, :magenta)
   circle(x3, y3, r3, :green)
   circle(-x3, y3, r3, :green)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r0+r2, " r0+r2", :magenta, :left, :vcenter)
       point(0, r0 + 2r2 - r1, " r0+2r2-r1")
       point(x3, y3, "(x3,y3)")
       point(0, r0, "r0 ", :red, :right, :top, delta=-delta)
   else
       plot!(showaxis=false)
   end
end;

算額（その513）

算額（その513）

岩手県盛岡市 八幡宮 文化6年(1869)
一関市博物館>>和算に挑戦>>平成27年度出題問題>>平成27年度出題問題(2)[中級問題]&解答例
https://www.city.ichinoseki.iwate.jp/museum/wasan/h27/normal.html

大円内に 8 個の小円が入っている。大円の直径が 25 寸のとき，小円の直径を求めよ。

大円の直径と中心座標を R, (0, 0)
小円の直径と中心座標を r, (r, 2√3r), (0, √3r), (r, 0)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive

eq = r^2 + (2sqrt(Sym(3))r)^2 - (R - r)^2
solve(eq, r)[1] |> println

   R*(-1 + sqrt(13))/12

小円の半径は大円の半径の (√13 - 1)/12 倍である。

   R = 12.5;  r = 2.71412;  小円の直径 = 5.42823

大円の直径が 25 寸のとき小円の直径は約 5.42823 寸である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 25/2
   r = R*(√13 - 1)/12
   @printf("R = %g;  r = %g;  小円の直径 = %g\n", R, r, 2r)
   plot()
   circle(0, 0, R)
   circle4(r, 2√3r, r, :blue)
   circle(0, √3r, r, :blue)
   circle(0, -√3r, r, :blue)
   circle(r, 0, r, :blue)
   circle(-r, 0, r, :blue)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(r, 0, "r", :blue, :center, :bottom, delta=delta/2)
       point(0, √3r, " √3r", :blue, :left, :vcenter)
       point(r, 2√3r, " 小円:r\n (r,2√3r)", :blue, :center, :top, delta=-delta)
       point(R, 0, "R ", :red, :right, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;

算額（その512）

算額（その512）

兵庫県上郡町 上郡天満神社 文化3年(1806)
https://kamigori-town-museum.jp/cultural-asetts/559/
藤井康生: 上郡町天満神社算額見学記，和算，第93号，2002年1月31日
http://www.wasan.jp/kinkikaisi/wasan093.pdf

外円内に弦と矢および甲円，乙円が入っている。甲円と乙円の直径の和が 1尺4分，外円の直径と矢の差が 1尺4寸4分のとき矢を求めよ。

弦と y 軸の交点座標を (0, a) とおく。
外円の直径と中心座標を r0, (0, 0)
甲円の直径と中心座標を r1, (0, a + r1)
乙円の直径と中心座標を r2, (0, a + r2)
として，以下の連立方程式を解く。
矢（の長さ）は r0 - a である。

include("julia-source.txt");

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

eq1 = 2r1 + 2r2 - 104//10
eq2 = 2r0 - (r0 - a) - 144//10
eq3 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq4 = r1^2 + (a + r1)^2 - (r0 - r1)^2
eq5 = x2^2 + (a + r2)^2 - (r0 - r2)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (a, r0, r1, r2, x2))

   1-element Vector{NTuple{5, Sym}}:
    (63/20, 45/4, 18/5, 8/5, 42/5)

   a = 3.15;  r0 = 11.25;  r1 = 3.6; r2 = 1.6;  x2 = 8.4

矢 = res[1][2] - res[1][1]
矢 |> println

   81/10

矢は 8寸1分 である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r0, r1, r2, x2) = res[1]
   @printf("a = %g;  r0 = %g;  r1 = %g; r2 = %g;  x2 = %g\n",
       a, r0, r1, r2, x2)
   @printf("矢 = %g\n", r0 - a)
   plot()
   circle(0, 0, r0)
   circle(r1, a + r1, r1, :blue)
   circle(x2, a + r2, r2, :green)
   segment(-sqrt(r0^2 - a^2), a, sqrt(r0^2 - a^2), a, :brown)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(r1, a + r1, " 甲円:r1\n (r1,a+r1)", :blue, :left, :vcenter)
       point(x2, a + r2, "乙円:r2,(x2,a+r2) ", :black, :right, :vcenter)
       point(0, a, "a ", :brown, :right, :bottom, delta=delta/2)
       point(0, r0, "r0 ", :red, :right, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;

算額（その511）

算額（その511）

長野県長野市安茂里正覚院 久保寺観音堂 文政13年(1830)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

半円内に2本の斜線を引き，甲円 2 個と乙円を入れる。乙円の径が 1 寸のとき，甲円の径はいかほどか。

各円の半径，中心座標を以下のようにする。
外円: R, (0, 0)
甲円: r1, (x1, y1) および (x2, r1)
乙円: r3, (x3, y3)
外円の円周上にある2本の斜線の交点座標を (x, y)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms x1, y1, r1, x2, r3, x3, y3, x, R

r3 = 1//2
y = sqrt(R^2 - x^2)  # (x, y) は外円の円周上にある
eq1 = x1^2 + y1^2 - (R - r1)^2  # 右の甲円が外円に内接する
eq2 = x3^2 + y3^2 - (R - r3)^2  # 乙円が外円に内接する
eq3 = sqrt((R + x)^2 + y^2) + sqrt((R - x)^2 + y^2) - 2R - 2r1  # 鉤股弦に内接する円の直径
eq4 = (y3/x3)*(y/(x + R)) + 1
eq5 = numerator(distance(-R, 0, x, y, x2, r1) - r1^2) |> simplify
eq6 = numerator(distance(R, 0, x, y, x1, y1) - r1^2) |> simplify
eq7 = ((x - R)/2 - x3)^2 + (y/2 - y3)^2 - r3^2
eq8 = (y1/x1) * (y/(x - R)) + 1;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (x1, y1, r1, x2, x3, y3, x, R))

using Printf
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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (x1, y1, r1, x2, x3, y3, x, R) = u
   t = sqrt(R^2 - x^2)
   return [
       x1^2 + y1^2 - (R - r1)^2,  # eq1
       x3^2 + y3^2 - (R - 1/2)^2,  # eq2
       -2*R - 2*r1 + sqrt(R^2 - x^2 + (R - x)^2) + sqrt(R^2 - x^2 + (R + x)^2),  # eq3
       1 + y3*sqrt(R^2 - x^2)/(x3*(R + x)),  # eq4
       (-4*R^2*r1^2 + (R*r1 - R*t + r1*x - x2*t)^2 + (R^2 - R*x + R*x2 - r1*t - x*x2)^2)/(4*R^2),  # eq5
       (-4*R^2*r1^2 + (R*y1 - R*t - x*y1 + x1*t)^2 + (-R^2 - R*x + R*x1 + x*x1 + y1*t)^2)/(4*R^2),  # eq6
       (-y3 + t/2)^2 + (-R/2 + x/2 - x3)^2 - 1/4,  # eq7
       1 + y1*t/(x1*(-R + x)),  # eq8
   ]
end;

iniv = BigFloat[1.7, 4.2, 2.1, -3.5, -5.8, 2.4, -4.6, 6.4]
res = nls(H, ini=iniv)
println(res)

   (BigFloat[3.461538461538461538461538509169459561367959236632964526959234465030510845593508, 8.307692307692307692307692429999766197094115669634049393814154165178838427207838, 4.000000000000000000000000058695639811160390953981958412762403947086964336791015, -7.000000000000000000000000103675105612121327124114588925459978536567327493943071, -11.53846153846153846153846171778119379511106203957672690108302039508364569817559, 4.807692307692307692307692375277671136103845336879392455246078657573883943485054, -9.153846153846153846153846290689476553117442078831806765600936112963486841120747, 13.00000000000000000000000019241303418103091110455672021661459202029874280726291], true)

   r3 = 0.5;  x1 = 3.46154;  y1 = 8.30769; r1 = 4;  x2 = -7;  x3 = -11.5385;  y3 = 4.80769;  x = -9.15385;  y = 9.23077;  R = 13
   甲円の径 = 2r1 = 8

甲円の半径は 4 である。元の単位では甲円の直径は 8 寸である。

using Plots
function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1//2
   (x1, y1, r1, x2, x3, y3, x, R) = res[1]
   y = sqrt(R^2 - x^2)
   @printf("r3 = %g;  x1 = %g;  y1 = %g; r1 = %g;  x2 = %g;  x3 = %g;  y3 = %g;  x = %g;  y = %g;  R = %g\n",
       r3, x1, y1, r1, x2, x3, y3, x, y, R)
   @printf("甲円の径 = 2r1 = %g\n", 2r1)
   plot([R, x, -R, R], [0, y, 0, 0], color=:brown, lw=0.5)
   circle(0, 0, R, beginangle=0, endangle=180)
   circle(x1, y1, r1, :blue)
   circle(x2, r1, r1, :green)
   circle(x3, y3, r3, :magenta)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x1, y1, " 甲円:r1\n (x1,y1)", :blue, :left, :vcenter)
       point(x2, r1, " 甲円:r1\n (x2,r1)", :blue, :left, :top, delta=-delta)
       point(x3, y3, "  乙円:r3,(x3,y3)", :magenta, :left, :vcenter)
       point(x, y, "(x,y) ", :brown, :right, :vcenter)
   else
       plot!(showaxis=false)
   end
end;