算額（その117）

算額（その117）

一関市博物館 > 和算に挑戦 > 平成18年度出題問題(2)[中級問題]&回答例
文政13年(1830)刊『算法新書』より
https://www.city.ichinoseki.iwate.jp/museum/wasan/h18/normal.html

台形の中に図のように円が内接している。ABの長さが6cm，DCの長さが2cm,∠ABC,∠BCD=90度のとき，円の直径はいかほどか。

SymPy を使うまでもないが，円の中心からADへの距離が半径 r に等しいという方程式をとき以下のよう解が求まる。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms r::positive
eq = distance(0, 0, 4, 2r, 6 - r, r) - r^2
solve(eq, r)[1] |> println  # r は半径

   3/2

直径は 3 である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 3/2
   plot([0, 6, 6, 4, 0], [0, 0, 2r, 2r, 0])
   circle(6 - r, r, r)
   if more == true
       println("直径は $(2r)")
       point(6 - r, r, "(6-r,r)")
       point(0, 0, " A", :blue, :left, :top)
       point(6, 0, " B", :blue, :left, :top)
       point(6, 2r, " C", :blue, :left, :bottom)
       point(4, 2r, " D", :blue, :left, :bottom)
       hline!([0], color=:black, lw=0.5, ylims=(-0.5, 3.5))
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その116）

算額（その116）

一関市博物館 > 和算に挑戦 > 令和4年度出題問題(3) [上級問題] (高校生・一般向き)
文政13年(1830)刊『算法新書』より
https://www.city.ichinoseki.iwate.jp/museum/wasan/r4/hard.html

三角形を線分（界斜）で分けた2つの三角形それぞれに内接する半径の等しい2つの円がある。三角形の辺の大きい順に 15 寸，11 寸，6 寸としたとき，界斜の長さはいかほどか。

下図のように記号を定め，方程式を解く。方針としては 2 つの円それぞれの中心から接線までの距離に関する条件と，三角形の面積 20√2 について「底辺×高さ/2」とヘロンの公式による面積に関しての方程式である。

例によって，方程式が複雑になると solve() では解けなくなってしまうので nlsolve() で数値解を得ることにする。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms x1, x2, x3, x4, r
h = 8*sqrt(Sym(2))/3
eq1 = distance(0, 0, x1, h, x3, r) - r^2
eq2 = distance(x2, 0, x1, h, x3, r) - r^2
eq3 = distance(x2, 0, x1, h, x4, r) - r^2
eq4 = distance(x1, h, 15, 0, x4, r) - r^2
a = sqrt(h^2 + (x1 - x2)^2)
eq5 = r * (16 + a) - 20sqrt(Sym(2));

# solve([eq1, eq2, eq3, eq4, eq5])

println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")
println(eq5, ",")

   -r^2 + (r - 8*(16*r + 3*sqrt(2)*x1*x3)/(9*x1^2 + 128))^2 + (-3*x1*(8*sqrt(2)*r + 3*x1*x3)/(9*x1^2 + 128) + x3)^2,
   -r^2 + (r - sqrt(2)*(64*sqrt(2)*r - 24*x1*x2 + 24*x1*x3 + 24*x2^2 - 24*x2*x3)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2 + (x3 - (24*sqrt(2)*r*x1 - 24*sqrt(2)*r*x2 + 9*x1^2*x3 - 18*x1*x2*x3 + 9*x2^2*x3 + 128*x2)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2,
   -r^2 + (r - sqrt(2)*(64*sqrt(2)*r - 24*x1*x2 + 24*x1*x4 + 24*x2^2 - 24*x2*x4)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2 + (x4 - (24*sqrt(2)*r*x1 - 24*sqrt(2)*r*x2 + 9*x1^2*x4 - 18*x1*x2*x4 + 9*x2^2*x4 + 128*x2)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2,
   -r^2 + (r - sqrt(2)*(64*sqrt(2)*r + 24*x1*x4 - 360*x1 - 360*x4 + 5400)/(9*x1^2 - 270*x1 + 2153))^2 + (x4 - 3*(8*sqrt(2)*r*x1 - 120*sqrt(2)*r + 3*x1^2*x4 - 90*x1*x4 + 675*x4 + 640)/(9*x1^2 - 270*x1 + 2153))^2,
   r*(sqrt((x1 - x2)^2 + 128/9) + 16) - 20*sqrt(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=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, x2, x3, x4, r) = u
   return [
       -r^2 + (r - 8*(16*r + 3*sqrt(2)*x1*x3)/(9*x1^2 + 128))^2 + (-3*x1*(8*sqrt(2)*r + 3*x1*x3)/(9*x1^2 + 128) + x3)^2,
       -r^2 + (r - 8*(16*r - 3*sqrt(2)*x1*x2 + 3*sqrt(2)*x1*x3 + 3*sqrt(2)*x2^2 - 3*sqrt(2)*x2*x3)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2 + (x3 - (24*sqrt(2)*r*x1 - 24*sqrt(2)*r*x2 + 9*x1^2*x3 - 18*x1*x2*x3 + 9*x2^2*x3 + 128*x2)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2,
       -r^2 + (r - sqrt(2)*(64*sqrt(2)*r - 24*x1*x2 + 24*x1*x4 + 24*x2^2 - 24*x2*x4)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2 + (x4 - (24*sqrt(2)*r*x1 - 24*sqrt(2)*r*x2 + 9*x1^2*x4 - 18*x1*x2*x4 + 9*x2^2*x4 + 128*x2)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2,
       -r^2 + (r - sqrt(2)*(64*sqrt(2)*r + 24*x1*x4 - 360*x1 - 360*x4 + 5400)/(9*x1^2 - 270*x1 + 2153))^2 + (x4 - 3*(8*sqrt(2)*r*x1 - 120*sqrt(2)*r + 3*x1^2*x4 - 90*x1*x4 + 675*x4 + 640)/(9*x1^2 - 270*x1 + 2153))^2,
       r*(sqrt((x1 - x2)^2 + 128/9) + 16) - 20*sqrt(2),
   ]
end;

iniv = [10.8, 9.5, 8.3, 11.7, 1.25]
res = nls(H, ini=iniv)
println(res)

   ([10.333333333333236, 8.99999999999995, 7.999999999999937, 10.999999999999943, 1.4142135623730963], false)

多分，x1 = 10.333333333333236 = 31/3，r = 1.4142135623651715 = √2 であろう。
界斜の長さは

   h = 8*sqrt(2)/3
   a = sqrt(h^2 + (x1 - x2)^2) = 4 である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
 plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, x2, x3, x4, r) = res[1]
   h = 8*sqrt(Sym(2))/3
   a = sqrt(h^2 + (x1 - x2)^2)
   plot([0, 15, x1, 0], [0, 0, h, 0])
   circle(x3, r, r)
   circle(x4, r, r)
   segment(x2, 0, x1, h)
   if more == true
       @printf("x1 = %.3f, h = %.3f, x2 = %.3f, x3 = %.3f, x4 = %.3f, r = %.5f\n", x1, h, x2, x3, x4, r)
       @printf("界斜の長さ = %.1f", a)
       point(x1, h, " (x1,h)", :blue, :left, :bottom)
       point(x2, 0, "x2", :blue)
       point(x3, r, "(x3,r)", :red, :top) 
       point(x4, r, "(x4,r)", :red, :top) 
       hline!([0], color=:black, lw=0.5, ylims=(-1, 5))
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end
draw(true)
savefig("fig1.png");

   x1 = 10.333, h = 3.771, x2 = 9.000, x3 = 8.000, x4 = 11.000, r = 1.41421
   界斜の長さ = 4.0

算額（その115）

算額（その115）

和算に挑戦 －平成 24 年度中級解答例－
埼玉県桶川市 氷川天満神社に明治 43 年(1910)に奉納された算額の問 題をもとにしました
https://www.city.ichinoseki.iwate.jp/museum/wasan/h24/images/normal.pdf

新潟県三島郡出雲崎町 滝谷薬師堂 明治2年(1869)4月
http://www.wasan.jp/niigata/izumozakiyakusi.html

一辺の長さが 2尺5寸 の正三角形の中に大円と小円 3 個が内接している。それぞれの径を求めよ。

正三角形の一辺の長さを 25 寸として，以下のように記号を定め方程式を解く。SymPy なしでも解けるが，横着して円の中心から右斜辺への距離に関しての方程式を作る。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms r1, r2, x;
eq1 = distance(25//2, 0, 0, 25sqrt(Sym(3))/2, 0, r1) - r1^2;
eq2 = distance(25//2, 0, 0, 25sqrt(Sym(3))/2, 0, 2r1 + r2) - r2^2;
eq3 = distance(25//2, 0, 0, 25sqrt(Sym(3))/2, x, r2) - r2^2;

res = solve([eq1, eq2, eq3], (r1, r2, x));

for i = 1:8
   (r1, r2, x) = res[i]
   if r1 > 0 && r2 > 0 && x < 25/2
       println("i = $i")
       println("r1 = $r1 = $(r1.evalf())")
       println("r2 = $r2 = $(r2.evalf())")
       println("x  = $x = $(x.evalf())")
   end
end

   i = 7
   r1 = 25*sqrt(3)/6 = 7.21687836487032
   r2 = 25*sqrt(3)/18 = 2.40562612162344
   x  = 25/3 = 8.33333333333333

大円，小円の直径は元の単位では 7.21687836487032 * 2/10 尺，2.40562612162344 * 2/10 尺である。

println("$(7.21687836487032 * 2/10) 尺, $(2.40562612162344 * 2/10) 尺")

   1.443375672974064 尺, 0.481125224324688 尺

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x) = res[7]
   plot([25/2, 0, -25/2, 25/2], [0, 25sqrt(3)/2, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1)
   circle(0, 2r1+r2, r2, :blue)
   circle(x, r2, r2, :blue)
   circle(-x, r2, r2, :blue)
   if more == true
       @printf("r1 = %.5f,  r2 = %.5f,  x = %.5f\n", r1, r2, x)
       @printf("大円の径 = %.5f 尺,  小円の径 = %.5f 尺\n", 2r1/10, 2r2/10)
       point(0, r1, " r1", :red)
       point(0, 2r1+r2, "2r1+r2", :blue)
       point(x, r2, "(x,r2)", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

   r1 = 7.21688,  r2 = 2.40563,  x = 8.33333
   大円の径 = 1.44338 尺,  小円の径 = 0.48113 尺

 

算額（その114）

算額（その114）

福島県田村郡三春町平沢担橋 諏訪神社 大正15年(1926)8月
http://www.wasan.jp/fukusima/miharusuwa.html

三角形の中に菱形，甲円，乙円が入っている。大斜（一番長い辺）が 9 寸，甲円の径が 2 寸，乙円の径が 1 寸のとき，菱形の辺の長さはいかほどか。

以下に，求めるものは異なるが同じ問題がある。

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（222)
長野県東筑摩郡筑北村青柳 碩水寺豊川稲荷社 慶応4年(1868)

方程式を書きやすくするために左右反転させて，大斜を 18，甲円，乙円の半径を 2, 1，菱形の辺の長さを a として下図のように記号を定める。

算額の問のように，「菱形の辺の長さ」だけを求めるならば，三角形の相似関係を使えば簡単に答えを得ることができる。⊿EFC，⊿ODE，⊿OBC は相似（相似比 1:2:3）

ここでは，解を確認するための図を描くのに必要な座標をすべて求めることにする。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms a, x, y, x1, x2, y2, x3, y3;
eq0 = (18 - a) + (18 - a)/2 - 18
eq1 = distance(0, 0, x, y, x1, 2) - 2^2;
eq2 = distance(0, 0, x, y, x2, y2) - 1^2;
eq3 = distance(x3, y3, 18 - a, 0, x1, 2) - 2^2;
eq4 = distance(x, y, 18, 0, x2, y2) - 1^2;
eq5 = y3^2 + (18 - a - x3)^2 - 6^2;
eq6 = 2x - 3x3;
eq7 = 2y - 3y3;

例により，solve() では解けないので，nlsolve() で数値解を求めることにする。

# res = solve([eq0, eq1, eq2, eq3, eq4, eq5, eq6, eq7]);

eq0 |> println
eq1 |> println
eq2 |> println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println

   9 - 3*a/2
   (-x*(x*x1 + 2*y)/(x^2 + y^2) + x1)^2 + (-y*(x*x1 + 2*y)/(x^2 + y^2) + 2)^2 - 4
   (-x*(x*x2 + y*y2)/(x^2 + y^2) + x2)^2 + (-y*(x*x2 + y*y2)/(x^2 + y^2) + y2)^2 - 1
   (x1 - (x1*x3^2 - 24*x1*x3 + 144*x1 + 2*x3*y3 + 12*y3^2 - 24*y3)/(x3^2 - 24*x3 + y3^2 + 144))^2 + (-y3*(x1*x3 - 12*x1 - 12*x3 + 2*y3 + 144)/(x3^2 - 24*x3 + y3^2 + 144) + 2)^2 - 4
   (x2 - (x2*(x^2 - 36*x + y^2 + 324) + y*(x*y2 - x2*y + 18*y - 18*y2))/(x^2 - 36*x + y^2 + 324))^2 + (-y*(x*x2 - 18*x - 18*x2 + y*y2 + 324)/(x^2 - 36*x + y^2 + 324) + y2)^2 - 1
   y3^2 + (12 - x3)^2 - 36
   2*x - 3*x3
   2*y - 3*y3

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)
   (a, x, y, x1, x2, y2, x3, y3) = u
   return [
       9 - 3*a/2,
       (-x*(x*x1 + 2*y)/(x^2 + y^2) + x1)^2 + (-y*(x*x1 + 2*y)/(x^2 + y^2) + 2)^2 - 4,
       (-x*(x*x2 + y*y2)/(x^2 + y^2) + x2)^2 + (-y*(x*x2 + y*y2)/(x^2 + y^2) + y2)^2 - 1,
       (x1 - (x1*x3^2 - 24*x1*x3 + 144*x1 + 2*x3*y3 + 12*y3^2 - 24*y3)/(x3^2 - 24*x3 + y3^2 + 144))^2 + (-y3*(x1*x3 - 12*x1 - 12*x3 + 2*y3 + 144)/(x3^2 - 24*x3 + y3^2 + 144) + 2)^2 - 4,
       (x2 - (x^2*x2 - 36*x*x2 + x*y*y2 + 324*x2 + 18*y^2 - 18*y*y2)/(x^2 - 36*x + y^2 + 324))^2 + (-y*(x*x2 - 18*x - 18*x2 + y*y2 + 324)/(x^2 - 36*x + y^2 + 324) + y2)^2 - 1,
       y3^2 + (12 - x3)^2 - 36,
       2*x - 3*x3,
       2*y - 3*y3,
   ]
end;

iniv = [5.0, 12.0, 7, 8, 12, 5, 8, 4]
res = nls(H, ini=iniv)
println(res)

   ([6.0, 12.125711439226585, 6.8185580517266695, 7.63711610345335, 11.902365677877736, 5.5457053678177815, 8.083807626151057, 4.54570536781778], true)

a = 6 となるが，元の単位では 3 寸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, x, y, x1, x2, y2, x3, y3) = res[1]
   @printf("a  = %8.5f,  x = %8.5f,  y = %8.5f, x1 = %8.5f\n", a, x, y, x1)
   @printf("x2 = %8.5f, y2 = %8.5f, x3 = %8.5f, y3 = %8.5f\n", x2, y2, x3, y3)
   plot([0, 18, x, 0], [0, 0, y, 0], color=:black, lw=0.5)
   circle(x1, 2, 2)
   circle(x2, y2, 1)
   plot!([x3+a, x3, 18-a], [y3, y3, 0], color=:black, lw=0.5)
   if more == true
       point(x1, 2, "(x1,2),2", :green, :top)
       point(x3, y3, "E:(x3,y3) ", :green, :right, :bottom)
       point(x3 + a, y3, "  F:(x3+a, y3)")
       point(x, y, "  C:(x,y),1", :green, :center, :bottom)
       point(x2, y2, "(x2,y2),1", :green, :top)
       point(18 - a, 0, " D:18-a")
       point(18, 0, "  B:18")
       point(0, 0, "O ", :green, :right)
       hline!([0], color=:black, lw=0.5, xlims=(-1, 20))
       vline!([0], color=:black, lw=0.5, ylims=(-1, 7))
   else
       plot!(showaxis=false)
   end
end;

   a  =  6.00000,  x = 12.12571,  y =  6.81856, x1 =  7.63712
   x2 = 11.90237, y2 =  5.54571, x3 =  8.08381, y3 =  4.54571

算額（その113）

算額（その113）

(16) 京都府京都市右京区山ノ内宮脇町 山王神社 明治23年(1890)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.

愛媛県松山市吉藤町 三島神社 明治13年(1880)7月
http://www.wasan.jp/ehime/misima.html

兵庫県上郡町 上郡天満神社 文化5年(1808)
https://kamigori-town-museum.jp/cultural-asetts/559/

鉤が 108間6合，股が144間8合の直角三角形の内部に円と正三角形が入っている。正三角形の一辺の長さはいかほどか?

鉤股を 1086, 1448 とする。
円の半径を r，正三角形の頂点座標を，x1, (x2, y2), x3 として以下のように記号を定め方程式を解く。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms r, x1, x2, y2;
eq1 = distance(0, 1086, 1448, 0, r, r) - r^2;
eq2 = (x2 - x1)^2 + y2^2  - 4(x2 - x1)^2 |> expand
eq3 = distance(x1, 0, x2, y2, r, r) - r^2
eq4 = 1448y2 - 1086*(1448 - x2)

res = solve([eq1, eq2, eq3, eq4], (r, x1, x2, y2));
res

   8-element Vector{NTuple{4, Sym}}:
    (362, 362 - 362*sqrt(3), 5792/13 - 1448*sqrt(3)/13, 1086*sqrt(3)/13 + 9774/13)
    (362, 362 - 362*sqrt(3)/3, -1448*sqrt(3)/3, 362*sqrt(3) + 1086)
    (362, 362 + 362*sqrt(3), 1448*sqrt(3)/13 + 5792/13, 9774/13 - 1086*sqrt(3)/13)
    (362, 362*sqrt(3)/3 + 362, 1448*sqrt(3)/3, 1086 - 362*sqrt(3))
    (2172, 2172 - 2172*sqrt(3), 4344 - 2896*sqrt(3), -2172 + 2172*sqrt(3))
    (2172, 2172 - 724*sqrt(3), 21720/13 - 8688*sqrt(3)/13, -2172/13 + 6516*sqrt(3)/13)
    (2172, 2172 + 2172*sqrt(3), 4344 + 2896*sqrt(3), -2172*sqrt(3) - 2172)
    (2172, 724*sqrt(3) + 2172, 8688*sqrt(3)/13 + 21720/13, -6516*sqrt(3)/13 - 2172/13)

for i = 1:8
   println("r  = $(res[i][1].evalf())")
   println("x1 = $(res[i][2].evalf())")
   println("x2 = $(res[i][3].evalf())")
   println("y2 = $(res[i][4].evalf())")
   println()
end

   r  = 362.000000000000
   x1 = -265.002392339934
   x2 = 252.614648510790
   y2 = 896.539013616908

   r  = 362.000000000000
   x1 = 152.999202553355
   x2 = -836.003189786578
   y2 = 1713.00239233993

   r  = 362.000000000000
   x1 = 989.002392339934
   x2 = 638.462274566133
   y2 = 607.153294075400

   r  = 362.000000000000
   x1 = 571.000797446645
   x2 = 836.003189786578
   y2 = 458.997607660066

   r  = 2172.00000000000
   x1 = -1590.01435403960
   x2 = -672.019138719469
   y2 = 1590.01435403960

   r  = 2172.00000000000
   x1 = 917.995215320133
   x2 = 513.226352603200
   y2 = 701.080235547600

   r  = 2172.00000000000
   x1 = 5934.01435403960
   x2 = 9360.01913871947
   y2 = -5934.01435403960

   r  = 2172.00000000000
   x1 = 3426.00478467987
   x2 = 2828.31210893526
   y2 = -1035.23408170145

4 番目の解が妥当である。

正三角形の一辺の長さは x1 + 2(x2-x1) = 530.004784679867，もとの単位では 53間強である。

144間8合 × (√3 - 1)/2 = 53.00047846798671 間

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
  plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x1, x2, y2) = res[4]
   x3 = x1 + 2(x2 - x1)
   plot([0, 1448, 0, 0], [0, 0, 1086, 0], color=:black, lw=0.5)
   circle(r, r, r)
   segment(x1, 0, x2, y2, :green)
   segment(x2, y2, x3, 0, :green)
   if more == true
       println("x = $(sqrt((x2-x1)^2 + y2^2)) = $((sqrt((x2-x1)^2 + y2^2)).evalf())")
       println("x3 - x1 = $(x3 - x1) = $((x3 - x1).evalf())")
       point(r, r, "(r,r)", :green, :center)
       point(x1, 0, " x1", :green, :left, :bottom)
       point(x3, 0, " x3", :green, :left, :bottom)
       point(x2, y2, " (x2,y2)", :green, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   x = sqrt((-362 + 362*sqrt(3))^2 + (1086 - 362*sqrt(3))^2) = 530.004784679867
   x3 - x1 = -724 + 724*sqrt(3) = 530.004784679867

算額（その112）

算額（その112）

愛媛県松山市 伊佐爾波神社 嘉永3年(1850)5月
http://www.wasan.jp/ehime/isaniwa14.html

円内に 2 つの甲円，3 つの乙円が入っている。それぞれの円の半径を求めよ。

算額（その111）に似ているのだが，写真が不鮮明であるが，下の 2 つの乙円と甲円がそれぞれ直線に接している。

外円の半径を 1 とする。甲円，乙円の中心座標と半径を (r1, y1, r1)，(r2, y2, r2) とおき，以下の方程式を立てて解く。
算額（その111）と違うのは eq4 である。しょうもない条件であるが，これでは以下のように無意味な解しか得られない。

using SymPy
@syms r1, r2, y1, y2;
eq3 = r1^2 + (1 - r2 - y1)^2 - (r1 + r2)^2 |> expand
eq4 = y1 - (y2 + r2) - r1
eq5 = r1^2 + y1^2 - (1 - r1)^2 |> expand
eq6 = r2^2 + y2^2 - (1 - r2)^2 |> expand;

solve([eq3, eq4, eq5, eq6], (r1, r2, y1, y2))

   1-element Vector{NTuple{4, Sym}}:
    (0, 0, 1, 1)

そこで， nlsolve() で数値解を求めることにする。

eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println

   -2*r1*r2 + 2*r2*y1 - 2*r2 + y1^2 - 2*y1 + 1
   -r1 - r2 + y1 - y2
   2*r1 + y1^2 - 1
   2*r2 + y2^2 - 1

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)
  (r1, r2, y1, y2) = u
  return [
       -2*r1*r2 + 2*r2*y1 - 2*r2 + y1^2 - 2*y1 + 1,
       -r1 - r2 + y1 - y2,
       2*r1 + y1^2 - 1,
       2*r2 + y2^2 - 1
   ]
end;

iniv = [0.5, 0.3, 0.1, -0.6]
res = nls(H, ini=iniv)
println(res)

ちゃんと解が求まった。

   ([0.49309632821478017, 0.28307747532588534, 0.11750465339908704, -0.6586691501415785], true)

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, y1, y2) = [0.49309632821478017, 0.28307747532588534, 0.11750465339908704, -0.6586691501415785]
   plot()
   circle(0, 0, 1)
   circle(r2, y2, r2, :blue)
   circle(-r2, y2, r2, :blue)
   circle(0, 1 - r2, r2, :blue)
   circle(r1, y1, r1, :green)
   circle(-r1, y1, r1, :green)
   hline!([y2+r2], color=:black, lw=0.5)
   if more == true
       @printf("r1 = %.5f; r2 = %.5f;  y1 = %.5f;  y2 = %.5f", r1, r2, y1, y2)
       point(r1, y1, "(r1,y1)", :green, :center)
       point(0, 1 - r2, "1-r2 ", :blue, :right)
       point(r2, y2, "(r2,y2)", :blue, :center)
       point(0, y2+r2, "y2+r2", :blue, :center)
       point(r1, y2+r2, "(r1, y2+r2)")
       hline!([0, y2+r2], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.49310; r2 = 0.28308;  y1 = 0.11750;  y2 = -0.65867

 

算額（その111）

算額（その111）

愛媛県松山市 伊佐爾波神社 嘉永3年(1850)正月
http://www.wasan.jp/ehime/isaniwa7.html

円内に 2 つの甲円，3 つの乙円が入っている。それぞれの円の半径を求めよ。

外円の半径を 1 とする。甲円，乙円の中心座標と半径を (r1, y1, r1)，(r2, y2, r2) とおき，以下の方程式を立てて解く。

using SymPy
@syms r1, r2, y1, y2;
eq3 = r1^2 + (1 - r2 - y1)^2 - (r1 + r2)^2
eq4 = (r1 - r2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq5 = r1^2 + y1^2 - (1 - r1)^2
eq6 = r2^2 + y2^2 - (1 - r2)^2;

res = solve([eq3, eq4, eq5, eq6], (r1, r2, y1, y2))

   5-element Vector{NTuple{4, Sym}}:
    (0, 0, 1, 1)
    ((-142*sqrt(2) - 197 + 65*sqrt(-13 + 16*sqrt(2)) + 46*sqrt(-26 + 32*sqrt(2)))/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -16*sqrt(2)/49 + 2*sqrt(-13 + 16*sqrt(2))/49 + 13/49 + 4*sqrt(-26 + 32*sqrt(2))/49, (-10*sqrt(-26 + 32*sqrt(2)) - 13*sqrt(-13 + 16*sqrt(2)) + 30*sqrt(2) + 43)/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -4*sqrt(2)/7 - 2/7 + sqrt(-13 + 16*sqrt(2))/7)
    (-(197 + 142*sqrt(2) + 65*sqrt(-13 + 16*sqrt(2)) + 46*sqrt(-26 + 32*sqrt(2)))/(sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -16*sqrt(2)/49 - 4*sqrt(-26 + 32*sqrt(2))/49 - 2*sqrt(-13 + 16*sqrt(2))/49 + 13/49, (13*sqrt(-13 + 16*sqrt(2)) + 30*sqrt(2) + 43 + 10*sqrt(-26 + 32*sqrt(2)))/(sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -4*sqrt(2)/7 - sqrt(-13 + 16*sqrt(2))/7 - 2/7)
    ((-65*sqrt(13 + 16*sqrt(2)) + 46*sqrt(26 + 32*sqrt(2)) - 142*sqrt(2)*I + 197*I)/(sqrt(13 + 16*sqrt(2)) - 9*I + 4*sqrt(2)*I), 13/49 + 16*sqrt(2)/49 - 2*sqrt(-16*sqrt(2) - 13)/49 + 4*sqrt(-32*sqrt(2) - 26)/49, (-10*sqrt(26 + 32*sqrt(2)) + 13*sqrt(13 + 16*sqrt(2)) - 43*I + 30*sqrt(2)*I)/(sqrt(13 + 16*sqrt(2)) - 9*I + 4*sqrt(2)*I), -2/7 + 4*sqrt(2)/7 - I*sqrt(13 + 16*sqrt(2))/7)
    ((-65*sqrt(13 + 16*sqrt(2)) + 46*sqrt(26 + 32*sqrt(2)) - 197*I + 142*sqrt(2)*I)/(sqrt(13 + 16*sqrt(2)) - 4*sqrt(2)*I + 9*I), 13/49 + 16*sqrt(2)/49 - 4*sqrt(-32*sqrt(2) - 26)/49 + 2*sqrt(-16*sqrt(2) - 13)/49, (-10*sqrt(26 + 32*sqrt(2)) + 13*sqrt(13 + 16*sqrt(2)) - 30*sqrt(2)*I + 43*I)/(sqrt(13 + 16*sqrt(2)) - 4*sqrt(2)*I + 9*I), -2/7 + 4*sqrt(2)/7 + I*sqrt(13 + 16*sqrt(2))/7)

5 組の解が得られるが，2番目のものだけが適切である。なお，@syms において各変数を ::positive と宣言すると解が求まらないので注意。

for i = 1:5
   println("r1 = $(res[i][1].evalf())")
   println("r2 = $(res[i][2].evalf())")
   println("y1 = $(res[i][3].evalf())")
   println("y2 = $(res[i][4].evalf())")
   println()
end

   r1 = 0
   r2 = 0
   y1 = 1.00000000000000
   y2 = 1.00000000000000

   r1 = 0.494520180083009
   r2 = 0.288374102478171
   y1 = 0.104688298457764
   y2 = -0.650578046850383

   r1 = -45.1219371780525
   r2 = -0.681329898313661
   y1 = 9.55216595103462
   y2 = -1.53709459586173

   r1 = 0.31370849898476 - 0.463999436044365*I
   r2 = 0.727090142815704 + 0.445454898991966*I
   y1 = -0.82842712474619 - 0.560096865715887*I
   y2 = 0.522407749927483 - 0.852695809075959*I

   r1 = 0.31370849898476 + 0.463999436044365*I
   r2 = 0.727090142815704 - 0.445454898991966*I
   y1 = -0.82842712474619 + 0.560096865715887*I
   y2 = 0.522407749927483 + 0.852695809075959*I

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")

   (r1, r2, y1, y2) = (
       (-142*sqrt(2) - 197 + 65*sqrt(-13 + 16*sqrt(2)) + 46*sqrt(-26 + 32*sqrt(2)))/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9),
       -16*sqrt(2)/49 + 2*sqrt(-13 + 16*sqrt(2))/49 + 13/49 + 4*sqrt(-26 + 32*sqrt(2))/49,
       (-10*sqrt(-26 + 32*sqrt(2)) - 13*sqrt(-13 + 16*sqrt(2)) + 30*sqrt(2) + 43)/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9),
       -4*sqrt(2)/7 - 2/7 + sqrt(-13 + 16*sqrt(2))/7)

   plot()
   circle(0, 0, 1)
   circle(r2, y2, r2, :blue)
   circle(-r2, y2, r2, :blue)
   circle(0, 1 - r2, r2, :blue)
   circle(r1, y1, r1, :green)
   circle(-r1, y1, r1, :green)
   if more == true
       @printf("r1 = %.5f; r2 = %.5f;  y1 = %.5f;  y2 = %.5f", r1, r2, y1, y2)
       point(r1, y1, "(r1,y1)", :green, :center)
       point(0, 1 - r2, "1-r2 ", :blue, :right)
       point(r2, y2, "(r2,y2)", :blue, :center)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.49452; r2 = 0.28837;  y1 = 0.10469;  y2 = -0.65058

算額（その110）

算額（その110）

愛媛県松山市 伊佐爾波神社 嘉永3年(1850)正月
http://www.wasan.jp/ehime/isaniwa7.html

正三角形の底辺に内接し，隣同士の円も接する3つの円が，その上にあるもう1つの円とも接するような場合，大円，中円，小円の半径を求めよ。

算額（その4）に似ている。

正三角形の一辺の長さを 2 とし，大円，中円，小円の半径を r1, r2, r3 とし，中円の中心座標を (x2, r2) とする。x2 = 1 - √3*r2 である。

using SymPy
@syms r1::positive, r2::positive, r3::positive;

x2 = 1 - √Sym(3)r2
eq1 = (r2 - r3)^2 + x2^2 - (r2 + r3)^2;
eq2 = sqrt(Sym(3)) - r1 - 2r3 - 2r1;
eq3 = (r1 + 2r3 - r2)^2 + x2^2 - (r1 + r2)^2;

res = solve([eq1, eq2, eq3], (r1, r2, r3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-20/33 - sqrt(-113 + 72*sqrt(3))/33 + 6*sqrt(3)/11, sqrt(-113 + 72*sqrt(3))/11 + 2*sqrt(-339 + 216*sqrt(3))/33 + 8*sqrt(3)/33 + 6/11, -sqrt(3)*(-100*sqrt(3)/33 - 10*sqrt(-113/1452 + 6*sqrt(3)/121) + 35/11)/10)
    (-20/33 + sqrt(-113 + 72*sqrt(3))/33 + 6*sqrt(3)/11, -2*sqrt(-339 + 216*sqrt(3))/33 - sqrt(-113 + 72*sqrt(3))/11 + 8*sqrt(3)/33 + 6/11, -sqrt(3)*(-100*sqrt(3)/33 + 10*sqrt(-113/1452 + 6*sqrt(3)/121) + 35/11)/10)

(-20/33 - sqrt(-113 + 72*sqrt(3))/33 + 6*sqrt(3)/11, sqrt(-113 + 72*sqrt(3))/11 + 2*sqrt(-339 + 216*sqrt(3))/33 + 8*sqrt(3)/33 + 6/11, -sqrt(3)*(-100*sqrt(3)/33 - 10*sqrt(-113/1452 + 6*sqrt(3)/121) + 35/11)/10)

   (0.23500815165114197, 1.6355839657205662, 0.5135131763077257)

(-20/33 + sqrt(-113 + 72*sqrt(3))/33 + 6*sqrt(3)/11, -2*sqrt(-339 + 216*sqrt(3))/33 - sqrt(-113 + 72*sqrt(3))/11 + 8*sqrt(3)/33 + 6/11, -sqrt(3)*(-100*sqrt(3)/33 + 10*sqrt(-113/1452 + 6*sqrt(3)/121) + 35/11)/10)

   (0.4423806081209667, 0.2951073349188893, 0.20245449160298862)

r1 > r2 > r3 ゆえ，2 番目の解が適切である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = res[2]
   (r1, r2, r3) = (r1.evalf(), r2.evalf(), r3.evalf())
   x2 = 1-√3 * r2
   plot([-1, 1, 0, -1], [0, 0, √3, 0])
   circle(0, r1 + 2r3, r1)
   circle(0, r3, r3, :blue)
   circle(x2, r2, r2, :green)
   circle(-x2, r2, r2, :green)
   if more == true
       @printf("r1 = %.5f,  r2 = %.5f,  r3 = %.5f, x2 = %.5f", r1, r2, r3, x2)
       point(0, r1 + 2r3, "r1+2r3 ", :red, :right)
       point(0, r3, "r3 ", :blue, :right)
       point(x2, r2, "(x2,r2)", :green, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.44238,  r2 = 0.29511,  r3 = 0.20245, x2 = 0.48886

算額（その109）

算額（その109）

東京都府中市宮町 大國魂神社 明治18年
大國魂神社の算額（下段第5問）についての考察
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/ok1.pdf

外円内に三角形，甲円，乙円，丙円が入っている。甲円，乙円，丙円の径をそれぞれ 49 寸，28 寸，17 寸のとき，外円の径を求めよ。

図のように記号を定め，方程式を立てる。
三角形の頂点 A, B, C の座標をそれぞれ (-x0, 98 - R), (x0, 98 - R), (x4, y4) とする。
甲円，乙円，丙円の中心座標と半径をそれぞれ (0, 49 -R; R), (x2, 41, 28), (x3, y3, 17) とする。

方程式は solve() では解けないので，nlsolve() で数値解を求める。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms R, x2, x3, y3, x0, x4, y4;

eq1 = x3^2 + y3^2 - (R - 17)^2 |> expand;
eq2 = R^2 - (98 - R)^2 - x0^2 |> expand;
eq3 = distance(x0, 98 - R, -x0, 98 - R, x2, 41) - 28^2 |> expand;
eq4 = x4^2 + y4^2 - R^2 |> expand;
eq5 = distance(-x0, 98 - R, x4, y4, x2, 41) - 28^2;
eq6 = distance( x0, 98 - R, x4, y4, x2, 41) - 28^2;
eq7 = distance(-x0, 98 - R, x4, y4, x3, y3) - 17^2;

eq1 |> println
eq2 |> println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println

   -R^2 + 34*R + x3^2 + y3^2 - 289
   196*R - x0^2 - 9604
   R^2 - 114*R + 2465
   -R^2 + x4^2 + y4^2
   (41 - (41*R^2 + 82*R*y4 - 8036*R + 41*x0^2 + 82*x0*x4 + 41*x4^2 + 41*y4^2 - 8036*y4 + (x0 + x4)*(R*x2 - R*x4 + x0*y4 - 41*x0 + x2*y4 - 98*x2 + 57*x4) + 393764)/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (x2 - (x2*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) - (R + y4 - 98)*(R*x2 - R*x4 + x0*y4 - 41*x0 + x2*y4 - 98*x2 + 57*x4))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 - 784
   (x2 - (R^2*x4 + R*x0*y4 - 41*R*x0 + R*x4*y4 - 155*R*x4 + x0^2*x2 - 2*x0*x2*x4 + x0*y4^2 - 139*x0*y4 + 4018*x0 + x2*x4^2 - 57*x4*y4 + 5586*x4)/(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (-((98 - R)*(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) + (R + y4 - 98)*(R^2 + R*y4 - 155*R + x0^2 - x0*x2 - x0*x4 + x2*x4 - 57*y4 + 5586))/(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) + 41)^2 - 784
   (x3 - (x3*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) - (R + y4 - 98)*(R*x3 - R*x4 - x0*y3 + x0*y4 + x3*y4 - 98*x3 - x4*y3 + 98*x4))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (y3 - (y3*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) + (x0 + x4)*(R*x3 - R*x4 - x0*y3 + x0*y4 + x3*y4 - 98*x3 - x4*y3 + 98*x4))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 - 289

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10])

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)
   (R, x2, x3, y3, x0, x4, y4) = u
   return [
       -R^2 + 34*R + x3^2 + y3^2 - 289,
       196*R - x0^2 - 9604,
       R^2 - 114*R + 2465,
       -R^2 + x4^2 + y4^2,
       (41 - (41*R^2 + 82*R*y4 - 8036*R + 41*x0^2 + 82*x0*x4 + 41*x4^2 + 41*y4^2 - 8036*y4 + (x0 + x4)*(R*x2 - R*x4 + x0*y4 - 41*x0 + x2*y4 - 98*x2 + 57*x4) + 393764)/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (x2 - (x2*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) - (R + y4 - 98)*(R*x2 - R*x4 + x0*y4 - 41*x0 + x2*y4 - 98*x2 + 57*x4))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 - 784,
       (41 - (41*R^2 + 82*R*y4 - 8036*R + 41*x0^2 - 82*x0*x4 + 41*x4^2 + 41*y4^2 - 8036*y4 - (x0 - x4)*(R*x2 - R*x4 - x0*y4 + 41*x0 + x2*y4 - 98*x2 + 57*x4) + 393764)/(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (x2 - (x2*(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) - (R + y4 - 98)*(R*x2 - R*x4 - x0*y4 + 41*x0 + x2*y4 - 98*x2 + 57*x4))/(R^2 + 2*R*y4 - 196*R + x0^2 - 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 - 784,
       (x3 - (-x0*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) + (x0 + x4)*(R^2 + R*y3 + R*y4 - 196*R + x0^2 + x0*x3 + x0*x4 + x3*x4 + y3*y4 - 98*y3 - 98*y4 + 9604))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 + (y3 - ((98 - R)*(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604) + (R + y4 - 98)*(R^2 + R*y3 + R*y4 - 196*R + x0^2 + x0*x3 + x0*x4 + x3*x4 + y3*y4 - 98*y3 - 98*y4 + 9604))/(R^2 + 2*R*y4 - 196*R + x0^2 + 2*x0*x4 + x4^2 + y4^2 - 196*y4 + 9604))^2 - 289,
   ]
end;

iniv = [big"80.0", 30, -30, 55, 80, 30, 75]
res = nls(H, ini=iniv)
println(res)

   (BigFloat[85.0, 28.00000000000000000000000000000000000000000000000000000000000000000000000000055, -31.99999993545069734104462106203685900316672771867706686845499772119474317529376, 60.00000003442629475144286876691367519831107855003889767015735701621077260821287, 84.0, 35.99999999999999999999999999999999999999999999999999999999999999999999999999945, 77.0], true)

外円の半径は 85.0 である。元の単位での径（直径）は 85 寸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, x2, x3, y3, x0, x4, y4) = res[1]
   println("R = $R")
   plot()
   circle(0, 0, R)
   circle(0, 49-R, 49, :blue)
   circle(x2, 41, 28, :magenta)
   circle(x3, y3, 17, :green)
   plot!([x0, x4, -x0, x0], [98 - R, y4, 98 - R, 98 - R], color=:brown, lw=0.5)
    if more
       @printf("x0 = %.3f, x2 = %.3f, x3 = %.3f, y3 = %.3f, x4 = %.3f, y4 = %.3f", x0, x2, x3, y3, x4, y4)
       point(0, 49 - R, " 甲(0,49-R;49)", :blue, :left)
       point(x2, 41, "乙(x2,41;28)", :magenta, :top)
       point(x3, y3, "丙(x3,y3;17)", :green, :top)
       point(0, 0, " O")
       point(-x0, 98 - R, " A(-x0,98-R)", :brown)
       point(x0, 98 - R, "B(x0,98-R) ", :brown, :right)
       point(x4, y4, " C(x4,y4)", :brown, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
    end
end;

 

 

算額（その108）

算額（その108）

東京都府中市宮町 大國魂神社 明治18年
大國魂神社の算額（上段第2問）についての考察
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/ok2.pdf

外円内に甲円，乙円，丙円が内接している。更に弦にも接している。甲円，乙円，丙円の径をそれぞれ 4 寸，13 寸，8 寸 4 分 5 厘のとき，外円の径を求めよ。

図のように記号を定め，方程式を立てる。
外円の半径を R とする。甲円，乙円，丙円の中心座標と半径を (0, 400 - R, 400)，(x4, y4, 1300), (x, y, 845) とする。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms x1, x3, y3, x4, y4, x, y, R, r, x0, y0;

eq1 = x^2 + y^2 - (R - 845)^2
eq2 = distance(x1, (800 - R), x3, y3, x4, y4) - 1300^2
eq3 = (x - x0)^2 + (y - y0)^2 - 845^2
eq4 = distance(x1, (800 - R), x3, y3, x, y) - 845^2;
eq5 = R + y4 - 2100
eq6 = x4^2 + y4^2 - (R - 1300)^2
eq7 = x3^2 + y3^2 - R^2
eq8 = x1^2 + (R - 800)^2 - R^2
eq9 = x0 - (x1 + x3) / 2
eq10 = y0 - ((800 - R) + y3) / 2;

方程式は solve() では解けないので，nlsolve() で数値解を求める。

eq1 |> println
eq2 |>  println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println
eq8 |> println
eq9 |> println
eq10 |> println

   x^2 + y^2 - (R - 845)^2
   (x4 - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y4 - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 1690000
   (x - x0)^2 + (y - y0)^2 - 714025
   (x - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 714025
   R + y4 - 2100
   x4^2 + y4^2 - (R - 1300)^2
   -R^2 + x3^2 + y3^2
   -R^2 + x1^2 + (R - 800)^2
   x0 - x1/2 - x3/2
   R/2 + y0 - y3/2 - 400

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10])

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, x3, y3, x, y, R, x4, y4, x0, y0) = u
   return [
       x^2 + y^2 - (R - 845)^2,
       (x4 - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y4 - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 1690000,
       (x - x0)^2 + (y - y0)^2 - 714025,
       (x - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 714025,
       R + y4 - 2100,
       x4^2 + y4^2 - (R - 1300)^2,
       -R^2 + x3^2 + y3^2,
       -R^2 + x1^2 + (R - 800)^2,
       x0 - x1/2 - x3/2,
       R/2 + y0 - y3/2 - 400,
   ]
end;

iniv = [-1500.0, 900, 2000, -1000, 800, 2200, 900, -100, -400, 300]
res = nls(H, ini=iniv)
println(res)

   ([-1700.0, 873.9999999999999, 2025.75, -1089.0000042590011, 816.7499943213318, 2206.25, 900.0, -106.24999999999997, -413.00000000000006, 309.74999999999994], true)

外円の半径は 2206.25 である。元の単位での径（直径）は 22 寸 0 分 6 厘 2 毛 5 糸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, x3, y3, x, y, R, x4, y4, x0, y0) = res[1]
   println("R = $R")
   plot()
   circle(0, 0, R)
   circle(x, y, 845, :blue)
   circle(x4, y4, 1300, :magenta)
   circle(0, 400-R, 400, :green)
   segment(x1, 800 - R, -x1, 800 - R)
   segment(x1, 800 - R, x3, y3)
    if more
       point(x4, y4, "(x4,y4)")
       point(0, 400-R, "400-R")
       point(x, y, "(x,y)")
       point(x0, y0, "(x0,y0)")
       point(x1, 800 - R, "(x1,800-R)")
       point(x3, y3, "(x3,y3)")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
    end
end;

算額（その107）

算額（その107）

埼玉県行田市 河原神社 万延2年(1861)3月
http://www.wasan.jp/saitama/kawara.html

https://yamabukiwasan.sakura.ne.jp/ymbk27.pdf

正方形の一辺の長さが 31 寸，甲円の直径が 18 寸のとき，乙円の直径を求めなさい。

乙円の半径を r1 とし，図のように記号を定め甲円，乙円が外接するという方程式を解く。

using SymPy

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

eq = (31//2 - r1)^2 + r1^2 - (r1 + 9)^2
solve(eq)[1] |> println

   7/2

乙円の半径は 7/2 である。もとの単位での直径は 7 寸である。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 7/2
   plot()
   circle(0, 0, 9)
   circle(31/2 - r1, r1, r1, :blue)
   circle(31/2 - r1, -r1, r1, :blue)
   circle(-(31/2 - r1), r1, r1, :blue)
   circle(-(31/2 - r1), -r1, r1, :blue)

   circle(r1, 31/2 - r1, r1, :blue)
   circle(-r1, 31/2 - r1, r1, :blue)
   circle(r1, -(31/2 - r1), r1, :blue)
   circle(-r1, -(31/2 - r1), r1, :blue)

   a = 31/2
   plot!([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:black, lw=0.5)
    if more
       point(9, 0, "9 ", :red, :bottom, :right)
       point(31/2 - r1, r1, "(31/2-r1,r1)", :blue, :center)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
    end
end;

 

算額（その106）

算額（その106）

埼玉県行田市 河原神社 万延2年(1861)3月
http://www.wasan.jp/saitama/kawara.html

https://yamabukiwasan.sakura.ne.jp/ymbk27.pdf

菱形の中に大小三つの円があります。菱形の長軸が20寸、短軸が15寸の時、小円の直径は何寸ですか？

大円，小円の半径を r1, r2 とする。下図のように，三角形 OAB, AaO, Abc が相似であることを使えば，SymPy など使わなくても解ける。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))
end;

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

原点から直線 AB までの距離が r1 である。

eq1 = distance(0, 15/2, 20/2, 0, 0, 0) - r1;
solve(eq1)[1] |> println

   6

同じように「点 c から直線 AB までの距離が r2 である」を解こうとすると，'NotImplementedError' になる。

eq2 = distance(0, 15/2, 20/2, 0, r1 + r2, 0) - r2;
eq2 |> println
1.　solve(eq2(r1 => 6))

   -r2 + sqrt((9*r1/25 + 9*r2/25 - 18/5)^2 + (12*r1/25 + 12*r2/25 - 24/5)^2)

原因は sqrt() の中が複雑なためのようだ。

二乗距離にすれば解ける。

eq3 = distance(0, 15/2, 20/2, 0, r1 + r2, 0)^2 - r2^2
eq3 |> println

   -r2^2 + (9*r1/25 + 9*r2/25 - 18/5)^2 + (12*r1/25 + 12*r2/25 - 24/5)^2

小円の半径は 3/2 である。もとの単位での直径は 3 寸である。

solve(eq3(r1 => 6))[1] |> println

   3/2

@syms r1::positive, r2::positive;
solve([eq1, eq3], (r1, r2))

   1-element Vector{Tuple{Sym, Sym}}:
    (6, 3/2)

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (6, 3/2)
   plot()
   circle(0, 0, r1)
   circle(15/2, 0, 3/2, :blue)
   circle(-15/2, 0, 3/2, :blue)
   plot!([20/2, 0, -20/2, 0, 20/2], [0, 15/2, 0, -15/2, 0], color=:black, lw=0.5)
    if more
       cosine = 4/5
       sine = 3/5
       (ax, ay) = (20/2*sine^2, 20/2*cosine*sine)
       (bx, by) = ((20/2 - r1 - r2)sine^2 + r1 + r2, (20/2 - r1 - r2)cosine*sine)
       point(ax, ay, " a", :red, :left, :bottom) 
       point(bx, by, " b", :blue, :left, :bottom) 
       point(r1 + r2, 0, "c", :blue)
       point(0, 0, " O", :black)
       point(20/2, 0, " A", :black)
       point(0, 15/2, " B", :black, :left, :bottom)
       segment(ax, ay, 0, 0)
       segment(bx, by, r1 + r2, 0)
       point(ax/2, ay/2, "r1", :red, mark=false)
       point((bx + r1 + r2)/2, by/2, "r2", :blue, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
    end
end;

算額（その105）

算額（その105）

福島県船引町 蚕養国神社 明治24年(1891)2月
http://www.wasan.jp/fukusima/kogaikuni.html

直線の上に 2 種類の円が 3 個ずつある。下の小円の中心間の距離（子）を 1 とすると，大円の径はいくつか。

図のように記号を定め，方程式を解く。方程式中に分数が出てくると計算時間がかかりすぎるので，子を 2 とし，各円の半径について方程式を立てる。

using SymPy
@syms r1::positive, r2::positive, y1::positive;

eq1 = r1^2 + (2r1 + r2 - y1)^2 - (r1 + r2)^2;
eq2 = (r1 - 1)^2 + (y1 - r2)^2 - (r1 + r2)^2;
eq3 = 1 + (r1 - r2)^2 - (r1 + r2)^2;

solve([eq1, eq2, eq3], (r1, r2, y1))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (1, 1/4, 3/2)

大円の半径は 1，直径は 2 である。子と大円の直径が同じである。子が 1 なら，大円の直径は 1 である。

算額では答えとして，「大円径4寸」とあるが，その左の「術」では「子が大円の径」と書いてある。「術」のほうが正しい。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, y1) = (1, 1/4, 3/2)
   plot()
   circle(0, r1, r1, :blue)
   circle(-r1, y1, r1, :blue)
   circle(r1, y1, r1, :blue)
   circle(0, 2r1 + r2, r2)
   circle(1, r2, r2)
   circle(-1, r2, r2)
   hline!([0], color=:black, lw=0.5)
    if more
       point(0, r1, " r1", :blue)
       point(r1, y1, "(r1,y1)", :blue)
       point(0, y1, " y1", :blue)
       point(0, 2r1 + r2, " 2r1+r2", :red)
       point(1, r2, " (1,r2)", :red)
       vline!([0], color=:black, lw=0.5)
    end
end;

 

算額（その104）

算額（その104）

愛媛県四国中央市川之江町 大岡八幡神社 文化元年(1804)4月
http://www.wasan.jp/ehime/kawanoe2.html

外円の中に甲円，乙円が1個ずつ，丙円，丁円，戊円が2つずつ互いに接して入っている。
甲円，乙円の直径をそれぞれ66寸，44寸としたとき，戊円の直径はいくつか。

甲円，乙円の半径をそれぞれ 33，22 とする。丙円，丁円，戊円の中心座標と半径を (x1, y1, r1), (x2, y2, r2), (x3, y3, r3) として，以下の9本の方程式を解く。

条件式と未知数の数が等しいので，解けるはずだが solve() では解けない。

using SymPy

@syms x1::positive, y1::positive, r1::positive, x2::positive, y2::positive, r2::positive, x3::positive, y3::positive, r3::positive;
@syms x1, y1, r1, x2, y2, r2, x3, y3, r3
eq1 = x3^2 + (y3 - 33)^2 - (33 + r3)^2 |> expand
eq2 = x1^2 + (y1 - 33)^2 - (33 + r1)^2 |> expand
eq3 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2 |> expand
eq4 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2 |> expand
eq5 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand
eq6 = x2^2 + (y2 - 88)^2 - (22 + r2)^2 |> expand
eq7 = x3^2 + (88 - y3)^2 - (22 + r3)^2 |> expand
eq8 = x2^2 + (y2 - 55)^2 - (55 - r2)^2 |> expand
eq9 = x1^2 + (y1 - 55)^2 - (55 - r1)^2 |> expand
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9]);

今までの経験で，nlsolve() で解くことにする。

eq1 |> println
eq2 |> println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println
eq8 |> println
eq9 |> println

   -r3^2 - 66*r3 + x3^2 + y3^2 - 66*y3
   -r1^2 - 66*r1 + x1^2 + y1^2 - 66*y1
   -r1^2 - 2*r1*r3 - r3^2 + x1^2 - 2*x1*x3 + x3^2 + y1^2 - 2*y1*y3 + y3^2
   -r1^2 - 2*r1*r2 - r2^2 + x1^2 - 2*x1*x2 + x2^2 + y1^2 - 2*y1*y2 + y2^2
   -r2^2 - 2*r2*r3 - r3^2 + x2^2 - 2*x2*x3 + x3^2 + y2^2 - 2*y2*y3 + y3^2
   -r2^2 - 44*r2 + x2^2 + y2^2 - 176*y2 + 7260
   -r3^2 - 44*r3 + x3^2 + y3^2 - 176*y3 + 7260
   -r2^2 + 110*r2 + x2^2 + y2^2 - 110*y2
   -r1^2 + 110*r1 + x1^2 + y1^2 - 110*y1

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, y2, r2, x3, y3, r3) = u
return [
       -r3^2 - 66*r3 + x3^2 + y3^2 - 66*y3,
       -r1^2 - 66*r1 + x1^2 + y1^2 - 66*y1,
       -r1^2 - 2*r1*r3 - r3^2 + x1^2 - 2*x1*x3 + x3^2 + y1^2 - 2*y1*y3 + y3^2,
       -r1^2 - 2*r1*r2 - r2^2 + x1^2 - 2*x1*x2 + x2^2 + y1^2 - 2*y1*y2 + y2^2,
       -r2^2 - 2*r2*r3 - r3^2 + x2^2 - 2*x2*x3 + x3^2 + y2^2 - 2*y2*y3 + y3^2,
       -r2^2 - 44*r2 + x2^2 + y2^2 - 176*y2 + 7260,
       -r3^2 - 44*r3 + x3^2 + y3^2 - 176*y3 + 7260,
       -r2^2 + 110*r2 + x2^2 + y2^2 - 110*y2,
       -r1^2 + 110*r1 + x1^2 + y1^2 - 110*y1
     ];
end;

iniv = [40.0, 57, 14,
       33,    82, 12,
       20,    68, 8];

res = nls(H, ini=iniv)

   ([40.58178048548881, 57.391304347826086, 14.347826086956518, 33.33503397022296, 82.5, 11.785714285714292, 21.21320343559642, 67.5, 7.499999999999996], false)

収束していないが，ftol=1e-14 が厳しすぎるのである。図に描いてみると十分な精度である。
戊円の直径は 7.5*2 = 15 である。元の単位で言えば 15 寸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, y1, r1, x2, y2, r2, x3, y3, r3) = [40.58178048548881, 57.391304347826086, 14.347826086956518, 33.33503397022296, 82.5, 11.785714285714292, 21.21320343559642, 67.5, 7.499999999999996]
   plot()
   circle(0, 55, 55)
   circle(0, 33, 33)
   circle(0, 88, 22)
   circle(x1, y1, r1, :green)
   circle(-x1, y1, r1, :green)
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   circle(x3, y3, r3, :magenta)
   circle(-x3, y3, r3, :magenta)
    if more
       a = @sprintf("丙円 (%.1f, %.1f, %.1f)\n", x1, y1, r1)
       b = @sprintf("丁円 (%.1f, %.1f, %.1f)\n", x2, y2, r2)
       c = @sprintf("戊円 (%.1f, %.1f, %.1f)", x3, y3, r3)
       point(-20, 22, a * b * c, :black, mark=false)
       point(0, 33, " 甲 (0,33,33)", :red)
       point(0, 88, " 乙 (0,88,22)", :red)
       point(x1, y1, "    丙\n    (x1,y1;r1)", :green, :top)
       point(x2, y2, "    丁\n    (x2,y2;r2)", :blue, :top)
       point(x3, y3, "    戊\n    (x3,y3;r3)", :magenta, :top)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
    end
end;

 

算額（その103）

算額（その103）

七〇 加須市大字外野 棘脱地蔵堂 明治6年(1873)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

東京都中央区 福徳神社 令和2年(2020)
http://www.wasan.jp/tokyo/hukutoku.html

大円の中に斜線を隔てて甲円 2 個，乙円 2 個，丙円 1 個を入れる。乙円の直径が1のとき，丙円の直径はいかほどか。

大円の半径を R として，以下のように記号を定める。

⊿ABC と ⊿ADE が相似で，甲円の半径 DE = 1 なので，丙円の半径 BC は DE/2=0.5。∠BAC を θ とすれば，cosθ = √8/3, sinθ = sqrt(1-cosθ^2) により，B，D，F の座標を求めることができる。

include("julia-source.txt")

using SymPy

@syms R, r1, x2, y2, r2

R = 2r1
r3 = r1/2  # ⊿ABC と ⊿ADE が相似，DE = 1
AD = sqrt(Sym(8))  # sqrt(AE^2 - ED^2) = sqrt(3^2 - 1)
cosθ = sqrt(Sym(8))/3  # AD/3, θ = ∠BAC
sinθ = sqrt(1 - cosθ^2)
Dx = AD * sinθ
Dy = R - AD * cosθ;

次に，乙円の中心座標を (x2, y2)，半径を r2 として以下の方程式を立て解く。

eq1 = x2^2 + (R - r1 - y2)^2 - (r1+ r2)^2 |> expand
eq2 = x2^2 + y2^2 - (R - r2)^2 |> expand
eq3 = dist2(0, R, Dx, Dy, x2, y2, r2)
res = solve([eq1, eq2, eq3], (r1, x2, y2))[4]  # 4 of 4

   (25*r2/16, 3*sqrt(2)*r2/2, r2/8)

甲円の半径 r1 は乙円の半径 r2 の 25/16 倍である。
さらに，丙円の半径は乙円の半径の 1/2 である。
よって，丙円の直径は 乙円の直径の 25/32 倍である。

乙円の直径が 1 寸のとき，丙円の直径は 0.78125 寸である。

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

   x2 = 1.06066;  y2 = 0.0625
   r2 = 0.5;  r3 = 0.390625
   Bx = 0.368285;  By = 0.520833;  Dx = 0.73657;  Dy = -0.520833

function draw(r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x2, y2) = (25*r2/16, 3*sqrt(2)*r2/2, r2/8)
   R = 2r1
   r3 = r1/2  # ⊿ABC と ⊿ADE が相似，DE = 1
   AD = √8r1  # sqrt(AE^2 - ED^2) = sqrt(3^2 - 1)
   cosθ = √8/3  # AD/3, θ = ∠BAC
   sinθ = sqrt(1 - cosθ^2)
   AB = √2r1  # sqrt((R - r3)^2 -r3^2)
   Bx = AB * sinθ
   By = R - AB * cosθ
   Dx = AD * sinθ
   Dy = R - AD * cosθ
   AF = 4AD / 3                                                                         
   Fx = AF * sinθ
   Fy = R - AF * cosθ
   @printf("x2 = %g;  y2 = %g\n", x2, y2)
   @printf("r2 = %g;  r3 = %g\n", r2, r3)
   @printf("Bx = %g;  By = %g;  Dx = %g;  Dy = %g\n", Bx, By, Dx, Dy)
   plot()
   circle(0, 0, R)
   circle22(0, r1, r1)
   circle(0, r3, r3, :magenta)
   circle2(x2, y2, r2, :blue)
   segment(0, R, Fx, Fy)
   segment(0, R, -Fx, Fy)
    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)
       segment(0, R, 0, -R)
       segment(0, -R, Fx, Fy)
       segment(0, r3, Bx, By)
       segment(0, -r1, Dx, Dy)
       point(0, 0, "0 ", :red, :right, delta=-delta/2)
       point(0, R, "A ", :red, :right, :bottom, delta=delta/2)
       point(0, r3, "C ", :magenta, :right, :bottom)
       point(Bx, By, " B", :magenta)
       point(Dx, Dy, " D", :red)
       point(Fx, Fy, " F", :red)
       point(0, -R, " G", :red)
       point(0, -r1, "E ", :red, :right, :bottom)
       point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, :bottom, delta=delta/2)
       point(0, r1, "甲円:r1\n(0,r1)", :red, :center, delta=-delta/2)
       point(0, r3, "丙円:r3\n(0,r3)", :magenta, :center, delta=-delta/2)
    end
end;