算額（その772）

算額（その772）

宮城県白石市小原 小原温泉薬師堂 大正5年(1916)
徳竹亜紀子，谷垣美保：宮城県白石市小原地区の算額調査，仙台高等専門学校名取キャンパス研究紀要，第57号，2021.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2021/04/No57_2.pdf

外円の中に大円，中円，小円，甲円，乙円が入っている。中円，小円，甲円の直径がそれぞれ 9 寸，6 寸，3 寸のとき，乙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, R - r1)
中円の半径と中心座標を r2, (x2, y2)
小円の半径と中心座標を r3, (x3, y3)
甲円の半径と中心座標を r4, (x4, y4)
乙円の半径と中心座標を r5, (x5, y5)
とおき，以下の連立方程式を解く。

算額にある図とは大幅に異なる結果になる。そもそも「大円」は「中円」はおろか，「小円」よりも小さい。

include("julia-source.txt");

using SymPy
@syms R, r1, r2, x2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq3 = x4^2 + (R - r1 - y4)^2 - (r1 + r4)^2
eq4 = x5^2 + (R - r1 - y5)^2 - (r1 + r5)^2
eq5 = x2^2 + y2^2 - (R - r2)^2
eq6 = x3^2 + y3^2 - (R - r3)^2
eq7 = x4^2 + y4^2 - (R - r4)^2
eq8 = x5^2 + y5^2 - (R - r5)^2
eq9 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq10 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq11 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2;

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

function H(u)
   (R, r1, x2, y2, x3, y3, x4, y4, r5, x5, y5) = u
   return [
       x2^2 - (r1 + r2)^2 + (R - r1 - y2)^2,  # eq1
       x3^2 - (r1 + r3)^2 + (R - r1 - y3)^2,  # eq2
       x4^2 - (r1 + r4)^2 + (R - r1 - y4)^2,  # eq3
       x5^2 - (r1 + r5)^2 + (R - r1 - y5)^2,  # eq4
       x2^2 + y2^2 - (R - r2)^2,  # eq5
       x3^2 + y3^2 - (R - r3)^2,  # eq6
       x4^2 + y4^2 - (R - r4)^2,  # eq7
       x5^2 + y5^2 - (R - r5)^2,  # eq8
       -(r2 + r3)^2 + (x2 - x3)^2 + (y2 - y3)^2,  # eq9
       -(r2 + r4)^2 + (x2 - x4)^2 + (y2 - y4)^2,  # eq10
       -(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq11
   ]
end;

(r2, r3, r4) = (9, 6, 3) .// 2
iniv = BigFloat[7.5, 2.43, 2.7, -1.31, -4.2, 1.62, 3.9, 4.56, 1.0, -3.4, 5.54]
res = nls(H, ini=iniv)

   ([7.5031904642363045, 2.4327679290250375, 2.7, -1.3149726097831362, -4.2, 1.6244150815566774, 3.9, 4.563802772896491, 1.0, -3.4, 5.543598670009763], true)

乙円の半径は 1 寸である（直径は 2 寸）。

その他のパラメータは以下のとおりである。

   R = 7.50319;  r1 = 2.43277;  x2 = 2.7;  y2 = -1.31497;  x3 = -4.2;  y3 = 1.62442;  x4 = 3.9;  y4 = 4.5638;  r5 = 1;  x5 = -3.4;  y5 =5.5436

function draw(more=false)
   pyplot(size=(400, 400), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, r4) = (9, 6, 3) .// 2
   (R, r1, x2, y2, x3, y3, x4, y4, r5, x5, y5) = res[1]
   @printf("r2 = %g;  r3 = %g;  r4 = %g のとき，乙円の直径は %g\n", r2, r3, r4, 2r5)
   @printf("R = %g;  r1 = %g;  x2 = %g;  y2 = %g;  x3 = %g;  y3 = %g;  x4 = %g;  y4 = %g;  r5 = %g;  x5 = %g;  y5 =%g\n", R, r1, x2, y2, x3, y3, x4, y4, r5, x5, y5)
   plot()
   circle(0, 0, R)
   circle(0, R - r1, r1, :blue)
   circle(x2, y2, r2, :green)
   circle(x3, y3, r3, :magenta)
   circle(x4, y4, r4, :purple)
   circle(x5, y5, r5, :orange)
   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, R - r1, "大円:r1\n(0,R-r1)", :blue, :center, :bottom, delta=delta)
       point(x2, y2, "中円:r2\n(x2,y2)", :green, :center, :bottom, delta=delta)
       point(x3, y3, "小円:r3\n(x3,y3)", :magenta, :center, :bottom, delta=delta)
       point(x4, y4, "甲円:r4,(x4,y4)", :purple, :left, delta=-delta)
       point(x5, y5, "乙円:r5,(x5,y5)", :black, :right, delta=-delta)
       point(R, 0, " R", :red, :left, :bottom, delta=delta)
   end
end;

算額（その771）

算額（その771）

神壁算法 天明5年乙巳 關流四傅藤田貞資六弟子 筑後州久留米 佐田甚助方之
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

和算問題あれこれ2 令和６年４月の問題-No.３
https://gunmawasan.web.fc2.com/k-n-mondai.html

円弧の弦の上に，底辺が弦の上で頂点が円弧上にある甲，乙，丙の正三角形が載っている。一辺の長さがそれぞれ 206 寸 1 分，229 寸，183 寸 2 分のとき，弦の長さはいかほどか。

正三角形の底辺の 2 つの x 座標，頂点の x 座標を，左から（丙，乙，甲の順で），a, b, c, d, e, f, g とする。
弦と y 軸の交点の y 座標を h とする。
頂点の y 座標 は h + √3(b - a)，h + √3(d - c)，h + √3(f - e) とする。
A = g - e，B = e - c，C = c - a とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     c::positive, d::positive,
     e::positive, f::positive,
     g::negative, h::positive,
     R::positive
@syms a, b, c, d, e, f, g, h, R, A, B, C
s3 = sqrt(Sym(3))
eq1 = g - e - A
eq2 = e - c - B
eq3 = c - a - C
eq4 = b - (a + c)/2
eq5 = d - (c + e)/2
eq6 = f - (e + g)/2
eq7 = (h + (b - a)s3)^2 + b^2 - R^2
eq8 = (h + (d - c)s3)^2 + d^2 - R^2
eq9 = (h + (f - e)s3)^2 + f^2 - R^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9],
   (R, h, a, b, c, d, e, f, g))

   2-element Vector{NTuple{9, Sym{PyCall.PyObject}}}:
    (-sqrt(3)*sqrt((A^2 - A*B + B^2)*(B^2 - B*C + C^2)*(A^2 + A*B - A*C + B^2 + B*C + C^2))/(3*(A*C - B^2)), sqrt(3)*(-A^2*B - A^2*C - A*B*C - A*C^2 + B^3 - B*C^2)/(6*(A*C - B^2)), (A^2*B - A^2*C - A*B*C - A*C^2 + B^3 + 2*B^2*C - B*C^2)/(2*(A*C - B^2)), (A^2*B - A^2*C - A*B*C + B^3 + B^2*C - B*C^2)/(2*(A*C - B^2)), (B - C)*(A^2 - A*C + B^2 + B*C)/(2*(A*C - B^2)), (A - C)*(A*B - A*C + B*C)/(2*(A*C - B^2)), (A - B)*(A*B - A*C + B^2 + C^2)/(2*(A*C - B^2)), (A^2*B - A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)), (A^2*B + A^2*C - 2*A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)))
    (sqrt(3)*sqrt((A^2 - A*B + B^2)*(B^2 - B*C + C^2)*(A^2 + A*B - A*C + B^2 + B*C + C^2))/(3*(A*C - B^2)), sqrt(3)*(-A^2*B - A^2*C - A*B*C - A*C^2 + B^3 - B*C^2)/(6*(A*C - B^2)), (A^2*B - A^2*C - A*B*C - A*C^2 + B^3 + 2*B^2*C - B*C^2)/(2*(A*C - B^2)), (A^2*B - A^2*C - A*B*C + B^3 + B^2*C - B*C^2)/(2*(A*C - B^2)), (B - C)*(A^2 - A*C + B^2 + B*C)/(2*(A*C - B^2)), (A - C)*(A*B - A*C + B*C)/(2*(A*C - B^2)), (A - B)*(A*B - A*C + B^2 + C^2)/(2*(A*C - B^2)), (A^2*B - A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)), (A^2*B + A^2*C - 2*A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)))

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

弦の長さは 2sqrt(R^2 - h^2) である。
正三角形の一辺の長さがそれぞれ 206 寸 1 分，229 寸，183 寸 2 分のとき，弦の長さは 1029.80004005632 寸である。

2sqrt(res[2][1]^2 - res[2][2]^2)(A => 206.1, B => 229.0, C => 183.2) |> println

   1029.80004005632

その他のパラメータは以下のとおりである。
R = 764.582;  h = 565.211;  a = -337.775;  b = -246.175;  c = -154.575;  d = -40.075;  e = 74.425;  f = 177.475;  g = 280.525

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (A, B, C) = (206.1, 229.0, 183.2)
   (R, h, a, b, c, d, e, f, g) = (-sqrt(3)*sqrt((A^2 - A*B + B^2)*(B^2 - B*C + C^2)*(A^2 + A*B - A*C + B^2 + B*C + C^2))/(3*(A*C - B^2)), sqrt(3)*(-A^2*B - A^2*C - A*B*C - A*C^2 + B^3 - B*C^2)/(6*(A*C - B^2)), (A^2*B - A^2*C - A*B*C - A*C^2 + B^3 + 2*B^2*C - B*C^2)/(2*(A*C - B^2)), (A^2*B - A^2*C - A*B*C + B^3 + B^2*C - B*C^2)/(2*(A*C - B^2)), (B - C)*(A^2 - A*C + B^2 + B*C)/(2*(A*C - B^2)), (A - C)*(A*B - A*C + B*C)/(2*(A*C - B^2)), (A - B)*(A*B - A*C + B^2 + C^2)/(2*(A*C - B^2)), (A^2*B - A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)), (A^2*B + A^2*C - 2*A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)))
   @printf("R = %g;  h = %g;  a = %g;  b = %g;  c = %g;  d = %g;  e = %g;  f = %g;  g = %g\n", R, h, a, b, c, d, e, f, g)
   s3 = √3
   plot([a, b, c, d, e, f, g], h .+ [0, (b - a)s3, 0, (d - c)s3, 0, (f - e)s3, 0], color=:blue, lw=0.5)
   θ = asind(h/R)
   circle(0, 0, R, beginangle=θ, endangle=180 - θ, n=1000)
   x = R*cosd(θ)
   segment(-x, h, x, h, :red)
   plot!([-x, 0, x], [h, 0, h], color=:gray50, lw=0.5)
   @printf("弦の長さ = %g\n", 2R*cosd(θ))
   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.([a, b, c, d, e, f, g], h, Char["abcdefg"...], :blue, :center, delta=-delta)
       point(x, h, "(x,h) ", :red, :right, delta=-delta)
       point(x/2+4delta, h/2, "R", :gray50, :center, :vcenter, mark=false)
       point.([b, d, f], h + 75, ["丙三角", "乙三角", "甲三角"], :blue, :center, delta=-delta, mark=false)
   end
end;

算額（その770）

算額（その770）

神壁算法 天明8年戊申 二月 關流藤田貞資門人 細川能登守家士 渡邉金八源耀
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

正五角形に5本の斜線を入れ，内部に小さな正五角形を作り，できた 6 つの区画にそれぞれ同じ直径を持つ円を内接させる。
正五角形の一辺の長さが 533 寸 6 分のとき，等円の直径はいかほどか。

正五角形の一辺の長さを L，正五角形が内接する円の半径を R とする。
内側の正五角形が内接する円の半径を c，その一辺の長さを l とする。
等円の半径を r1 とする。
斜線と正五角形の一辺で作られる三角形 BAR において第二余弦定理を使う。
AR, AB, BR の長さを AR, AB, BR，また，A, B, R の座標を (x, y), (x0, y0), (0, R) とする。
ポイントは，斜線は外側の五角形の一辺と平行ではないということである。

include("julia-source.txt");

using SymPy
@syms L::positive, R::positive,
     x0::positive, y0::positive,
     x1::positive, y1::positive,
     x::negative, y::positive,
     AR::positive, l::positive, r1::positive
d18 = Sym(18)

R = L/2sind(2d18)
x0 = R*cosd(d18)
y0 = R*sind(d18)
c = sqrt(x^2 + y^2)

eq1 = (x0 - x)^2 + (y0 - y)^2 - (AR + l)^2
eq2 = x^2 + (R - y)^2 - AR^2
eq3 = AR^2 + (AR + l)^2 - 2AR*(AR + l)*cosd(4d18) - L^2
eq4 = dist(x0, y0, R*cosd(-3d18), R*sind(-3d18), x1, y1) - r1^2
eq5 = dist(x0, y0, x, y, x1, y1) - r1^2
eq6 = dist(x0, y0, x, y, 0, 0) - r1^2
eq7 = dist(c*cosd(atand(y, x) - 4d18), c*sind(atand(y, x) - 4d18), R*cosd(-3d18), R*sind(-3d18), x1, y1) - r1^2;

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

function H(par)
    (AR, l, x, y, x1, y1, r1) = par
    d18 = 18
    R = L/2sind(2d18)
    x0 = R*cosd(d18)
    y0 = R*sind(d18)
    c = sqrt(x^2 + y^2)
    return [
        (x0 - x)^2 + (y0 - y)^2 - (AR + l)^2,
        x^2 + (R - y)^2 - AR^2,
        AR^2 + (AR + l)^2 - 2AR*(AR + l)*cosd(4d18) - L^2,
        dist(x0, y0, R*cosd(-3d18), R*sind(-3d18), x1, y1) - r1^2,
        dist(x0, y0, x, y, x1, y1) - r1^2,
        dist(x0, y0, x, y, 0, 0) - r1^2,
        dist(c*cosd(atand(y, x) - 4d18), c*sind(atand(y, x) - 4d18), R*cosd(-3d18), R*sind(-3d18), x1, y1) - r1^2
    ]
end;

L = 5336//10
iniv = BigFloat[347.1, 179.5, -94.4, 119.9, 259.5, 10.0, 123.5]
res = nls(H, ini=iniv)

   ([347.0551765561487, 179.4560303257959, -94.42748155686841, 119.94507841468496, 259.51771893633804, 10.023199995813941, 123.50001782997607], true)

正五角形の一辺の長さが 533 寸 6 分のとき，等円の直径は 247 寸である。

その他のパラメータは以下のとおりである。

AR = 347.055;  x = -94.4275;  y = 119.945;  x1 = 259.518;  y1 = 10.0232;  r1 = 123.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   L = 5336//10
   R = L/2sind(2*18)
   x0 = R*cosd(18)
   y0 = R*sind(18)
   (AR, l, x, y, x1, y1, r1) = res[1]
   @printf("等円の直径 = %g\n", 2r1)
   @printf("AR = %g;  x = %g;  y = %g;  x1 = %g;  y1 = %g;  r1 = %g\n", AR, x, y, x1, y1, r1)
   r = sqrt(x^2 + y^2)
   θ = atand(y, x)
   ix = [r*cosd(72i + θ) for i in 0:5]
   iy = [r*sind(72i + θ) for i in 0:5]
   ox = [R*cosd(72i + 18) for i in 0:5]
   oy = [R*sind(72i + 18) for i in 0:5]
   plot()
   circle(0, 0, R, :gray80)
   circle(0, 0, r, :gray90)
   circle(0, 0, r1)
   rotate(x1, y1, r1, angle=72)
   for i in 1:5
       segment(ox[i], oy[i], ix[i], iy[i], :blue)
       segment(ox[i], oy[i], ox[i + 1], oy[i + 1], :blue)
   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(x, y, "A:(x,y) ", :blue, :right, :vcenter)
       point(x0, y0, " B:(x0,y0)", :blue, :left, :vcenter)
       point(x1, y1, "等円:r1\n(x1,y1)", :red, :center, :bottom, delta=delta)
       point(0, R, " R:(0,R)", :black, :left, :bottom, delta=delta)
       point(r, 0, " r", :gray50, :left, :bottom, delta=delta)
   end
end;