九十六 岩手県大船渡市立根町 気仙安養寺稲荷堂 文政5年(1822)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html
キーワード：円7個，楕円，正方形
#Julia #SymPy #算額 #和算 #数学
正方形の中に楕円，元円，利円，貞円，享円を容れる。貞円の直径が 1 寸のとき，元円の直径はいかほどか。
注：享は正しくは亨である。⇒ 元亨利貞
注：利円，享円の条件については何も述べられていない。実際，制限はあるもののどんな値も取りうる。
正方形の一辺の長さを \(2a\)
楕円の長半径，短半径，中心座標を \(a,\ b,\ (0,\ b);\ b = a - r_1\)
元円の半径と中心座標を \(r_1,\ (0,\ 2a - r_1),\ (r_1,\ a),\ (r_1,\ a - 2r_1)\)
貞円の半径と中心座標を \(r_2,\ (a - r_3,\ r_3)\)
楕円と元円の接点座標を \( (x_1,\ y_1)\)
楕円と貞円の接点座標を \( (x_2,\ y_2)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, b::positive,
r1::positive, x1::positive, y1::positive,
r2::positive, x2::positive, y2::positive
b = a - r1
eq1 = x2^2/a^2 + (y2 - a + r1)^2/b^2 - 1
eq2 = (x2 - a + r2)^2 + (y2 - r2)^2 - r2^2
eq3 = -b^2*x2/(a^2*(y2 - a + r1)) - (a - r2 - x2)/(y2 - r2)
eq4 = x1^2/a^2 + (y1 - a + r1)^2/b^2 - 1
eq5 = (x1 - r1)^2 + (y1 - a)^2 - r1^2
eq6 = -b^2*x1/(a^2*(y1 - a + r1)) - (r1 - x1)/(y1 - a);
function H(u)
(a, r1, x1, y1, x2, y2) = u
return [
-1 + (-a + r1 + y2)^2/(a - r1)^2 + x2^2/a^2, # eq1
-r2^2 + (-r2 + y2)^2 + (-a + r2 + x2)^2, # eq2
-(a - r2 - x2)/(-r2 + y2) - x2*(a - r1)^2/(a^2*(-a + r1 + y2)), # eq3
-1 + (-a + r1 + y1)^2/(a - r1)^2 + x1^2/a^2, # eq4
-r1^2 + (-a + y1)^2 + (-r1 + x1)^2, # eq5
-(r1 - x1)/(-a + y1) - x1*(a - r1)^2/(a^2*(-a + r1 + y1)), # eq6
]
end;
r2 = 1/2
iniv = BigFloat[2.1, 0.63, 0.83, 2.6, 1.5, 0.52]
res = nls(H, ini=iniv)
([3.5742972321311055, 1.1362300699229417, 1.470601088978576, 4.6602138004791485, 2.755899335723288, 0.8855162507120506], true)
元円の直径は 2*1.1362300699229417 = 2.2724601398458835 である。「答」では 2.243 寸有奇となっており，かなり異なった結論になった。
ちなみに，楕円の長半径，短半径は 3.5743，\(b\) = 2.43807 とこれまた，あまりきれいな数ではない。
2*1.1362300699229417
2.2724601398458835
その他のパラメータは以下のとおりである。
\(r_2 = 0.5;\ a = 3.5743;\ b = 2.43807;\ r_1 = 1.13623;\ r_2 = 0.5\)
\(x_1 = 1.4706;\ y_1 = 4.66021;\ x_2 = 2.7559;\ y_2 = 0.885516\)
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 1/2
(a, r1, x1, y1, x2, y2) = res[1]
b = a - r1
@printf("貞円の直径が %g のとき，元円の直径は %g である。\n", 2r2, 2r1)
@printf("r2 = %g; a = %g; b = %g; r1 = %g; r2 = %g; x1 = %g; y1 = %g; x2 = %g; y2 = %g\n",
r2, a, b, r1, r2, x1, y1, x2, y2)
plot([a, a, -a, -a, a], [0, 2a, 2a, 0, 0], color=:green, lw=0.5)
ellipse(0, a - r1, a, b, color=:blue)
circle(0, 2a - r1, r1)
circle2(r1, a, r1)
circle2(r1, a - 2r1, r1)
circle2(a - r2, r2, r2, :green)
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(x2, y2, "(x2,y2) ", :green, :right, :bottom, delta=2delta, deltax=4delta)
point(x1, y1, "(x1,y1)", :green, :left, :bottom, delta=delta/2)
point(r1, a, "元円:r1,(r1, a)", :red, :center, delta=-delta/2)
point(r1, a - 2r1, "元円:r1,(r1, a-2r1)", :red, :center, :bottom, delta=delta/2)
point(0, 2a - r1, "元円:r1,(0,2a-r1)", :red, :center, delta=-delta/2)
point(a - r2, r2, "貞円:r3,(a-r3,r3)", :green, :center, delta=-delta/2)
point(0, 2a - 2r1, "2a-2r1", :blue, :center, delta=-delta/2)
point(0, a - r1, "a-r1", :blue, :center, delta=-delta/2)
end
end;
以下のアイコンをクリックして応援してください
