神壁算法 〇二 上州白雲山 寛政9年(1797)
藤田貞資門人 上州板鼻驛 小野十五郎榮重
藤田嘉言(1807)：続神壁算法
http://www.wasan.jp/jinpeki/jinpeki.html
長野県長野市松代町西条 西楽寺 弘化5年(1848)
中村信弥「改訂増補 長野県の算額」県内の算額(P.176)
http://www.wasan.jp/zoho/zoho.html
キーワード：外円，累円
#Julia #SymPy #算額 #和算 #数学
甲円は外円に内接し，乙円，丙円...は隣の円と甲円に外接し，外円に内接している。外円，甲円，乙円の直径がそれぞれ 168 寸，88 寸，3 寸のとき，乙円から始まる累円のうち最も大きい円の直径を求めよ。
甲円，乙円のパラメータを初期値として，丙円，丁円...のパラメータは順次計算できる。
デカルトの円定理で直径だけを求めることはできるが，図を描くためには累円の中心座標も必要なので，逐次的に数値解を求める。
外円の半径と中心座標を \(R,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (0,\ r_1 - R)\)
乙円の半径と中心座標を \(r_2,\ (x_2,\ y_2)\)
として，以下の連立方程式で \( (x_2,\ y_2)\) を求める。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive, y2::negativez
(R, r1, r2) = (168, 88, 3) .// 2
eq1 = x2^2 + (r1 - R - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
res1 = solve([eq1, eq2], (x2, y2))[1]
(231/10, -396/5)
現在の円のパラメータ（初期状態では \(r_2,\ x_2,\ y_2\)）から，次円のパラメータを求める関数は以下のようになる。
解は 2 通り求まる。一方は累円が大きくなる方向，もう一方は累円が小さくなる方向である。
累円が一番大きくなると解は振動する。
次に，最大の累円のパラメータを再設定して今度は累円が小さくなる方向で累円のパラメータを求める。
@syms R::positive, r1::positive, r2::positive, x2::positive, y2
@syms r3, x3, y3
function nextcircle(R, r1, r2, x2, y2)
eq3 = x3^2 + (r1 - R - y3)^2 - (r1 + r3)^2
eq4 = x3^2 + y3^2 - (R - r3)^2
eq5 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
res2 = solve([eq3, eq4, eq5], (r3, x3, y3))
end;
実際に計算すると，累円が最大になるのは甲円から数えて 10 番目の円である。
円の半径は 77/2，中心座標は (-231/10, 196/5) である。
n = 2 (3//2, 231//10, -396//5) ## 乙円
n = 3 (231/118, 15477/590, -22932/295)
n = 4 (77/29, 4389/145, -10948/145)
n = 5 (231/61, 10857/305, -21924/305)
n = 6 (231/40, 8547/200, -1638/25)
n = 7 (77/8, 2079/40, -266/5)
n = 8 (231/13, 3927/65, -1764/65)
n = 9 (33, 231/5, 108/5)
n = 10 (77/2, -231/10, 196/5)
n = 11 (231/10, -3003/50, -252/25)
n = 12 (231/19, -5313/95, -4284/95)
n = 13 (7, -231/5, -308/5)
n = 14 (231/52, -9933/260, -4536/65)
n = 15 (231/76, -12243/380, -7056/95)
n = 16 (11/5, -693/25, -1924/25)
n = 17 (231/139, -16863/695, -54684/695)
n = 18 (231/178, -19173/890, -35532/445)
n = 19 (77/74, -7161/370, -14924/185)
n = 20 (231/271, -23793/1355, -110124/1355)
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r1) = (168, 88) .// 2
plot()
circle(0, 0, R)
circlef(0, r1 - R, r1, :green)
point(0, r1 - R, " 1", :white, :left, :vcenter, mark=false)
(r2, x2, y2) = (3//2, 231//10, -396//5)
for n = 2:10
println("n = ", n, " ", (r2, x2, y2))
circlef(x2, y2, r2, n)
point(x2, y2, string(n), :white, :center, :vcenter, mark=false)
(r2, x2, y2) = nextcircle(168//2, 88//2, r2, x2, y2)[2]
end
(r2, x2, y2) = (77//2, -231//10, 196//5)
for n = 11:20
(r2, x2, y2) = nextcircle(168//2, 88//2, r2, x2, y2)[1]
point(x2, y2, string(n), :white, :center, :vcenter, mark=false)
@printf("n = %d; r2 = %.15g; x2 = %.15g; y2 = %.15g\n", n, r2, x2, y2)
circlef(x2, y2, r2, n)
end
end;
以下のアイコンをクリックして応援してください
