十七 群馬県高崎市八幡町 八幡宮 文化7年(1810)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円5個，菱形
#Julia #SymPy #算額 #和算 #数学
菱形の中に，大円 2 個，中円 2 個，小円 1 個を容れる。菱形の対角線の長い方が 12 寸，短い方が 5 寸のとき，小円の直径はいかほどか。
菱形の対角線を \(2a,\ 2b;\ a > b\)
大円の半径と中心座標を \(r_1,\ (r_3 + r_1,\ 0)\)
中円の半径と中心座標を \(r_2,\ (0,\ r_3 + r_2)\)
小円の半径と中心座標を \(r_3,\ (0,\ 0)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, b::positive,
r1::positive, r2::positive, r3::positive;
sinθ = b/sqrt(a^2 + b^2)
cosθ = a/sqrt(a^2 + b^2)
eq1 = (r3 + r1)^2 + (r3 + r2)^2 - (r1 + r2)^2 |> expand
eq2 = (a - r3 - r1)*sinθ - r1
eq3 = (b - r3 - r2)*cosθ - r2;
SymPy では能力的に一度に解けないので，逐次解いていく。
まず，eq1 を解いて \(r_1\) を求める。
ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println
r3*(r2 + r3)/(r2 - r3)
eq2 の \(r_1\) に，上で求めた \(r_1\) を代入し，\(r_2\) を求める。
eq12 = eq2(r1 => ans_r1)
ans_r2 = solve(eq12, r2)[1]
ans_r2 |> println
r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))
eq3 に，上で求めた \(r_1,\ r_2\) を代入し，\(r_3\) を求める。
eq13 = eq3(r1 => ans_r1, r2 => ans_r2) |> simplify |> numerator
ans_r3 = solve(eq13, r3)[2] |> factor
ans_r3 |> println
-a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2
解として得られた \(r_3\) を 代入して \(r_2,\ r_1\) を確定する。
a = 12/2
b = 5/2
r3 = -a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2
r2 = r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))
r1 = r3*(r2 + r3)/(r2 - r3)
(r1, r2, r3) |> println
(1.5296184351192643, 0.9631806558860881, 0.49337363357064845)
菱形の対角線の長い方が 12 寸，短い方が 5 寸のとき，小円の直径は 2*0.49337363357064845 = 0.9867472671412969 寸である。
「答」では「小円の径一寸三分二厘半微弱」としているが，正しいとは思えない。
すべてのパラメータは以下のとおりである。
\( a = 6;\ b = 2.5;\ r_1 = 1.52962;\ r_2 = 0.963181;\ r_3 = 0.493374\)
描画関数プログラムのソースを見る
function draw(a, b, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r3 = -a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2
r2 = r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))
r1 = r3*(r2 + r3)/(r2 - r3)
@printf("菱形の対角線の長い方が %g，短い方が %g のとき，小円の直径は %g である。\n", 2a, 2b, 2r3)
@printf("a = %g; b = %g; r1 = %g; r2 = %g; r3 = %g\n", a, b, r1, r2, r3)
plot([a , 0, -a, 0, a], [0, b, 0, -b, 0], color=:green, lw=0.5)
circle2(r3 + r1, 0, r1)
circle22(0, r3 + r2, r2, :blue)
circle(0, 0, r3, :magenta)
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", :green, :left, :vcenter)
point(0, b, "b", :green, :center, :bottom, delta=delta)
point(r3 + r1, 0, "大円:r1\n(r3+r1,0)", :red, :center, delta=-delta)
point(0, r3 + r2, "中円:r2\n(0,r3+r2)", :blue, :center, :bottom, delta=delta)
point(0, 0, "小円:r3,(0,0)", :magenta, :right, :vcenter, deltax=-7delta)
end
end;
draw(12/2, 5/2, true)
以下のアイコンをクリックして応援してください
