https://sangaku0418.hatenablog.com/entry/2024/11/30/160614

福島市立子山篠葉沢稲荷明治10年(1877)

街角の数学 Street Wasan　街角の数学～落書き帳「○△□」～ 68．ウルトラマン・ピーチ
http://streetwasan.web.fc2.com/math15.8.19.html
キーワード：円4個，外円，斜線
#Julia #SymPy #算額 #和算 #数学

外円の中に 2 本の斜線（等斜）を引き，区画された領域に 4 個の等円を容れる。外円の直径が 13 寸のとき，等円の直径はいかほどか。

外円の半径と中心座標を \(R, (0, 0)\)
等円の半径と中心座標を \(r, (x_1, y_1), (r, y_2)\)
斜線と円の交点座標を\( \left (x_0, -\sqrt{r^2 - x_0^2}\right )\)
とおき，以下の連立方程式を解く。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy
@syms R::positive, r::positive, x1::positive, y1::positive, y2::positive, x0::positive
t = sqrt(R^2 - x0^2)
eq1 = x1^2 + y1^2 - (R - r)^2
eq2 = y1* (t + R) - x0*x1 # y1/x1 - x0/(t + R)
eq3 = dist2(0, R, x0, -t, x1, y1, r) |> numerator
eq4 = dist2(0, R, x0, -t, r, y2, r) |> numerator
eq5 = sqrt(x0^2 + (R - t)^2) + sqrt(x0^2 + (R + t)^2) - 2R - 2r;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r, x1, y1, y2, x0))

function H(u)
(r, x1, y1, y2, x0) = u
t = sqrt(R^2 - x0^2)
return [
x1^2 + y1^2 - (R - r)^2, # eq1
-x0*x1 + y1*(R + sqrt(R^2 - x0^2)), # eq2
-4*R^3*r^2 + R^3*x0^2 - 4*R^3*x0*x1 + 4*R^3*x1^2 - 4*R^2*r^2*sqrt(R^2 - x0^2) - 2*R^2*x0^2*y1 + R^2*x0^2*sqrt(R^2 - x0^2) + 4*R^2*x0*x1*y1 - 4*R^2*x0*x1*sqrt(R^2 - x0^2) + 4*R^2*x1^2*sqrt(R^2 - x0^2) + 2*R*r^2*x0^2 + 2*R*x0^3*x1 - 3*R*x0^2*x1^2 + R*x0^2*y1^2 - 2*R*x0^2*y1*sqrt(R^2 - x0^2) + 4*R*x0*x1*y1*sqrt(R^2 - x0^2) - 2*x0^3*x1*y1 - x0^2*x1^2*sqrt(R^2 - x0^2) + x0^2*y1^2*sqrt(R^2 - x0^2), # eq3
x0*(-4*R^3*r + R^3*x0 + 4*R^2*r*y2 - 4*R^2*r*sqrt(R^2 - x0^2) - 2*R^2*x0*y2 + R^2*x0*sqrt(R^2 - x0^2) - R*r^2*x0 + 2*R*r*x0^2 + 4*R*r*y2*sqrt(R^2 - x0^2) + R*x0*y2^2 - 2*R*x0*y2*sqrt(R^2 - x0^2) - r^2*x0*sqrt(R^2 - x0^2) - 2*r*x0^2*y2 + x0*y2^2*sqrt(R^2 - x0^2)), # eq4
-2*R - 2*r + sqrt(x0^2 + (R - sqrt(R^2 - x0^2))^2) + sqrt(x0^2 + (R + sqrt(R^2 - x0^2))^2), # eq5
]
end;

R = 13/2
iniv = BigFloat[1.9, 4.1, 1.7, -3.4, 4.6]
res = nls(H, ini=iniv)

([2.0, 4.153846153846154, 1.7307692307692308, -3.5, 4.615384615384615], true)

外円の直径が 13 寸のとき，等円の直径は 4 寸である。

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

function draw(R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r, x1, y1, y2, x0) = [2.0, 4.153846153846154, 1.7307692307692308, -3.5, 4.615384615384615]
t = sqrt(R^2 - x0^2)
@printf("外円の直径が %g のとき，等円の直径は %g である。\n", 2R, 2r)
@printf("R = %g; r = %g; x1 = %g; y1 = %g; y2 = %g; x0 = %g\n", R, r, x1, y1, y2, x0)
plot()
circle(0, 0, R)
circle2(x1, y1, r, :blue)
circle2(r, y2, r, :blue)
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)
plot!([0, x0, 0, -x0, 0, 0], [R, -t, -R, -t, R, -R], color=:green, lw=0.7)
point(x1, y1, "等円:r,(x1,y1)", :blue, :center, delta=-delta/2)
point(r, y2, "等円:r,(r,y2)", :blue, :center, delta=-delta/2)
point(x0, -t, " (x0,sqrt(R^2-x0^2)", :green, :left, :vcenter)
point(0, R, "R", :red, :center, :bottom, delta=delta/2)
end
end;