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

福島県白河市明神境明神万延元年(1860)

和算の館
http://wasan.jp//fukusima/sakai1.html
http://wasan.jp//fukusima/sakai1-2.png
http://wasan.jp//fukusima/sakai1kaisetu.png
キーワード：円4個，半円，円弧，二等辺三角形
#Julia #SymPy #算額 #和算 #数学

2 個の半円が交わり，中に二等辺三角形と円弧があり，等円 4 個が接している。「子」が 1 寸のとき，「矢」はいかほどか。

半円の半径と中心座標を \(s,\ (0,\ 0),\ (0,\ s)\)
円弧の半径と中心座標を \(R,\ (0,\ y)\)
等円の半径と中心座標を \(r,\ (x_1,\ y_1),\ (x_2,\ s - r)\)
とおき，以下の連立方程式を解く。

SymPy の能力的には eq1〜eq6 をまとめて解くことが難しい。後半に示す数値解を吟味すると，\(R = (2 + \sqrt{2})子\) であることがわかる。

include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms 子::positive, R::positive, x1::positive, y1::positive, x2::positive, r::positive, y::positive
s = 子/(√Sym(2) - 1)
R = (2 + √Sym(2))*子;

その結果 eq6 を解くことができ，\(y = 2(1 + \sqrt{2})子\) であることもわかる。

eq6 = (y - s)^2 + s^2 - R^2
ans_y = solve(eq6, y)[1] |> simplify
@show(ans_y)

ans_y = 2*子*(1 + sqrt(2))

\(2 子 \left(1 + \sqrt{2}\right)\)

\(R,\ y\) がわかれば，「矢」を求めることができる。

矢 = R - (ans_y - s) |> simplify

\(子\)

「矢」は「子」に等しい。

図を描くために必要なすべてのパラメータは，eq1〜eq4 を解くことにより得られる。

using SymPy
@syms 子::positive, R::positive, x1::positive, y1::positive, x2::positive, r::positive, y::positive
s = 子/(√Sym(2) - 1)
R = (2 + √Sym(2))*子
y = 2子*(1 + √Sym(2))
eq1 = x1^2 + y1^2 - (s + r)^2
eq2 = x1^2 + (s - y1)^2 - (s - r)^2
eq3 = x2^2 + (s - r)^2 - (s +r)^2
eq4 = x1^2 + (y - y1)^2 - (R + r)^2
eq5 = x2^2 + (y - (s - r))^2 - (R - r)^2;
res = solve([eq1, eq2, eq3, eq4], (r, x1, y1, x2))[1]

(子/14 + sqrt(2)*子/7, sqrt(141*sqrt(2)*子^2/98 + 207*子^2/98), 9*子/14 + 11*sqrt(2)*子/14, sqrt(6*sqrt(2)*子^2/7 + 10*子^2/7))

# r
ans_r = res[1] |> simplify
@show(ans_r)

ans_r = 子*(1 + 2*sqrt(2))/14

\(\displaystyle \frac{子 \left(1 + 2 \sqrt{2}\right)}{14}\)

ans_r(子 => 1).evalf()

\(0.273459080339014\)

# x1
ans_x1 = res[2] |> simplify
@show(ans_x1)

ans_x1 = 子*sqrt(282*sqrt(2) + 414)/14

\(\displaystyle \frac{子 \sqrt{282 \sqrt{2} + 414}}{14}\)

ans_x1(子 => 1).evalf()

\(2.03641369512682\)

# y1
ans_y1 = res[3] |> simplify
@show(ans_y1)

ans_y1 = 子*(9 + 11*sqrt(2))/14

\(\displaystyle \frac{子 \left(9 + 11 \sqrt{2}\right)}{14}\)

ans_y1(子 => 1).evalf()

\(1.75402494186457\)

# y2
ans_y2 = res[4] |> simplify

\(\displaystyle \frac{子 \sqrt{42 \sqrt{2} + 70}}{7}\)

ans_y2(子 => 1).evalf()

\(1.62503984013749\)

矢(子 => 1).evalf()

\(1.0\)

以下に，図を描画するすべてのパラメータの数値解を求める例を示す。

function driver(子)
function H(u)
function parameters()
eq1 = x1^2 + y1^2 - (s + r)^2
eq2 = x1^2 + (s - y1)^2 - (s - r)^2
eq3 = x2^2 + (s - r)^2 - (s +r)^2
eq4 = x1^2 + (y - y1)^2 - (R + r)^2
eq5 = x2^2 + (y - (s - r))^2 - (R - r)^2
eq6 = (y - s)^2 + s^2 - R^2
return [eq1, eq2, eq3, eq4, eq5, eq6]
end;
(R, y, r, x1, y1, x2) = u
return parameters()
end;
s = 子/(√2 - 1)
iniv = s .* BigFloat[3.414, 4.828, 0.273, 2.036, 1.754, 1.625]
res = nls(H, ini=iniv)
res[2] || println("収束していない")
return res[1]
end;
println(driver(1))

[3.414213562373094, 4.828427124746189, 0.2734590803390135, 2.036413695126819, 1.7540249418645744, 1.6250398401374904]

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

function draw(子, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", showaxis=false, fontfamily="IPAMincho")
s = 子/(√2 - 1)
(R, y, r, x1, y1, x2) = driver(子)
# (R, y, r, x1, y1, x2) = 子 .* [3.414213562373094, 4.828427124746189, 0.2734590803390135, 2.036413695126819, 1.7540249418645744, 1.6250398401374904]
R = (2 + √2)*子
println( (R, y, r, x1, y1, x2))
矢 = s - (y - R)
@printf("子 = %g のとき，矢 = %g\n", 子, 矢)
θ = asind( (y - s)/R)
plot([s, s, -s, -s, s], [0, s, s, 0, 0], color=:green, lw=0.5)
plot!([-s, 0, s], [0, s, 0], color=:green, lw=0.5)
circle(0, y, R, beginangle=270 - θ, endangle=270 + θ)
circle(0, 0, s, :blue, beginangle=0, endangle=180)
circle(0, s, s, :blue, beginangle=180, endangle=360)
circle2(x1, y1, r, :magenta)
circle2(x2, s - r, r, :magenta)
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(s, s, "(s,s)", :green, :right, :bottom, delta=delta)
point(x1, y1, "r,(x1,y1)", :black, :center, delta=-delta/2, deltax=-2delta)
point(x2, s - r, "r,(x2,s-r)", :black, :center, :bottom, delta=delta, deltax=4delta)
point(0, 0, "0", :black, :center, :bottom, delta=delta)
end
dimension_line(0, s, 0, y - R, " 矢", :black, :left, lw=1, length=s/20)
dimension_line(s, 0, s - 子/√2, 子/√2, " 子", :black, :left, lw=1, length=s/20)
end;
draw(1, true)