算額（その2159） - 算額あれこれ

藤田貞資：精要算法（下巻）天明元年(1781)

和算の館
http://www.wasan.jp/seiyou/seiyou.html

キーワード：円3個，三角形，中鈎
#Julia #SymPy #算額 #和算 #数学

三角形の中に，全円と，中鈎を隔てて中円と小円を容れる。中円の直径が106寸2分，小円の直径が 88寸5分，また，中鈎と全円の直径の差と大斜の積は，長股と中斜の差と中鈎の積と等しい。このとき，大斜はいかほどか。

図形は算額（その0679）と同じであるが，条件と求めるものが違う。

三角形の大斜，中斜，中鈎，長股，短股を \(a,\ b,\ h，c, d\)
全円，中円，小円の半径を \(r_1,\ r_2,\ r_3\)
とおき，算法助術の公式57 を適用する。

一度に解けないので，eq1, eq2, eq3, eq4 から \(a,\ b,\ c,\ h\) を求め，eq5 に代入して更新した方程式を解いて \(r_1\) を求める。

include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r1, r2, r3, h, a, b, c
#(r2, r3) = (50//2, 48//2)
#(r2, r3) = (106.2/2, 88.5/2)
eq1 = h*(r2 + r3) - (2r2*r3 + r1*h)
eq2 = a*(h - 2r1) - h*(b - c)
eq3 = (c + h - b) - 2r2
eq4 = ( (a - c) + h - sqrt( (a - c)^2 + h^2)) - 2r3
eq5 = (a + b + sqrt( (a - c)^2 + h^2))*r1 - a*h
# √を外す
eq5 = (-a*h/r1 + a + b)^2 - (a^2 - 2*a*c + c^2 + h^2)

\(\displaystyle - a^{2} + 2 a c - c^{2} - h^{2} + \left(- \frac{a h}{r_{1}} + a + b\right)^{2}\)

eq5_2 = eq5 |> expand

\(\displaystyle \frac{a^{2} h^{2}}{r_{1}^{2}} - \frac{2 a^{2} h}{r_{1}} - \frac{2 a b h}{r_{1}} + 2 a b + 2 a c + b^{2} - c^{2} - h^{2}\)

\(- h r_{1} + h \left(r_{2} + r_{3}\right) - 2 r_{2} r_{3}\)

\(a \left(h - 2 r_{1}\right) - h \left(b - c\right)\)

\(- b + c + h - 2 r_{2}\)

\(a - c + h - 2 r_{3} - \sqrt{h^{2} + \left(a - c\right)^{2}}\)

\(\displaystyle - a^{2} + 2 a c - c^{2} - h^{2} + \left(- \frac{a h}{r_{1}} + a + b\right)^{2}\)

res = solve([eq1, eq2, eq3, eq4], (a, b, c, h))[1];

# a
ans_a = res[1]
@show(ans_a)

ans_a = 2*r2^2*r3/( (r1 - r3)*(-r1 + r2 + r3))

\(\displaystyle \frac{2 r_{2}^{2} r_{3}}{\left(r_{1} - r_{3}\right) \left(- r_{1} + r_{2} + r_{3}\right)}\)

# b
ans_b = res[2]
@show(ans_b)

ans_b = (-2*r1^2*r2 - r1^2*r3 + 2*r1*r2^2 + 2*r1*r2*r3 + 2*r1*r3^2 - 3*r2^2*r3 - r3^3)/(r1^2 - r1*r2 - 2*r1*r3 + r2*r3 + r3^2)

\(\displaystyle \frac{- 2 r_{1}^{2} r_{2} - r_{1}^{2} r_{3} + 2 r_{1} r_{2}^{2} + 2 r_{1} r_{2} r_{3} + 2 r_{1} r_{3}^{2} - 3 r_{2}^{2} r_{3} - r_{3}^{3}}{r_{1}^{2} - r_{1} r_{2} - 2 r_{1} r_{3} + r_{2} r_{3} + r_{3}^{2}}\)

# c
ans_c = res[3]
@show(ans_c)

ans_c = r3*(-r1^2 + 2*r1*r3 - r2^2 - r3^2)/(r1^2 - r1*r2 - 2*r1*r3 + r2*r3 + r3^2)

\(\displaystyle \frac{r_{3} \left(- r_{1}^{2} + 2 r_{1} r_{3} - r_{2}^{2} - r_{3}^{2}\right)}{r_{1}^{2} - r_{1} r_{2} - 2 r_{1} r_{3} + r_{2} r_{3} + r_{3}^{2}}\)

# h
ans_h = res[4]
@show(ans_h)

ans_h = 2*r2*r3/(-r1 + r2 + r3)

\(\displaystyle \frac{2 r_{2} r_{3}}{- r_{1} + r_{2} + r_{3}}\)

\(\left(r_{1} - r_{3}\right)^{3} \left(r_{1} - r_{2} - r_{3}\right)^{2} \left(r_{1}^{4} - 2 r_{1}^{3} r_{2} - 2 r_{1}^{3} r_{3} + r_{1}^{2} r_{2}^{2} + 4 r_{1}^{2} r_{2} r_{3} + r_{1}^{2} r_{3}^{2} - 2 r_{1} r_{2}^{2} r_{3} - 2 r_{1} r_{2} r_{3}^{2} + 2 r_{2}^{2} r_{3}^{2}\right) \left(r_{1}^{3} r_{2} + r_{1}^{3} r_{3} - 2 r_{1}^{2} r_{2}^{2} - 2 r_{1}^{2} r_{2} r_{3} - 2 r_{1}^{2} r_{3}^{2} + r_{1} r_{2}^{3} + 3 r_{1} r_{2}^{2} r_{3} + r_{1} r_{2} r_{3}^{2} + r_{1} r_{3}^{3} - 2 r_{2}^{3} r_{3}\right)\)

eq5_2 は 4 個の項に因数分解されるが，有効な方程式は以下の \(r_1\) に関する 3 次式である。

ParseError:
eq5_2 は 4 個の項に因数分解されるが，有効な方程式は以下の \(r_1\) に関する 3 次式である。
#└────────────────────────────────────────────────────────────────────────────────┘ ── extra tokens after end of expression

Stacktrace:

[1] top-level scope

@ In[10]:1

temp = eq5_3.args[4]

\(r_{1}^{3} r_{2} + r_{1}^{3} r_{3} - 2 r_{1}^{2} r_{2}^{2} - 2 r_{1}^{2} r_{2} r_{3} - 2 r_{1}^{2} r_{3}^{2} + r_{1} r_{2}^{3} + 3 r_{1} r_{2}^{2} r_{3} + r_{1} r_{2} r_{3}^{2} + r_{1} r_{3}^{3} - 2 r_{2}^{3} r_{3}\)

ans_r1 = solve(temp, r1)[1]

\(\displaystyle - \frac{- \frac{3 \left(r_{2}^{3} + 3 r_{2}^{2} r_{3} + r_{2} r_{3}^{2} + r_{3}^{3}\right)}{r_{2} + r_{3}} + \frac{\left(- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}\right)^{2}}{\left(r_{2} + r_{3}\right)^{2}}}{3 \sqrt[3]{- \frac{27 r_{2}^{3} r_{3}}{r_{2} + r_{3}} + \frac{\sqrt{- 4 \left(- \frac{3 \left(r_{2}^{3} + 3 r_{2}^{2} r_{3} + r_{2} r_{3}^{2} + r_{3}^{3}\right)}{r_{2} + r_{3}} + \frac{\left(- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}\right)^{2}}{\left(r_{2} + r_{3}\right)^{2}}\right)^{3} + \left(- \frac{54 r_{2}^{3} r_{3}}{r_{2} + r_{3}} - \frac{9 \left(- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}\right) \left(r_{2}^{3} + 3 r_{2}^{2} r_{3} + r_{2} r_{3}^{2} + r_{3}^{3}\right)}{\left(r_{2} + r_{3}\right)^{2}} + \frac{2 \left(- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}\right)^{3}}{\left(r_{2} + r_{3}\right)^{3}}\right)^{2}}}{2} - \frac{9 \left(- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}\right) \left(r_{2}^{3} + 3 r_{2}^{2} r_{3} + r_{2} r_{3}^{2} + r_{3}^{3}\right)}{2 \left(r_{2} + r_{3}\right)^{2}} + \frac{\left(- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}\right)^{3}}{\left(r_{2} + r_{3}\right)^{3}}}} - \frac{\sqrt[3]{- \frac{27 r_{2}^{3} r_{3}}{r_{2} + r_{3}} + \frac{\sqrt{- 4 \left(- \frac{3 \left(r_{2}^{3} + 3 r_{2}^{2} r_{3} + r_{2} r_{3}^{2} + r_{3}^{3}\right)}{r_{2} + r_{3}} + \frac{\left(- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}\right)^{2}}{\left(r_{2} + r_{3}\right)^{2}}\right)^{3} + \left(- \frac{54 r_{2}^{3} r_{3}}{r_{2} + r_{3}} - \frac{9 \left(- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}\right) \left(r_{2}^{3} + 3 r_{2}^{2} r_{3} + r_{2} r_{3}^{2} + r_{3}^{3}\right)}{\left(r_{2} + r_{3}\right)^{2}} + \frac{2 \left(- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}\right)^{3}}{\left(r_{2} + r_{3}\right)^{3}}\right)^{2}}}{2} - \frac{9 \left(- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}\right) \left(r_{2}^{3} + 3 r_{2}^{2} r_{3} + r_{2} r_{3}^{2} + r_{3}^{3}\right)}{2 \left(r_{2} + r_{3}\right)^{2}} + \frac{\left(- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}\right)^{3}}{\left(r_{2} + r_{3}\right)^{3}}}}{3} - \frac{- 2 r_{2}^{2} - 2 r_{2} r_{3} - 2 r_{3}^{2}}{3 \left(r_{2} + r_{3}\right)}\)

ans_r1(r2 => 106.2/2, r3 => 88.5/2)

\(67.6664037158923\)

# r1 は長い式であるが，共通項を整理するなどして短く表現できる
# ans_r1 |> println

x = (r2 + r3)
u = (-2*r2^2 - 2*r2*r3 - 2*r3^2)
s = u^2/x^2
t = u^3/x^3
w = (r2^3 + 3*r2^2*r3 + r2*r3^2 + r3^3)
v = (-27*r2^3*r3/x + sqrt(-4*(-3*w/x + s)^3 + (-54*r2^3*r3/x - 9*u*w/x^2 + 2*t)^2)/2 - 9*u*w/(2*x^2) + t)
ans_r1 = -(-3*w/x + s)/(3*v^(1/3)) - v^(1/3)/3 - u/(3*x);

\(r_2,\ r_3\) に実値を代入すれば \(r_1\)，更に\(r_1\) に代入して \(a\) を求めることができる。結果は「答曰」に述べられているものに一致する。

ans_r1(r2 => 1062//20, r3 => 885//20).evalf()

\(67.6664037158916\)

ans_a(r1 => ans_r1)(r2 => 1062//20, r3 => 885//20).evalf()

\(359.00095142033\)

本問題では，\(r_1\) の式は SymPy で解く限りは簡単ではなく，解析解を求める意味はあまりない。数値解は以下のように簡単に求めることができる。

include("julia-source.txt"); # julia-source.txt ソース
function driver(r2, r3)
function H(u)
function parameters()
eq1 = h*(r2 + r3) - (2r2*r3 + r1*h)
eq2 = a*(h - 2r1) - h*(b - c)
eq3 = (c + h - b) - 2r2
eq4 = ( (a - c) + h - sqrt( (a - c)^2 + h^2)) - 2r3
eq5 = (a + b + sqrt( (a - c)^2 + h^2))*r1 - a*h
return [eq1, eq2, eq3, eq4, eq5]
end;
(a, b, c, h, r1) = u
return parameters()
end;
iniv = BigFloat[105, 85, 64, 55, 23]
iniv = BigFloat[6.22, 4.94, 3.87, 3.07, 1.27].*r2
res = nls(H, ini=iniv)
res[2] || println("収束していない")
return res[1]
end;
driver(106.2/2, 88.5/2) |> println

[359.00095142033024, 266.522543869703, 214.40783176334432, 158.31471210635868, 67.66640371589159]

描画関数プログラムのソースを見る

function draw(r2, r3, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b, c, h, r1) = driver(r2, r3)
#(a, b, c, h, r1) = (359.000951420330, 266.522543869703, 214.407831763345, 158.314712106359, 67.6664037158916)
x = (r2 + r3)
u = (-2*r2^2 - 2*r2*r3 - 2*r3^2)
s = u^2/x^2
t = u^3/x^3
w = (r2^3 + 3*r2^2*r3 + r2*r3^2 + r3^3)
v = (-27*r2^3*r3/x + sqrt(-4*(-3*w/x + s)^3 + (-54*r2^3*r3/x - 9*u*w/x^2 + 2*t)^2)/2 - 9*u*w/(2*x^2) + t)
r1 = -(-3*w/x + s)/(3*v^(1/3)) - v^(1/3)/3 - u/(3*x)
a = 2*r2^2*r3/( (r1 - r3)*(-r1 + r2 + r3))
b1 = -2*r1^2*r2 - r1^2*r3 + 2*r1*r2^2 + 2*r1*r2*r3 + 2*r1*r3^2 - 3*r2^2*r3 - r3^3
b2 = r1^2 - r1*r2 - 2*r1*r3 + r2*r3 + r3^2
b = b1/b2
c = r3*(-r1^2 + 2*r1*r3 - r2^2 - r3^2)/(r1^2 - r1*r2 - 2*r1*r3 + r2*r3 + r3^2)
h = 2*r2*r3/(-r1 + r2 + r3)
plot([c, c, 0, a, c, 0], [0, h, 0, 0, h, 0], color=:green, lw=0.5)
circle(c + r3, r3, r3)
circle(c - r2, r2, r2, :blue)
circle(c - (r2 - r3), r1, r1, :magenta)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
vline!([0], color=:gray80, lw=0.5)
hline!([0], color=:gray80, lw=0.5)
point(c + r3, r3, "小円:r3\n(c+r3,r3)", :red, :center, delta=-delta)
point(c - r2, r2, "中円:r2\n(c-r2,r2)", :blue, :center, delta=-delta)
point(c - (r2 - r3), r1, "大円:r1\n(c-r2+r3,r1)", :magenta, :center, :bottom, delta=delta)
point(0, 0, "(0,0)", :green, :center, delta=-delta)
point(c, 0, "(c,0)", :green, :center, delta=-delta)
point(a, 0, "(a,0)", :green, :center, delta=-delta)
point(c, h, "(c,h)", :green, :center, :bottom, delta=delta)
point(a/2, 0, "大斜", :black, :center, delta=-delta, mark=false)
point(c/2, h/2, "中斜", :black, :right, :vcenter, deltax=-7delta, mark=false)
point( (a + c)/2, h/2, "小斜", :black, :left, :vcenter, deltax=4delta, mark=false)
point(c, 0.7h, "中鈎", :black, :left, deltax=delta, mark=false)
plot!(xlims=(-7delta, a + 5delta), ylims=(-8delta, h + 3delta))
end
end;
draw(106.2/2, 88.5/2, true)
#draw(50/2, 48/2, true)