福島県耶麻郡猪苗代町中小松西浜 猪苗代小平潟天満宮 明治14年(1881)
http://www.wasan.jp/fukusima/kohiragata3.html
街角の数学 Street Wasan ～落書き帳「○△□」～ 330．小平潟天満宮（その３） 明治14年(1881)
http://streetwasan.web.fc2.com/math17.11.1.html
キーワード：円8個，楕円，等脚台形
#Julia #SymPy #算額 #和算 #数学
等脚台形の中に大円が内接しており，大円に 3 個の楕円が内接している。楕円は隣の楕円と大円に一点で接する。台形の頂点近辺にある 6 個の小円の直径が 1 寸のとき，真ん中にある内円の直径はいかほどか。
「問」には「楕円は隣の楕円と大円に一点で接する」としか書いていないが，「大円と一点で接する『短径が最も長い楕円（曲率楕円）』」である。短径が短い（すなわち内円の直径が大きい）場合には解が定まらない。
『算法助術の公式86』で，直径 \(d\) の円に 1 点で内接する曲率楕円の長径 \(a\) と短径 \(b\) の関係が述べられている \(d = a^2/b\)
大円の中心を原点に置く。
大円の半径と中心座標を \(R,\ (0,\ 0)\)
台形の右上と右下の頂点の座標を \( (x_b,\ R),\ (x_a,\ -R)\)
上方の楕円の長径と短径と中心座標を \(a,\ b,\ (0,\ r_1 + b)\)
内円の半径と中心座標を \(r_1,\ (0,\ 0)\)
小円の半径と中心座標を \(r_2,\ (x_1,\ R - r_2),\ (x_2,\ y_2),\ (x_3,\ r_2)\)
とおき以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R, a, b, r1, x0, y0, xa::positive, xb::positive, r2::positive, x1::positive, x2::positive, y2::negative, x3::positive
@syms d
eq1 = 2b + r1 - R
eq2 = R - a^2/b # b = (R - r1)/2
eq3 = x0^2/a^2 + (y0 - r1 - b)^2/b^2 - 1
eq4 = -b^2*x0/(a^2*(y0 - r1 - b)) - 1/sqrt(Sym(3))
eq5 = x0/y0 - sqrt(Sym(3));
eq6 = dist(xa, -R, xb, R, x1, R - r2) - r2^2
eq7 = dist(xa, -R, xb, R, x2, y2) - r2^2
eq8 = dist(0, 0, xa, - R, x2, y2) - r2^2
eq9 = x3^2 + (r2 - R)^2 - (R + r2)^2;
eq10 = x2^2 + y2^2 - (R + r2)^2
eq11 = xa*xb - R^2;
eq12 = x1^2 + (R - r2)^2 - (R + r2)^2;
function H(u)
(R, r1, a, b, x0, y0, xa, xb, x1, x2, y2, x3) = u
return [
-R + 2*b + r1, # eq1
R - a^2/b, # eq2
-1 + (-b - r1 + y0)^2/b^2 + x0^2/a^2, # eq3
-sqrt(3)/3 - b^2*x0/(a^2*(-b - r1 + y0)), # eq4
x0/y0 - sqrt(3), # eq5
-r2^2 + (2*R - 2*R*(2*R*(2*R - r2) + (x1 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2) - r2)^2 + (x1 - xa - (-xa + xb)*(2*R*(2*R - r2) + (x1 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2))^2, # eq6
-r2^2 + (R - 2*R*(2*R*(R + y2) + (x2 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2) + y2)^2 + (x2 - xa - (-xa + xb)*(2*R*(R + y2) + (x2 - xa)*(-xa + xb))/(4*R^2 + (-xa + xb)^2))^2, # eq7
-r2^2 + (x2 - xa*(-R*y2 + x2*xa)/(R^2 + xa^2))^2 + (R*(-R*y2 + x2*xa)/(R^2 + xa^2) + y2)^2, # eq8
x3^2 + (-R + r2)^2 - (R + r2)^2, # eq9
x2^2 + y2^2 - (R + r2)^2, # eq10
-R^2 + xa*xb, # eq11
x1^2 + (R - r2)^2 - (R + r2)^2, # eq12
]
end;
r2 = 1/2
iniv = BigFloat[7.1, 1.1, 4.6, 3.1, 3.2, 1.8, 8, 6.1, 5.3, 6.6, -4.5, 3]
res = nls(H, ini=iniv)
([3.5, 0.5, 2.29128784747792, 1.5, 1.5155444566227676, 0.875, 3.968626966596886, 3.086709862908689, 2.6457513110645907, 3.307189138830738, -2.25, 2.6457513110645907], true)
小円の直径が 1 寸のとき，内円の直径も 1 寸である。
なお，大円の直径は小円の直径の 7 倍，楕円の短径は小円の直径の 3 倍である。
その他のパラメータは以下のとおりである。
\(R = 3.5; \ r_1 = 0.5;\ a = 2.29129;\ b = 1.5;\ x_0 = 1.51554;\ y_0 = 0.875\)
\(x_a = 3.96863;\ x_b = 3.08671;\ x_1 = 2.64575;\ x_2 = 3.30719\)
\(y_2 = -2.25;\ x_3 = 2.64575\)
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 1/2
(R, r1, a, b, x0, y0, xa, xb, x1, x2, y2, x3) = res[1]
@printf("R = %g; r1 = %g; a = %g; b = %g; x0 = %g; y0 = %g; xa = %g; xb = %g; x1 = %g; x2 = %g; y2 = %g; x3 = %g\n", R, r1, a, b, x0, y0, xa, xb, x1, x2, y2, x3)
plot()
circle(0, 0, R)
circle(0, 0, r1, :green)
l = r1 + b
ellipse(0, r1 + b, a, b, color=:blue)
x = l*cos(pi/6)
y = -l*sin(pi/6)
ellipse(x, y, a, b, color=:blue, φ= 240)
ellipse(-x, y, a, b, color=:blue, φ= 120)
circle(x3, r2 - R, r2)
circle(-x3, r2 - R, r2)
segment(-xa, -R, xa, -R)
segment(-xb, R, xb, R)
segment(xa, -R, xb, R)
segment(-xa, -R, -xb, R)
circle(x1, R - r2, r2)
circle(-x1, R - r2, r2)
circle(x2, y2, r2)
circle(-x2, y2, r2)
if more
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)
point(x0, y0, "(x0,y0)", :blue, :left, delta=-delta/2)
point(0, R, "R", :blue, :center, :bottom, delta=delta/2)
point(xa, -R, "(xa,-R)", :black, :right, :bottom, delta=delta/2)
point(xb, R, "(xb,R)", :black, :right, :bottom, delta=delta/2)
point(0, r1, " r1", :blue, :center, :bottom, delta=delta/2)
point(0, R - b, "R-b", :blue, :center, :bottom, delta=delta/2)
point(a, R - b, "(a,R-b) ", :blue, :right, :vcenter)
point(x1, R - r2, "(x1,R-r2)", :red, :right, :bottom, delta=delta)
point(x2, y2, "小円:r2,(x2,y2)", :red, :right, :vcenter)
point(x3, r2 - R, " (x3,r2-R)", :red, :left, :vcenter)
point(0, 0, "内円:r1", :green, :center, :vcenter, mark=false)
end
end;
以下のアイコンをクリックして応援してください
