埼玉県加須市 秋葉山本坊(秋葉宮) 文化5年(1808)
和算の館
http://www.wasan.jp/saitama/akibahonbo1.html
百四十三 群馬県榛名町榛名山 榛名神社 明治33年(1900)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
早坂平四郎: 算額の一考察
苫小牧工業専門学校紀要
https://www.tomakomai-ct.ac.jp/uploads/files/wp/kiyou/old/kiyou5-8.pdf
キーワード：円5個，外円，斜線
#Julia #SymPy #算額 #和算 #数学
円内に大円1個，中円2個，小円2個が入っているこれらは互いに内接・外接している他2,弦を共通接線としている。中円，小円の半径がわかっているとき，大円の半径を求めよ。
外円の半径と中心座標を \(r_0,\ (0,\ 0)\)
大円の半径と中心座標を \(r_1,\ (0,\ y_1)\)
中円の半径と中心座標を \(r_2,\ (r_2,\ y_2)\)
小円の半径と中心座標を \(r_3,\ (r_3,\ y_3)\)
弦の端点座標を \( (xa,\ ya),\ (xb,\ yb)\)
ただし，ya^2 \(= r_0^2 - xa^2,\ yb^2 = r_0^2 - xb^2\)
solve() では解ききらないので nlsolve() で数値解を求める。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r0::positive, r1::positive, y1::positive,
r2::positive, y2::negative, r3::positive, y3::positive,
xa::positive, ya::negative, xb::positive, yb::positive;
(r2, r3) = (15, 5)
eq1 = r3^2 + y3^2 - (r0 - r3)^2
eq2 = r2^2 + y2^2 - (r0 - r2)^2
eq3 = r3^2 + (y3 - y1)^2 - (r1 + r3)^2
eq4 = r2^2 + (y1 - y2)^2 - (r1 + r2)^2
eq5 = distance(xa, ya, xb, yb, 0, y1) - r1^2
eq6 = distance(xa, ya, xb, yb, r2, y2) - r2^2
eq7 = distance(xa, ya, xb, yb, r3, y3) - r3^2;
eq8 = xa^2 + ya^2 - r0^2
eq9 = xb^2 + yb^2 - r0^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (r0, r1, y1, y2, y3, xa, xb))
function H(u)
(r0, r1, y1, y2, y3, xa, ya, xb, yb) = u
return [
y3^2 - (r0 - 5)^2 + 25, # eq1
y2^2 - (r0 - 15)^2 + 225, # eq2
-(r1 + 5)^2 + (-y1 + y3)^2 + 25, # eq3
-(r1 + 15)^2 + (y1 - y2)^2 + 225, # eq4
-r1^2 + (y1 - (xa^2*yb - xa*xb*ya - xa*xb*yb + xb^2*ya + y1*ya^2 - 2*y1*ya*yb + y1*yb^2)/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2))^2 + (xa*y1*ya - xa*y1*yb - xa*ya*yb + xa*yb^2 - xb*y1*ya + xb*y1*yb + xb*ya^2 - xb*ya*yb)^2/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2)^2, # eq5
(y2 - (ya*(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) - (ya - yb)*(xa^2 - xa*xb - 15*xa + 15*xb - y2*ya + y2*yb + ya^2 - ya*yb))/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2))^2 + (-(xa*(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) - (xa - xb)*(xa^2 - xa*xb - 15*xa + 15*xb - y2*ya + y2*yb + ya^2 - ya*yb))/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) + 15)^2 - 225, # eq6
(y3 - (ya*(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) - (ya - yb)*(xa^2 - xa*xb - 5*xa + 5*xb - y3*ya + y3*yb + ya^2 - ya*yb))/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2))^2 + (-(xa*(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) - (xa - xb)*(xa^2 - xa*xb - 5*xa + 5*xb - y3*ya + y3*yb + ya^2 - ya*yb))/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) + 5)^2 - 25, # eq7
-r0^2 + xa^2 + ya^2, # eq8
-r0^2 + xb^2 + yb^2, # eq9
]
end;
(r2, r3) = (15, 5)
iniv = [big"40.0", 18, 10, -20.4, 34.1, 32, -24.4, 8.2, 39.1]
res = nls(H, ini=iniv);
println(res);
([38.763883748662835, 17.99038105676658, 10.951523807752837, -18.43155367352307, 33.39161340506353, 32.33934584554869, -21.37300618915924, 8.622548477793336, 37.79272867931013], true)
\(r_0 = 38.7639;\ r_1 = 17.9904;\ y_1 = 10.9515;\ y_2 = -18.4316\)
\(y_3 = 33.3916;\ xa = 32.3393;\ ya = -21.373;\ xb = 8.62255;\ yb = 37.7927\)
中円，小円の半径をそれぞれ 15, 5 としたとき，大円の半径は 17.99038105676658 である。
Float64(res[1][2])
17.99038105676658
術によれば，中円，小円の半径を a, b とすると，大円の「直径」は \( \displaystyle \frac{(a + b + 6\sqrt{a b})(a + b + 2\sqrt{a b})}{4a b}\) であるとしている。しかしこれは上述の大円の半径を2倍したものではない。何らかの不適切な式であると思われる。
(a, b) = (15, 5)
(a + b + 6√(a*b))*(a + b + 2√(a*b))/(4a*b)
8.952135486850342
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, r3) = (15, 5)
(r0, r1, y1, y2, y3, xa, ya, xb, yb) = res[1]
@printf("r0 = %g; r1 = %g; y1 = %g; y2 = %g; y3 = %g; xa = %g; ya = %g; xb = %g; yb = %g\n", r0, r1, y1, y2, y3, xa, ya, xb, yb)
plot()
circle(0, 0, r0, :black)
circle(0, y1, r1)
circle(r2, y2, r2, :orange)
circle(-r2, y2, r2, :orange)
circle(r3, y3, r3, :green)
circle(-r3, y3, r3, :green)
segment(xa, ya, xb, yb, :blue)
segment(-xa, ya, -xb, yb, :blue)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(r0, 0, " r0", :black, :left, :bottom, delta=delta)
point(0, y1, " y1 大円:r1", :red, :left, :vcenter)
point(r2, y2, " 中円:r2,(r2,y2)", :orange, :center, :top, delta=-delta)
point(r3, y3, " 小円:r3,(r3,y3)", :green, :left, :vcenter)
point(xa, ya, "(xa,ya)", :blue, :left, :top, delta=-delta)
point(xb, yb, "(xb,yb)", :blue, :left, :bottom, delta=delta)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
end
end;
以下のアイコンをクリックして応援してください
