算額（その575）

算額（その575）

福岡県久留米高良山者
藤田貞資(1789): 神壁算法
http://hyonemitsu.web.fc2.com/Korataisha.pdf

外円内に斜線を隔てて甲円 1 個，乙円 2 個，丙円 1 個がある。外円，甲円，丙円の直径が 60 寸，20 寸，28 寸のとき，乙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x, y)
丙円の半径と中心座標を r3, (0, r3 - R)
斜線と外円の交点座標を (x1, sqrt(R^2 - x1^2)), (x2, -sqrt(R^2 - x2^2))
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, x::positive, y::positive,
     r3::positive, x1::positive, x2::positive

(R, r1, r3) = (60, 20, 28) .// 2
eq1 = distance(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), 0, R - r1)  - r1^2
eq2 = distance(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), x, y) - r2^2
eq3 = distance(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), 0, r3 - R) - r3^2
eq4 = distance(-x2, -sqrt(R^2 - x2^2), x1, sqrt(R^2 - x1^2), x, y) - r2^2
eq5 = x^2 + y^2 - (R - r2)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r2, x2, y2, x1, x2))

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)
   (r2, x, y, x1, x2) = u
   return [
       (20 - (-x1^2*sqrt(900 - x2^2) - 20*x1^2 + x1*x2*sqrt(900 - x1^2) - x1*x2*sqrt(900 - x2^2) + x2^2*sqrt(900 - x1^2) - 20*x2^2 + 40*sqrt(900 - x1^2)*sqrt(900 - x2^2) + 36000)/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2 + (-x1*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) - 20*sqrt(900 - x1^2) - 20*sqrt(900 - x2^2) + 900)/2)^2/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)^2 - 100,  # eq1
       -r2^2 + (x - (-x1*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x*x1 + x*x2 + x1*x2 - y*sqrt(900 - x1^2) - y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)/2)/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))^2 + (y - (2*sqrt(900 - x1^2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) - (sqrt(900 - x1^2) + sqrt(900 - x2^2))*(x*x1 + x*x2 + x1*x2 - y*sqrt(900 - x1^2) - y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2,  # eq2
       (-16 - (-x1^2*sqrt(900 - x2^2) + 16*x1^2 + x1*x2*sqrt(900 - x1^2) - x1*x2*sqrt(900 - x2^2) + x2^2*sqrt(900 - x1^2) + 16*x2^2 - 32*sqrt(900 - x1^2)*sqrt(900 - x2^2) - 28800)/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2 + (-x1*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 16*sqrt(900 - x1^2) + 16*sqrt(900 - x2^2) + 900)/2)^2/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)^2 - 196,  # eq3
       -r2^2 + (x - (-x2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x*x1 + x*x2 + x1*x2 + y*sqrt(900 - x1^2) + y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)/2)/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))^2 + (y - (-2*sqrt(900 - x2^2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (sqrt(900 - x1^2) + sqrt(900 - x2^2))*(x*x1 + x*x2 + x1*x2 + y*sqrt(900 - x1^2) + y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2,  # eq4
       x^2 + y^2 - (30 - r2)^2,  # eq5
   ]
end;

iniv = BigFloat[29, 40, 12, 21, 28] .* (30/70)
res = nls(H, ini=iniv)

   (BigFloat[12.50000000000000000000000000000000000000000000000000000001913542531009364124636, 16.77050983124842272306880251548457176580463769708644293204781524730502299339952, 4.999999999999999999999999999999999999999999999999999999963738320809945557470043, 17.39163982499836430540468409013214849787147613031186674435362589343567963493901, 22.3606797749978969640917366873127623544061835961152572427143001079596068757832], true)

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

その他のパラメータは以下の通り。
乙円の直径 = 25;  r2 = 12.5;  x = 16.7705;  y = 5;  x1 = 17.3916;  x2 = 22.3607

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1, r3) = (60, 20, 28) .// 2
   (r2, x, y, x1, x2) = res[1]
   @printf("乙円の直径 = %g;  r2 = %g;  x = %g;  y = %g;  x1 = %g;  x2 = %g\n", 2r2, r2, x, y, x1, x2)
   plot()
   circle(0, 0, R)
   circle(0, R - r1, r1, :blue)
   circle(x, y, r2, :magenta)
   circle(-x, y, r2, :magenta)
   circle(0, r3 - R, r3, :green)
   segment(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), :gray)
   segment(x1, sqrt(R^2 - x1^2), -x2, -sqrt(R^2 - x2^2), :gray)
   if more
       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, R - r1, " 甲円:r1\n (0,R-r1)", :blue, :left, :vcenter)
       point(0, r3 - R, " 丙円:r1\n (0,r3-R)", :green, :left, :vcenter)
       point(x, y, "乙円:r2,(x,y)", :magenta, :center, :top, delta=-delta/2)
       point(x1, √(R^2 - x1^2), "(x1,√(R^2 - x1^2))", :black, :center, :bottom, delta=delta/2)
       point(x2, -√(R^2 - x2^2), "(x2,-√(R^2 - x2^2))", :black, :center, :bottom, delta=delta/2)
   end
end;

算額（その574）

算額（その574）

佐賀県神埼市 冠者神社（米光丁奉納） 平成21年
http://www.wasan.jp/saga/kanjya.html

大小の直角三角形の中に黃円，赤円が入っている。黃円の直径が 2 寸，甲斜が 4 寸，丁斜が 12 寸のとき，赤円の直径を求めよ。

乙斜，丙斜，戊斜の長さ を 「乙斜」，「丙斜」，「戊斜」
赤円の半径を r2
として，以下の eq1〜eq4 の連立方程式を解けば「乙斜，丙斜，戊斜, 赤円の半径」が求まる。この場合には座標軸の設定は不要である。
なお，図を描くために以下のように座標軸を設定し，追加パラメータとして，小さい方の直角三角形の頂点を (x, y)，黃円の中心座標を （x1，y1); y1 = 13 として eq1〜eq7 をまとめて解く。

算額の常套手段では一般解の式を求めて最後に具体的な数値を代入して解を求めるのであるが，3 変数を式に含めたままの一般式を求めるのはあまり現実的でないので，最初から具体的に数値を使用して連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 甲斜::positive, 丁斜::positive, 乙斜::positive, 丙斜::positive, 戊斜::positive, r2::positive,
     x::positive, y::positive, r1::positive, x1::positive

(r1, 甲斜, 丁斜) = (1, 4, 12)
eq1 = 乙斜 + 甲斜 - 丙斜 - 2r1
eq2 = 丙斜 + 丁斜 - 戊斜 - 2r2
eq3 = 乙斜^2 + 甲斜^2 - 丙斜^2
eq4 = 丙斜^2 + 丁斜^2 - 戊斜^2
eq5 = x^2 + (y - 丁斜)^2 - 乙斜^2
eq6 = (丙斜 - x)^2 + (y - 丁斜)^2 - 甲斜^2
eq7 = distance(0, 丁斜, x, y, x1, 丁斜 + r1) - r1^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (乙斜, 丙斜, 戊斜, r2, x, y, x1))

   2-element Vector{NTuple{7, Sym}}:
    (3, 5, 13, 2, 9/5, 48/5, 1/2)
    (3, 5, 13, 2, 9/5, 72/5, 2)

2 通りの解が得られるが，2 番目のものが適解である。
黃円の直径が 2 寸，甲斜が 4 寸，丁斜が 12 寸のとき，赤円の直径は 4 寸である。

その他のパラメータは以下の通り。
乙斜 = 3;  丙斜 = 5;  戊斜 = 13;  r2 = 2;  x = 1.8;  y = 14.4;  x1 = 2

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, 甲斜, 丁斜) = (1, 4, 12)
   (a, b, c, r2, x, y, x1) = (3, 5, 13, 2, 9/5, 72/5, 2)
   @printf("赤円の直径 = %g;  乙斜 = %g;  丙斜 = %g;  戊斜 = %g;  r2 = %g;  x = %g;  y = %g;  x1 = %g\n", 2r2, a, b, c, r2, x, y, x1)
   plot([b, x, 0, b, 0, 0], [12, y, 12, 12, 0, 12], color=:blue, lw=0.5)
   circle(r2, 12- r2, r2)
   circle(x1, 13, 1, :green)
   if more
       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(x, y, "(x,y)", :blue, :center, :bottom, delta=delta)
       point(b, 12, " (b,12)", :blue, :left, :vcenter)
       point(x1, 12 + r1, "黃円:r1,(x1,12+r1)", :green, :center, delta=-delta/2)
       point(r2, 12 - r2, "赤円:r2\n(r2,12-r2)", :red, :center, delta=-delta/2)
       plot!(xlims=(-0.5, 5.5))
       point(3, 14, "甲斜", :black, mark=false)
       point(0.1, 14, "乙斜", :black, mark=false)
       point(3.5, 11.8, "丙斜", :black, mark=false)
       point(0.2, 7, "丁斜", :black, mark=false)
       point(3, 7, "戊斜", :black, mark=false)
   end
end;

算額（その573）

算額（その573）

七十二 群馬県富岡市一ノ宮 貫前神社 嘉永2年(1849)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

群馬の算額 72-6 貫前神社 嘉永2年
http://takasakiwasan.web.fc2.com/gunnsann/g072-6.html

楕円の中に 2 個の合同な菱形と，2 個の合同な円を入れる。円は楕円の長径の端で楕円に接する最大のものとする。楕円の短径が 1 寸のとき，菱形の一辺の長さを求めよ。

楕円の長径と短径を a, b とする。
円は曲率円なので，その半径は b^2/a である。
円の半径と中心座標を r, (0, a - r); r = b^2/a
第一象限にある菱形の頂点座標を (x, y); y = b/2
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, x::positive, y::positive,
     r::positive, c::positive
r = b^2/a
y = b/2
eq1 = x^2/a^2 + y^2/b^2 - 1
eq2 = y/sqrt(x^2 + y^2) - r/(a - r);
res = solve([eq1, eq2], (a, x))

   1-element Vector{Tuple{Sym, Sym}}:
    (sqrt(5)*b, sqrt(15)*b/2)

長径は短径の √5 倍である。
短径が 1/2 寸のとき，長径は a = (1/2)√5 ≒ 1.11803，r = (1/2)^2 / ((1/2)√5) ≒ 0.223607，x = (1/2)*√15/2 ≒ 0.968246，y = 1/4 である。

a = (1/2)^2/√5
r = (1/2)^2 / ((1/2)√5)
x = (1/2)*√15/2
y = 1/4
(a, r, x, y)

   (0.11180339887498948, 0.22360679774997896, 0.9682458365518543, 0.25)

菱形の一辺の長さは sqrt(x^2 + y^2) = 1 である。

sqrt(x^2 + y^2)

   1.0

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = 1//2
   (a, x) = (sqrt(5)*b, sqrt(15)*b/2)
   y = b/2
   r = b^2/a
   c = sqrt(x^2 + y^2)
   @printf("a = %g;  r = %g;  x = %g;  y = %g;  c = %g\n", a, r, x, y, c)
   plot()
   ellipse(0, 0, a, b, color=:red)
   circle(a - r, 0, r, :green)
   circle(r - a, 0, r, :green)
   plot!([0, x, 0, -x, x, 0, -x, 0], [0, y, b, y, -y, -b, -y, 0], color=:blue, lw=0.5)
   if more
       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(x, y, " (x,y)", :green, :left, :vcenter)
       point(a - r, 0, "a-r", :green, :center, :bottom, delta=delta/2)
       point(a, 0, "a ", :green, :right, :bottom, delta=delta/2)
       point(0, b, " b", :green, :left, :bottom, delta=delta/2)
   end
end;

 

算額（その572）

算額（その572）

群馬の算額 19-2 榛名神社 文化8年
http://takasakiwasan.web.fc2.com/gunnsann/g019-2.html

正方形の中に 3 本の斜線を引き，区画された領域に大円，中円，小円を入れる。正方形の一辺の長さが 3 寸のとき，大円の直径を求めよ。

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

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive, r3::positive,
     b::positive, c::positive, d::positive, e::positive
eq1 = distance(d, 0, 0, e, r1, a - r1) - r1^2 |> expand |> simplify |> numerator
eq2 = distance(d, 0, 0, e, a - r2, r2) - r2^2 |> expand |> simplify |> numerator
eq3 = distance(d, 0, 0, e, r3, r3) - r3^2 |> expand |> simplify |> numerator
eq4 = distance(b, 0, a, c, r1, a - r1) -r1^2 |> expand |> simplify |> numerator
eq5 = distance(b, 0, a, c, a - r2, r2) -r2^2 |> expand |> simplify |> numerator
eq6 = distance(b, 0, a, c, r3, r3) -r3^2 |> expand |> simplify |> numerator
eq7 = distance(b, 0, a, c, a - r3, a - r3) -r3^2 |> expand |> simplify |> numerator;

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)
   (r1, r2, r3, b, c, d, e) = u
   return [
       a^2*d - 2*a*d*e - 2*a*d*r1 + 2*a*e*r1 + d*e^2 + 2*d*e*r1 - 2*e^2*r1 - 2*e*r1^2,  # eq1
       a^2*e - 2*a*d*e + 2*a*d*r2 - 2*a*e*r2 + d^2*e - 2*d^2*r2 + 2*d*e*r2 - 2*d*r2^2,  # eq2
       d*e - 2*d*r3 - 2*e*r3 + 2*r3^2,  # eq3
       a^4 - 2*a^3*b - 2*a^3*r1 + a^2*b^2 + 2*a^2*b*c + 4*a^2*b*r1 - 2*a^2*c*r1 - 2*a*b^2*c - 2*a*b^2*r1 + 2*a*c*r1^2 + b^2*c^2 + 2*b^2*c*r1 - 2*b*c^2*r1 - 2*b*c*r1^2,  # eq4
       a^2*c - 2*a^2*r2 - 2*a*b*c + 4*a*b*r2 - 2*a*c*r2 + 2*a*r2^2 + b^2*c - 2*b^2*r2 + 2*b*c*r2 - 2*b*r2^2,  # eq5
       2*a*b*r3 - 2*a*r3^2 + b^2*c - 2*b^2*r3 - 2*b*c*r3 + 2*b*r3^2,  # eq6
       a^4 - 2*a^3*b - 2*a^3*c - 2*a^3*r3 + a^2*b^2 + 4*a^2*b*c + 4*a^2*b*r3 + a^2*c^2 + 4*a^2*c*r3 - 2*a*b^2*c - 2*a*b^2*r3 - 2*a*b*c^2 - 6*a*b*c*r3 - 2*a*c^2*r3 - 2*a*c*r3^2 + b^2*c^2 + 2*b^2*c*r3 + 2*b*c^2*r3 + 2*b*c*r3^2,  # eq7
   ]
end;
a = 3
iniv = BigFloat[1.05, 0.7, 0.48, 0.7, 2.5, 2.0, 1.4]
res = nls(H, ini=iniv)
   (BigFloat[1.074953912118040074148850906669271748631882063090927359480694112255763934484609, 0.6824054007626747794460829207016340156591023107782284209892666737724326939646588, 0.4738479700017459009683547528030295288720148655840843288085599668979468150150135, 0.670122225679428548319738745786848395358404741740433798054273686137994090593031, 2.329877774320571451680261254213151604641595258259566201945726313849903192359496, 2.121320343559642573202533086314547117854507813065422109765019606987497365607011, 1.330325847824819753801482812391145847210905426272977123802622116991951746211069], true)

大円の直径は 2.14991 寸である。

その他のパラメータは以下の通り。

r1 = 1.07495;  r2 = 0.682405;  r3 = 0.473848;  b = 0.670122;  c = 2.32988;  d = 2.12132;  e = 1.33033

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, b, c, d, e) = [1.05, 0.7, 0.48, 0.7, 2.5, 2.0, 1.4]
   (r1, r2, r3, b, c, d, e) = res[1]
   @printf("大円の直径 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  b = %g;  c = %g;  d = %g;  e = %g\n",
       2r1, r1, r2, r3, b, c, d, e)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
   circle(r1, a - r1, r1)
   circle(a - r2, r2, r2, :green)
   circle(r3, r3, r3, :magenta)
   circle(a - r3, a - r3, r3, :magenta)
   segment(b, 0, a, c)
   segment(d, 0, 0, e)
   segment(a, a - d, a - e, a)
   if more
       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(b, 0, "b ", :black, :right, delta=3delta)
       point(d, 0, "d", :black, :left, delta=3delta)
       point(a, 0, "a ", :black, :right, delta=3delta)
       point(0, e, " e", :black, :left, delta=2delta)
       point(a, c, "(a,c) ", :black, :right, :vcenter)
       point(a, a - d, "(a,a-d) ", :black, :right, :vcenter)
       point(a - e, a, "(a-e,a)", :black, :center, :bottom, delta=delta/2)
       point(0, a, " a", :black, :left, :bottom, delta=delta/2)
       point(r1, a - r1, "大円:r1,(r1,a-r1)", :red, :center, delta=-delta/2)
       point(a - r2, r2, "中円:r2,(a-r2,r2)", :green, :center, delta=-delta/2)
       point(r3, r3, "小円:r3,(r3,r3)", :magenta, :center, delta=-delta/2)
       point(a - r3, a - r3, "小円:r3,(a-r3,a-r3)", :magenta, :center, delta=-delta/2)
   end
end;