算額(その799)
藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
等脚台形の中に斜線を引き,区画された領域に甲円と乙円を入れる。乙円の直径が 15 寸,上頭(上底)の長さが 12 寸,斜線の長さが 20 寸のとき,甲円の直径はいかほどか。
等脚台形の高さを h,下頭(下底),上頭の長さをそれぞれ 2a, 2b
甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (0, h - r2)
とおいて以下の連立方程式を解く。
なお,斜の長さについては算法助術の公式 39 を援用した。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, h::positive,
r1::negative, r2::positive, 斜::positive, d
eq1 = h^2 + (a - b)^2 - (斜 + a + b)^2
eq2 = dist(a, 0, b, h, 0, r1) - r1^2
eq2 = numerator(apart(eq2, d))
eq3 = dist(a, 0, b, h, 0, h - r2) - r2^2
eq3 = numerator(apart(eq3, d))
res = solve([eq1, eq2, eq3], (r1, a, h));
res[2]
(-b*斜/(2*r2) + r2 + r2*斜/(2*b), -斜/2 + r2^2/b + r2^2*斜/(2*b^2), r2*(2*b + 斜)/b)
2 組の解が得られるが,2 番目のものが適解である。
b = 12//2
斜 = Sym(20)
r2 = 15//2
(-b*斜/(2*r2) + r2 + r2*斜/(2*b), -斜/2 + r2^2/b + r2^2*斜/(2*b^2), r2*(2*b + 斜)/b)
(12, 15, 40)
甲円の半径の一般解は r2*(1 + 斜/2b) - b*斜/2r2 である
r2*(1 + 斜/2b) - b*斜/2r2 |> println
12
図を描くために斜と等脚辺の交点座標 (x1, y1), (x2, y2) を求める。
一般解は求めにくいようなので,問の条件下で交点座標を求める。
@syms x1::positive, y1::positive, x2::negative, y2::positive, d
(a, h, r1) = (15, 40, 12)
eq4 = dist(x1, y1, x2, y2, 0, r1) - r1^2
eq4 = numerator(apart(eq4, d))
eq5 = dist(x1, y1, x2, y2, 0, h - r2) - r2^2
eq5 = numerator(apart(eq5, d))
eq6 = y1*(b - a) - h*(x1 - a)
eq7 = y2*(a - b) - h*(x2 + a)
res2 = solve([eq4, eq5, eq6, eq7], (x1, y1, x2, y2));
res2[2]
(9*sqrt(10)/41 + 390/41, 1000/41 - 40*sqrt(10)/41, -390/41 + 9*sqrt(10)/41, 40*sqrt(10)/41 + 1000/41)
2 組の解が得られるが,2 番目のものが適解である。
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
b = 12//2
斜 = Sym(20)
r2 = 15//2
(a, h, r1) = (15, 40, 12)
@printf("甲円の直径 = %g; 下頭 = %g; 高さ = %g\n", 2r1, 2a, h)
plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:black, lw=0.5)
circle(0, r1, r1)
circle(0, h - r2, r2, :blue)
if more == true
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)
(x1, y1, x2, y2) = (9*sqrt(10)/41 + 390/41, 1000/41 - 40*sqrt(10)/41, -390/41 + 9*sqrt(10)/41, 40*sqrt(10)/41 + 1000/41)
segment(x1, y1, x2, y2, :green)
point(b, h, "(b,h)", :black, :right, :bottom, delta=delta/2)
point(a, 0, "(a,0) ", :black, :right, :bottom, delta=delta/2)
point(x1, y1, "(x1,y1)", :black, :right, delta=-delta)
point(x2, y2, "(x2,y2) ", :black, :left, delta=-2delta)
point(0, r1, " 甲円:r1,(0,r1)", :red, :left, :vcenter)
point(0, h - r2, " 乙円:r2\n (0,h-r2)", :blue, :left, :vcenter)
end
end;