算額（その420）

算額（その420）

茨城県笠間市福原 吾國山上頂神（吾国山上頂神前） 文化7年(1810)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

長径 8 寸，短径 6 寸の菱形内に，甲円 2 個，乙円 2 個が入っている。甲円の直径が 4 寸のとき，乙円の直径を求めよ。

菱形の中心を原点とする。
甲円の半径と中心座標を r1, (x1, 0)
乙円の半径と中心座標を r2, (0, y2)
とし，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x1::positive, y2::positive;

eq1 = x1^2 + y2^2 - (r1 + r2)^2
eq2 = distance(0, 3, 4, 0, x1, 0) - r1^2
eq3 = distance(0, 3, 4, 0, 0, y2) - r2^2;
res = solve([eq1, eq2, eq3], (r2, x1, y2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (4*(-16*sqrt(15)*r1^(5/2) - 900*sqrt(15)*sqrt(r1) + 16*r1^3 + 60*r1^2 + 900*r1 + 3375)/(9*(4*r1^2 + 225)), -(5*r1 - 12)/3, 20*sqrt(15)*sqrt(r1)/9 - 20*r1/9 - 16/3)
    (4*(4*r1 - 15)/9, -(5*r1 - 12)/3, 20*r1/9 - 16/3)

最初の組が適解である。

r1 = 4/2
(4*(-16*sqrt(15)*r1^(5/2) - 900*sqrt(15)*sqrt(r1) + 16*r1^3 + 60*r1^2 + 900*r1 + 3375)/(9*(4*r1^2 + 225)), -(5*r1 - 12)/3, 20*sqrt(15)*sqrt(r1)/9 - 20*r1/9 - 16/3)

   (0.48493231101926837, 0.6666666666666666, 2.3938346112259135)

r2 = 0.484932;  x1 = 0.666667;  y2 = 2.39383

乙円の直径は 0.969865 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 4//2
   (r2, x1, y2) = (4*(-16*sqrt(15)*r1^(5/2) - 900*sqrt(15)*sqrt(r1) + 16*r1^3 + 60*r1^2 + 900*r1 + 3375)/(9*(4*r1^2 + 225)), -(5*r1 - 12)/3, 20*sqrt(15)*sqrt(r1)/9 - 20*r1/9 - 16/3)
   @printf("r2 = %g;  x1 = %g;  y2 = %g\n", r2, x1, y2)
   @printf("乙円の直径 = %g\n", 2r2)
   plot([4, 0, -4, 0, 4], [0, 3, 0, -3, 0], color=:black, lw=0.5)
   circle(x1, 0, r1)
   circle(-x1, 0, r1)
   circle(0, y2, r2, :blue)
   circle(0, -y2, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, 0, " x1", :red)
       point(0, y2, " y2", :blue)
       point(4, 0, " 4", :black, :left, :bottom)
       point(0, 3, " 3", :black, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その419）

算額（その419）

東京都江東区木場（深川） 大丸屋内稲荷絵馬堂 文化7年(1810)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

弓形の中に正三角形 1 個，大円 2 個，中円 4 個，小円 4 個が入っている。大円の直径が 1 尺のとき，小円の直径を求めよ。

弓形を構成する外円の中心を原点とする。
弓形を構成する弦と外円の中心の距離を a とする（図参照）
外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (x1, a + r1)
中円の半径と中心座標を r2, (x21, a + r2), (x22, y22)
小円の半径と中心座標を r3, (x31, a + r3), (x32, y32)
として，以下の連立方程式の数値解を求める。

未知数は r0, a, x1, r2, x21, x22, y22, r3, x31, x32, y32 の11個であるが，条件式（方程式）は 12 個になった。
以下の eq4 〜 eq9 のうちどれか一つを使わないようにする（eq9 を除いた）。

include("julia-source.txt");

using SymPy

@syms r0::positive, a::positive, r1::positive, x1::positive,
     r2::positive, x21::positive, x22::positive, y22::positive,
     r3::positive, x31::positive, x32::positive, y32::positive;

eq1 = x1^2 + (a + r1)^2 - (r0 - r1)^2
eq2 = x21^2 + (a + r2)^2 - (r0 - r2)^2
eq3 = x22^2 + y22^2 - (r0 - r2)^2
eq4 = (x21 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq5 = (x1 - x22)^2 + (y22 - a - r1)^2 - (r1 + r2)^2
eq6 = (x31 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq7 = (x1 - x32)^2 + (y32 - a - r1)^2 - (r1 + r3)^2
eq8 = (x21 - x31)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq9 = (x22 - x32)^2 + (y22 - y32)^2 - (r2 + r3)^2
eq10 = distance(0, r0, (r0 - a)/sqrt(Sym(3)), a, x1, a + r1) - r1^2
eq11 = distance(0, r0, (r0 - a)/sqrt(Sym(3)), a, x22, y22) - r2^2
eq12 = distance(0, r0, (r0 - a)/sqrt(Sym(3)), a, x32, y32) - r3^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 v, r.f_converged
end;

function H(u)
   (r0, a, x1, r2, x21, x22, y22, r3, x31, x32, y32) = u
   return [
       x1^2 + (a + r1)^2 - (r0 - r1)^2,  # eq1
       x21^2 + (a + r2)^2 - (r0 - r2)^2,  # eq2
       x22^2 + y22^2 - (r0 - r2)^2,  # eq3
       (r1 - r2)^2 - (r1 + r2)^2 + (-x1 + x21)^2,  # eq4
       -(r1 + r2)^2 + (x1 - x22)^2 + (-a - r1 + y22)^2,  # eq5
       (r1 - r3)^2 - (r1 + r3)^2 + (-x1 + x31)^2,  # eq6
       -(r1 + r3)^2 + (x1 - x32)^2 + (-a - r1 + y32)^2,  # eq7
       (r2 - r3)^2 - (r2 + r3)^2 + (x21 - x31)^2,  # eq8
       #-(r2 + r3)^2 + (x22 - x32)^2 + (y22 - y32)^2,  # eq9
       -r1^2 + (a/4 - r0/4 + r1/4 + sqrt(3)*x1/4)^2 + (sqrt(3)*a/4 - sqrt(3)*r0/4 + sqrt(3)*r1/4 + 3*x1/4)^2,  # eq10
       -r2^2 + (-r0/4 + sqrt(3)*x22/4 + y22/4)^2 + (-sqrt(3)*r0/4 + 3*x22/4 + sqrt(3)*y22/4)^2,  # eq11
       -r3^2 + (-r0/4 + sqrt(3)*x32/4 + y32/4)^2 + (-sqrt(3)*r0/4 + 3*x32/4 + sqrt(3)*y32/4)^2,  # eq12
   ]
end;
r1 = 10//2
iniv = [big"100.0", 60, 32, 8, 53, 15, 100, 3, 47, 12, 95] ./ 3.2
iniv = [big"31.25", 18.75, 10.0, 2.5, 16.5625, 4.6875, 31.25, 0.9375, 14.6875, 3.75, 29.6875]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[25.49038105676657998080848135617869247549515697170824483702794809282177979870005, 12.74519052838329013851885544688169255263494103177866816560880000647679820897323, 10.24519052838328991096467576691497724140415363411739473165687289639408890927506, 2.470438416976302806278212382647506529008215409707534746382983461692078901157198, 17.27432762617317838252967445536201528319013903197721182970015560305202172736248, 4.539957388152644446167925002781799424476871763565147607613821577005623652425758, 22.56782103024110685639605573917644547722212244239117758536157640732734540923325, 0.851900255445560589676464780553194462776632582898376305125453038963129343031564, 14.37290237744824272406559296990917376944931384289685950439663092059990005842643, 4.58897484763492825949410623312361446943111492490923485552995147479476165882802, 19.24584397689835787931127610394921161788065249923856835686360465925916856690406], true)

r0 = 25.4904;  a = 12.7452;  x1 = 10.2452;  r2 = 2.47044;  x21 = 17.2743;  x22 = 4.53996;  y22 = 22.5678;  r3 = 0.8519;  x31 = 14.3729;  x32 = 4.58897;  y32 = 19.2458

小円の直径は 1.7038 寸である。

using Plots

function draw(more=false)
   pyplot(size=(1000, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 10//2
   (r0, a, x1, r2, x21, x22, y22, r3, x31, x32, y32) = res[1]
   @printf("r0 = %g;  a = %g;  x1 = %g;  r2 = %g;  x21 = %g;  x22 = %g;  y22 = %g;  r3 = %g;  x31 = %g;  x32 = %g;  y32 = %g\n", r0, a, x1, r2, x21, x22, y22, r3, x31, x32, y32)
   @printf("小円の直径 = %g\n", 2r3)
   plot()
   circle(0, 0, r0, beginangle=0, endangle=180)
   circle(x1, a + r1, r1, :blue)
   circle(-x1, a + r1, r1, :blue)
   circle(x21, a + r2, r2, :green)
   circle(-x21, a + r2, r2, :green)
   circle(x22, y22, r2, :green)
   circle(-x22, y22, r2, :green)
   circle(x31, a + r3, r3, :magenta)
   circle(-x31, a + r3, r3, :magenta)
   circle(x32, y32, r3, :magenta)
   circle(-x32, y32, r3, :magenta)
   segment(0, r0, (r0 - a)/√3, a)
   segment(0, r0, -(r0 - a)/√3, a)
   segment(-sqrt(r0^2 - a^2), a, sqrt(r0^2 - a^2), a)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, a + r1, " 大円:r1,(x1,a+r1)", :blue, :left, :vcenter)
       point(x21, a + r2, "中円:r2,(x21,a+r2) ", :green, :right, :vcenter)
       point(x22, y22, " 中円:r2,(x22,y22)", :green, :left, :vcenter)
       point(x31, a + r3, "小円:r3,(x31,a+r3)  ", :magenta, :right, :vcenter)
       point(x32, y32, " 小円:r3,(x32,y32)", :magenta, :left, :vcenter)
       point(0, a, " a", :black, :left, :top, delta=-delta)
       point(√(r0^2 - a^2), a, "(√(r0^2-a^2),a)", :black, :right, :top, delta=-delta)
       point(0, r0, " r0", :red, :left, :bottom)
       point(-(r0 - a)/√3, a, "(-(r0-a)/√3, a)", :black, :center, :top, delta=-delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その418）

算額（その418）

福島県田村市船引町北移東鳥堂 東鳥堂（東鳥神社）観音堂 文久2年(1862)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

団扇の中に弦の中点を通る2線分がある。3 個の等円は外円に接し，弦と2線分にも接している。弦と矢の和が 16 寸，弦と矢の積が 48 平方寸のとき，等円の直径を求めよ。

等円の半径と中心座標を r, (x, a + r) とし，以下の方程式の数値解を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r::positive, x::positive,
     a::positive, b::positive, 矢::positive, 弦::positive;

矢 = r0 - a
弦 = 2sqrt(r0^2 - a^2)
eq1 = 弦 + 矢 - 16
eq2 = 弦*矢 - 48
eq3 = distance(0, a, b, sqrt(r0^2 - b^2), x, a + r) - r^2
eq4 = distance(0, a, b, sqrt(r0^2 - b^2), 0, r0 - r) - r^2
eq5 = x^2 + (a + r)^2 - (r0 - r)^2;
# solve([eq1, eq2, eq3, eq4, eq5], (r, r0, x, a, b))

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 v, r.f_converged
end;

function H(u)
   (r, r0, x, a, b) = u
   return [
       -a + r0 + 2*sqrt(-a^2 + r0^2) - 16,  # eq1
       2*(-a + r0)*sqrt(-a^2 + r0^2) - 48,  # eq2
       -r^2 + (-b*(-a*r + b*x + r*sqrt(-b^2 + r0^2))/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2) + x)^2 + (a + r - (a*(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2) - (a - sqrt(-b^2 + r0^2))*(-a*r + b*x + r*sqrt(-b^2 + r0^2)))/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2))^2,  # eq3
       b^2*(a^2 + a*r - a*r0 - a*sqrt(-b^2 + r0^2) - r*sqrt(-b^2 + r0^2) + r0*sqrt(-b^2 + r0^2))^2/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2)^2 - r^2 + (-r + r0 - (-a^2*r + a^2*r0 + a*b^2 + 2*a*r*sqrt(-b^2 + r0^2) - 2*a*r0*sqrt(-b^2 + r0^2) + b^2*r - b^2*r0 - r*r0^2 + r0^3)/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2))^2,  # eq4
       x^2 + (a + r)^2 - (-r + r0)^2,  # eq5
   ]
end;

iniv = [big"1.6", 5, 3, 2.5, 2.5]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[1.499999999999999999999999999999999999999999999999999999999999999999999999999948, 6.499999999999999999999999999999999999999999999999999999999999999999999999999862, 3.000000000000000000000000000000000000000000000000000000000000000000000000000138, 2.499999999999999999999999999999999999999999999999999999999999999999999999999862, 2.594733192202055198398672253319262240463466167190260192229005823351113326367367], true)

r = 1.5;  r0 = 6.5;  x = 3;  a = 2.5;  b = 2.59473
矢 = 4;  弦 = 12
等円の直径 = 3

弦と矢の和が 16 寸，弦と矢の積が 48 平方寸のとき，弦は 12 寸，矢は 4 寸，等円の直径は 3 寸である。

算額に描かれている図は，「弦は 12 寸，矢は 4 寸」に対応していない。実際に図を描いてみると，とても「団扇」とは言えないものである。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, r0, x, a, b) = res[1]
   矢 = r0 - a
   弦 = 2sqrt(r0^2 - a^2)
   @printf("r = %g;  r0 = %g;  x = %g;  a = %g;  b = %g\n", r, r0, x, a, b)
   @printf("矢 = %g;  弦 = %g\n", 矢, 弦)
   @printf("等円の直径 = %g\n", 2r)
   plot()
   circle(0, 0, r0)
   circle(0, r0 - r, r, :blue)
   circle(x, a + r, r, :blue)
   circle(-x, a + r, r, :blue)
   segment(-sqrt(r0^2 - a^2), a, sqrt(r0^2 - a^2), a)
   segment(0, a, b, sqrt(r0^2 - b^2))
   segment(0, a, -b, sqrt(r0^2 - b^2))
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x, a + r, "r,(x,a+r)", :blue, :center, delta=-delta)
       point(0, a, " a", :black)
       point(b, √(r0^2-b^2), "(b,√(r0^2-b^2))", :red, :left, :bottom, delta=delta)
       point(√(r0^2-a^2), a, "(√(r0^2-a^2),a)", :red, :right, :top, delta=-delta)
       point(0, r0 - r, " r0-r", :blue, :left, :vcenter)
       point(0, r0, " r0", :red, :left, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;