算額（その1095）改訂版 - 算額あれこれ

63 岩手県花泉町花泉字東鹿野天満社文政13年(1830)

安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円9個，楕円，正方形
#Julia, #SymPy, #算額, #和算

正方形の中に，大円 1 個，等円 8 個，楕円 1 個を容れる。正方形の一辺の長さが 1 寸のとき，楕円の短径はいかほどか。

正方形の一辺の長さを \(s\)
大円の半径と中心座標を \(r_1, (0, 0)\)
等円の半径と中心座標を \(r_2, (r_2, 0), (0, r_1 - r_2), (x - r_2, s/2 - r_2)\)
楕円の長半径と短半径，中心座標を \(a, b, (0, 0); a = r_1, b = r_1 - 2r_2\)
とおき，以下の連立方程式を解く。
しかし，（なぜか）SymPy で連立方程式を解くと不適切解しか得られないので，まず eq1 を解いて \(r_1\) を求め，しかる後に eq2 から \(r_2\) を求める。

include("julia-source.txt"); # julia-source.txt ソース

@syms s::positive, a::positive, b::positive, r1::positive, r2::positive
a = r1
b = r1 - 2r2
# 算法助術-公式84
eq1 = (2r2)^2 - 4(a^2 - b^2)*(b^2 - r2^2)/b^2
# 大円と右上の等円が外接する
eq2 = 2(s/2 - r2)^2 - (r1 + r2)^2;
# res = solve([eq1, eq2], (r1, r2))

# r1: 大円の半径 -- この時点ではまだ r2 に対する倍数
ans_r1 = solve(eq1, r1)[1]

\(\displaystyle r_{2} \left(\frac{19}{48 \sqrt[3]{\frac{\sqrt{87}}{36} + \frac{23}{64}}} + \sqrt[3]{\frac{\sqrt{87}}{36} + \frac{23}{64}} + \frac{7}{4}\right)\)

ans_r1.evalf() |> println

3.06659283332063*r2

# r2：等円の半径 ## 長い式になるので式の表示は割愛
ans_r2 = solve(eq2(r1 => ans_r1), r2)[1];

ans_r2.evalf() |> println

0.129015099263883*s

# r1: 大円の半径　-- s に対する倍数
ans_r1(r2 => ans_r2).evalf() |> println

0.395636778792771*s

# b: 楕円の短半径 ## 長い式になるので式の表示は割愛
b = (ans_r1 - 2ans_r2)(r2 => ans_r2);
b.evalf() |> println

0.137606580265006*s

楕円の短半径は，正方形の一辺の長さの 0.1376065802650062 倍である。
正方形の一辺の長さが 1 寸のとき，楕円の短径は 2*0.1376065802650062 = 0.27521316053001 寸である。

描画関数プログラムのソースを見る

function draw(s, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, a, b) = (s/2) .* [0.791273557585543, 0.258030198527765, 0.791273557585543, 0.275213160530012]
plot([1, 1, -1, -1, 1].*(s/2), [-1, 1, 1, -1, -1].*(s/2), color=:black, lw=0.5)
circle(0, 0, r1, :blue)
circle4(s/2 - r2, s/2 - r2, r2)
circle2(r2, 0, r2)
circle22(0, r1 - r2, r2)
ellipse(0, 0, a, b, color=:orange)
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(a, 0, "a ", :orange, :right, :bottom, delta=delta/2)
point(0, b, "b", :orange, :left, delta=-delta/2, deltax=delta/2)
point(0, r1, "r1", :blue, :center, :bottom, delta=delta/2)
point(r2, 0, "r2", :red, :center, delta=-delta/2)
point(0, r1 - r2, "r1-r2", :red, :center, delta=-delta/2)
point(s/2 - r2, s/2 - r2, "(s/2-r2,s/2-r2)", :red, :center, delta=-delta/2)
point(s/2, s/2, "(s/2,s/2)", :black, :right, :bottom, delta=delta/2)
end
end;