神壁算法 東都本郷春木町一丁目 佐野富次郎安行 天明七年
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
キーワード：円5個，楕円
#Julia #SymPy #算額 #和算 #数学
楕円の中に大円 2 個，中円 2 個，小円 1 個が入っている。
楕円の長径，短径がそれぞれ 2168 寸，1084 寸のとき，大円の直径はいかほどか。
楕円の長半径と短半径を \(a,\ b\)
大円の半径と中心座標を \(r_1,\ (x_1,\ 0);\ x_1 = r_1 + r_3\)
中円の半径と中心座標を \(r_2,\ (0,\ b - r_2)\)
小円の半径と中心座標を \(r_3,\ (0,\ 0)\)
楕円と大円の共通接点の座標を \( (x_0,\ y_0)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, b::positive, x0::positive, y0::positive,
r1::positive, x1::positive,r2::positive, r3::positive
(a, b) = (2168, 1084) .// 2
x1 = r1 + r3
eq1 = (x0 - x1)^2 + y0^2 - r1^2
eq2 = -b^2*x0/(a^2*y0) - (x1 - x0)/y0
eq3 = x0^2/a^2 + y0^2/b^2 - 1
eq4 = r3 + 2r2 - b
eq5 = x1^2 + (b - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (x0, y0, r1, r2, r3))[1]
(-1626/7 + 2710*sqrt(57)/21, sqrt(6315926/147 + 734410*sqrt(57)/49), 271*sqrt(57)/14 + 4065/14, 3523/7 - 271*sqrt(57)/7, -3252/7 + 542*sqrt(57)/7)
大円の直径は (271√57 + 4065)/7 = 873.000733136910 寸である。
ちなみに，中円，小円の直径はそれぞれ 421.998533726179 寸，240.002932547642 寸である。
res[3].evalf()*2 |> println
res[4].evalf()*2 |> println
res[5].evalf()*2 |> println
873.000733136910
421.998533726179
240.002932547642
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b) = (2168, 1084) .// 2
(x0, y0, r1, r2, r3) = (-1626/7 + 2710*sqrt(57)/21, sqrt(6315926/147 + 734410*sqrt(57)/49), 271*sqrt(57)/14 + 4065/14, 3523/7 - 271*sqrt(57)/7, -3252/7 + 542*sqrt(57)/7)
x1 = r1 + r3
plot()
ellipse(0, 0, a, b, color=:orange)
circle(x1, 0, r1)
circle(-x1, 0, r1)
circle(0, 0, r3, :green)
circle(0, b - r2, r2, :blue)
circle(0, r2 - b, r2, :blue)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
point(a, 0, " a", :black, :left, :bottom, delta=delta)
point(0, b, " b", :black, :left, :bottom, delta=delta)
point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta)
point(x1, 0, "大円:r1,(x1,y1)", :red, :center, delta=-delta)
point(0, b - r2, " 中円:r2,(0,b-r2)", :black, :left, :vcenter)
point(0, r3, " r3", :black, :left, :bottom, delta=delta)
point(0, r3, " 小円", :black, :left, :top, delta=-delta)
end
end;
以下のアイコンをクリックして応援してください
