算額（その876）

算額（その876）

新潟県長岡市 蒼柴神社 亨和元年(1801)3月
http://www.wasan.jp/niigata/aoshi.html

涌田和芳，外川一仁: 長岡蒼柴神社の算額，長岡工業高等専門学校研究紀要，第42巻，第2号(2006)
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_41-45/vol_42_2/42_2_1wakuta.pdf

外円内に，弦と 2 本の斜線を引き，区分された領域に甲円と乙円を 2 個ずつ入れる。大円の直径が 521 寸で，弦と矢の長さの差を最大にするとき，乙円の直径を求めよ。
注：矢（し）とは，弓形の孤の中点から弦におろした垂線
外円の半径と中心座標を R, (0, 0)
弦，矢の長さおよびその差を（そのまま）弦，矢，弦矢差
甲円の半径と中心座標を r1, (0, R - r1), (0, r - 矢 + r1)
乙円の半径と中心座標を r2, (x2, y2)
とおき以下の連立方程式を解く。

まず，弦と矢の長さの差が最大になるときの矢と弦を決定する。

include("julia-source.txt");

using SymPy
@syms R::positive, 弦::positive, 矢::positive, x::positive,
     r1::positive, r2::positive, x2::positive, y2::positive, d
R = 521//2
弦 = 2sqrt(R^2 - (R - 矢)^2)
弦矢差 = 弦 - 矢;

pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(弦矢差, xlims=(9, R), xlabel="矢", ylabel="弦 - 矢")

矢が 150 前後のときに弦矢差が最大になることがわかる。

最大値を取るときの矢の値を求めるには，導関数を求め，導関数が 0 になるときの値を求めればよい。

g = diff(弦矢差, 矢);
g |> println

   2*(521/2 - 矢)/sqrt(271441/4 - (521/2 - 矢)^2) - 1

mx_a = solve(g, 矢)[1]
mx_a |> factor |> println
mx_a.evalf() |> println

   -521*(-5 + sqrt(5))/10
   144.000858372261

弦矢差(矢 => mx_a.evalf()) |> println

   321.995708138695

弦(矢 => mx_a.evalf()) |> println

   465.996566510956

矢が 521/2 - 521*sqrt(5)/10 = 144.000858372261 のとき，弦と矢の差が最大値 321.995708138695 になる。ちなみにそのときの弦は 465.996566510956 である。

465.996566510956 - 144.000858372261, 321.995708138695

   (321.99570813869497, 321.995708138695)

弦と矢の差が最大となるときの矢，弦が決まったので，その状況における甲円と乙円の半径を求める。

using SymPy
@syms R::positive, 弦::positive, 矢::positive, x::positive,
     r1::positive, r2::positive, x2::positive, y2::positive, d
矢 = (5 - sqrt(Sym(5)))R/5
弦 = 2sqrt(R^2 - (R - 矢)^2)
eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = dist(x, sqrt(R^2 - x^2), -弦/2, R - 矢, 0, R - r1) - r1^2
eq3 = dist(x, sqrt(R^2 - x^2), -弦/2, R - 矢, 0, R - 矢 + r1) - r1^2
eq4 = dist(x, sqrt(R^2 - x^2), -弦/2, R - 矢, x2, y2) - r2^2
eq5 = dist(-x, sqrt(R^2 - x^2), 弦/2, R - 矢, x2, y2) - r2^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)
   (r1, r2, x2, y2, x) = u
   return [
x2^2 + y2^2 - (R - r2)^2,  # eq1
-r1^2 + (-x - (-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))*(-x*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (R - r1 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2 + (R - r1 - sqrt(R^2 - x^2) - (-x*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (R - r1 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2,  # eq2
-r1^2 + (-x - (-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))*(-x*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R + r1 - sqrt(R^2 - x^2)))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2 + (-R*(5 - sqrt(5))/5 + R + r1 - sqrt(R^2 - x^2) - (-x*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R + r1 - sqrt(R^2 - x^2)))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2,  # eq3
-r2^2 + (-x + x2 - (-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))*((-x + x2)*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (y2 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2 + (y2 - sqrt(R^2 - x^2) - ((-x + x2)*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (y2 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2,  # eq4
-r2^2 + (x + x2 - (x + sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))*((x + x2)*(x + sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (y2 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))/((x + sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2 + (y2 - sqrt(R^2 - x^2) - ((x + x2)*(x + sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (y2 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))/((x + sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2,  # eq5
   ]
end;

R = 521//2
iniv = BigFloat[35, 37, 122, 188, 127]
res = nls(H, ini=iniv)

   ([35.179523802100135, 36.00021459306524, 121.93467698217249, 188.49957081386952, 126.65410684822777], true)

外円の直径が 521 寸のとき，乙円の直径は 72 寸有奇である。

その他のパラメータは以下のとおりである。
r1 = 35.1795;  r2 = 36.0002;  x2 = 121.935;  y2 = 188.5;  x = 126.654

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 521//2
   矢 = (5 - √5)R/5
   弦 = 2sqrt(R^2 - (R - 矢)^2)
   (r1, r2, x2, y2, x) = res[1]
   @printf("弦 = %g;  矢 = %g;  弦と矢の差 = %g\n", 弦, 矢, 弦 - 矢)
   @printf("乙円の直径 = %.15g\n", 2r2)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  x = %g\n", r1, r2, x2, y2, x)
   plot()
   circle(0, 0, R)
   segment(-弦/2, R - 矢, 弦/2, R - 矢, :blue)
   circle2(x2, y2, r2, :green)
   circle(0, R - r1, r1, :magenta)
   circle(0, R - 矢 + r1, r1, :magenta)
   y = sqrt(R^2 - x^2)
   segment(x, y, -弦/2, R - 矢)
   segment(-x, y, 弦/2, R - 矢)
   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 - 矢, "R-矢", :blue, :center, delta=-delta/2)
       point(x, y, "(x,y)", :green, :left, :bottom, delta=delta/2)
       point(弦/2, R - 矢, "弦/2 ", :blue, :right, delta=-delta/2)
       point(0, R - r1, "甲円:r1,(0,R-r1)", :black, :center, :bottom, delta=delta/2)
       point(0, R - 矢/2, "R-矢/2   ", :black, :right, :vcenter)
       point(0, R - 矢 + r1, "甲円:r1,(0,R-矢+r1)", :black, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :black, :center, delta=-delta/2)
       point(0, R, "R", :red, :center, :bottom, delta=delta/2)
 end
end;

算額（その875）

算額（その875）

新潟県長岡市 蒼柴神社 亨和元年(1801)3月
http://www.wasan.jp/niigata/aoshi.html

涌田和芳，外川一仁: 長岡蒼柴神社の算額，長岡工業高等専門学校研究紀要，第42巻，第2号(2006)
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_41-45/vol_42_2/42_2_1wakuta.pdf

角錐台の中に大球 1 個，小球 4 個が入っている。上底，下底の正方形の一辺の長さがそれぞれ 2 寸，6 寸のとき，大球の直径はいかほどか。

角錐と角錐台の高さを h, h2
上底，下底の正方形の一辺の長さを 2a, 2b
大球の半径と中心座標を r1, (0, 0, h2 - r1)
小球の半径と中心座標を r2, (r2, r2, r2)
と置く。

eq1: x 軸 に対して 45° 方向の負の方向から投影すると，大円と小円は互いに外接している。

eq2, eq3, eq4: x 軸の正の方向から y-z 平面に投影された図形は下図のように大円，小円は角錐台の面に内接している。

これらの連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive, h::positive, h2::positive, d
(a, b) = (1, 3)
eq1 = 2r2^2 + (h2 - r1 - r2)^2 - (r1 + r2)^2 |> expand
eq2 = (b + h - 2r2)^2 - (h^2 + b^2)  # 3 + h - sqrt(h^2 + 9) - 2r2
eq3 = b/h - a/(h - h2)
eq4 = numerator(apart(dist(0, h, b, 0, 0, h2 - r1) - r1^2, d))
res = solve([eq1, eq2, eq3, eq4], (r1, r2, h, h2));
res[1]  # 1 of 3

   (3/2, 6/5, 36/5, 24/5)

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

上底，下底の正方形の一辺の長さがそれぞれ 2 寸，6 寸のとき，大球の直径は 3 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, h, h2) = (3/2, 6/5, 36/5, 24/5)
   plot([b, 0, -b, b], [0, h, 0, 0], color=:blue, lw=0.5)
   circle(0, h2-r1,r1, :green)
   circle(r2, r2, r2)
   circle(-r2, r2, r2, :lightpink)
   circle(0, r2, r2, :lightpink)
   segment(-a, h2, a, h2, :blue)
   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, h2-r1, "大球:r1,(0,0,h2-r1)", :green, :center, delta=-delta/2)
       point(r2, r2, "小球:r2,(r1,r1,r1)", :red, :center, delta=-delta/2)
       point(0, h2, "h2 ", :blue, :right, :bottom, delta=delta/2)
       point(a, h2, "(a,h2) ", :blue, :right, :bottom, delta=delta/2)
       point(b, 0, "b ", :blue, :right, :bottom, delta=delta/2)
       point(0, h, "h ", :blue, :right, :bottom, delta=delta/2)
   end
end;

function draw2(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, h, h2) = (3/2, 6/5, 36/5, 24/5)
   plot(√2 .* [b, a, -a, -b, b], [0, h2, h2, 0, 0], color=:blue, lw=0.5)
   circle(√2r2, r2, r2)
   circle(-√2r2, r2, r2, :lightpink)
   circle(0, r2, r2, :lightpink)
   circle(0, h2 - r1, r1, :green)
   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, h2-r1, "大球:r1,(0,0,h2-r1)", :green, :center, delta=-delta/2)
       point(√2r2, r2, "小球:r2,(√2r1,r1,r1)", :red, :center, delta=-delta/2)
       point(0, h2, "h2 ", :blue, :right, :bottom, delta=delta/2)
       point(√2a, h2, "(√2a,h2) ", :blue, :right, :bottom, delta=delta/2)
       point(√2b, 0, "b ", :blue, :right, :bottom, delta=delta/2)
       #point(0, h, "h ", :blue, :right, :bottom, delta=delta/2)
   end
end;

算額（その874）

算額（その874）

秋田県大仙市角間川町 角間川前寺

矢野充志：いろいろな算額，平成30年度 微分積分Ⅰ (2S・2I・2C)，2019/02/25.
https://www.libe.nara-k.ac.jp/~yano/biseki1_2018/meisatsu2.pdf

正方形の内部に斜線，四分円，半円，大円，小円を入れる。小円の直径が与えられたとき，大円の直径を求めよ。

正方形の一辺の長さを 2r1
四分円の半径と中心座標を 2r1, (0, 0)
半円の半径と中心座標を r1, (2r1, r1)
大円の半径と中心座標を r2, (r2, r2)
小円の半径と中心座標を r3, (2r1 - r3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, y3::positive, x::positive
eq1 = (2r1 - r3)^2 + y3^2 - (2r1 + r3)^2
eq2 = r3^2 + (y3 - r1)^2 - (r1 - r3)^2
eq3 = dist(0, 2r1, x, 0, r2, r2) - r2^2
eq4 = dist(0, 2r1, x, 0, 2r1, r1) - r1^2
res = solve([eq1, eq2, eq3, eq4], (r1, r2, y3, x))

   2-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (25*r3/8, 25*r3/16, 5*r3, 75*r3/16)
    (25*r3/8, 75*r3/8, 5*r3, 75*r3/16)

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

大円の半径は，小円の半径の 25/16 倍である。
小円の直径が 16 のとき，大円の直径は 25 である

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

   r3 = 8;  r1 = 25;  r2 = 12.5;  y3 = 40;  x= 37.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 8
   (r1, r2, y3, x) = r3 .* (25/8, 25/16, 5, 75/16)
   @printf("小円の直径が %g のとき，大円の直径は %g である\n", 2r3, 2r2)
   @printf("r3 = %g;  r1 = %g;  r2 = %g;  y3 = %g;  x= %g\n", r3, r1, r2, y3, x)
   plot([0, 2r1, 2r1, 0, 0], [0, 0, 2r1, 2r1, 0], color=:blue, lw=0.5)
   circle(0, 0, 2r1, beginangle=0, endangle=90)
   circle(2r1, r1, r1, :green, beginangle=90, endangle=270)
   circle(r2, r2, r2, :magenta)
   circle(2r1 - r3, y3, r3, :orange)
   segment(0, 2r1, x, 0, :brown)
   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(2r1, 0, "2r1 ", :blue, :right, :bottom, delta=delta/2)
       point(x, 0, " x", :blue, :left, :bottom, delta=delta/2)
       point(0, 2r1, " 2r1", :blue, :left, :bottom, delta=delta/2)
       point(r2, r2, "大円:r2,(r2,r2)", :magenta, :center, delta=-delta/2)
       point(2r1 - r3, y3, "小円:r3,(2r1-r3,y3)", :orange, :center, delta=-delta/2)
       point(2r1, r1, "半円:r1\n(2r1,r1)", :green, :right, :vcenter)
   end
end;