算額(その904)
一〇〇 桶川市小針領家 氷川諏訪神社 明治30年(1897)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.
外円内に甲乙丙の 3 円と斜線(弦)が入っている。外円,甲円の直径がそれぞれ 3 寸,2 寸のとき,乙円の直径はいかほどか。
注:外円,甲円,丙円は一直線上にあり,弦は外円と甲円の接点を通る。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (0, r3 - R)
弦と外円の交点座標を (x, y)
とおき,以下の連立方程式を解く。
SymPy の性能上,一度に解くことができないので,逐次的に解いてゆく。
乙円の半径は eq1 だけで求めることができる。
include("julia-source.txt");
using SymPy
@syms R::positive, r1::positive,
r2::positive, x2::positive, y2::positive,
r3::positive, x::positive, y::positive
r3 = R - r1
eq1 = r3/(2R - r3) - (R - 2r2)/R
ans_r2 = solve(eq1, r2)[1]
ans_r2 |> println
ans_r2(R => 3//2, r1 => 2/2)*2 |> println
R*r1/(R + r1)
1.20000000000000
乙円の半径は 外円,甲円の直径に基づいて計算できる。術の通り,外円と甲円の直径の積を外円と甲円の直径の和で割ればよい。
外円,甲円の直径がそれぞれ 3 寸,2 寸のとき,乙円の直径は 1 寸 2 分である。
算額の問題を解くということだけならば,ここまでで終了。
図を描くために,その他のパラメータを求める。
次に,r2 が決定されたという条件下で,弦の交点座標を求める。
r2 = R*r1/(R + r1)
eq2 = x^2 + y^2 - R^2 |> expand
eq3 = (R + y)*y2 - x2*x |> expand;
res2 = solve([eq2, eq3], (x, y))[1]
(y2*(R + (R*x2^2 - R*y2^2)/(x2^2 + y2^2))/x2, (R*x2^2 - R*y2^2)/(x2^2 + y2^2))
最後に,r2, x, y が決定されたという条件下で,乙円の中心座標 x2, y2 を求める。
eq4 = x2^2 + y2^2 - (R - r2)^2 |> expand
eq5 = dist2(0, R, x, -y, x2, y2, r2) |> expand;
eq14 = eq4(x => res2[1], y => res2[2])
ans_x2 = solve(eq14, x2)[1]
ans_x2 |> println
sqrt(R^4 - R^2*y2^2 - 2*R*r1*y2^2 - r1^2*y2^2)/(R + r1)
eq15 = eq5(x => res2[1], y => res2[2]) |> simplify |> numerator
eq25 = eq15(x2 => ans_x2);
ans_y2 = solve(eq25, y2)[2]
ans_y2 |> println
R^2*(R - r1)/(R^2 + 2*R*r1 + r1^2)
y2 の値も,外円の半径と甲円の半径から求めることができる。
その後は,逆順にたどると全てのパラメータの値が求まる。
(R, r1) = (1.5, 1)
r2 = R*r1/(R + r1)
y2 = R^2*(R - r1)/(R^2 + 2*R*r1 + r1^2)
x2 = sqrt(R^4 - R^2*y2^2 - 2*R*r1*y2^2 - r1^2*y2^2)/(R + r1)
x = y2*(R + (R*x2^2 - R*y2^2)/(x2^2 + y2^2))/x2
y = (R*x2^2 - R*y2^2)/(x2^2 + y2^2)
(x2, y2, x, y) |> println
(0.8818163074019442, 0.18, 0.5878775382679626, 1.3800000000000001)
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r1) = (1.5, 1)
r2 = R*r1/(R + r1)
y2 = R^2*(R - r1)/(R^2 + 2*R*r1 + r1^2)
x2 = sqrt(R^4 - R^2*y2^2 - 2*R*r1*y2^2 - r1^2*y2^2)/(R + r1)
x = y2*(R + (R*x2^2 - R*y2^2)/(x2^2 + y2^2))/x2
y = (R*x2^2 - R*y2^2)/(x2^2 + y2^2)
r3 = R - r1
@printf("r2 = %g; x2 = %g; y2 = %g; r3 = %g; x = %g; y = %g\n", r2, x2, y2, r3, x, y)
plot()
circle(0, 0, R)
circle(0, R - r1, r1, :blue)
circle(0, r3 - R, r3, :green)
circle(x2, y2, r2, :magenta)
segment(0, R, x, -y, :orange)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:black, lw=1)
vline!([0], color=:gray80, lw=0.5)
point(0, R, " R", :red, :left, :bottom, delta=delta/2)
point(x, -y, "(x,y)", :red, :left, delta=-delta/2)
point(0, R - r1, "甲円:r1 \n(0,R-r1) ", :blue, :right, :vcenter)
point(x2, y2, "乙円:r2 \n(x2,y2) ", :magenta, :right, :vcenter)
point(0, r3 - R, "丙円:r3,(0,r3-R) ", :green, :center, delta=-delta/2)
point(0, r3 - R)
end
end;
※コメント投稿者のブログIDはブログ作成者のみに通知されます