算額（その2132） - 算額あれこれ

長野県駒ヶ根市大字赤穂光前寺明治37年(1904)

「算額」第四集全国調査，香川県算額研究会.（香川県立図書館蔵）

中村信弥「改訂増補長野県の算額」県内の算額(P.264)
http://www.wasan.jp/zoho/zoho.html

キーワード：3次元，球11個
#Julia, #SymPy, #算額, #和算，#数学

大球 1 個，中球 2 個，小球 8 個がある。中球は上下に接し，小球は互いに外接し上下の中球にも接している大球は下の中球と小球が内接する。大球の直径が 8 寸のとき，小球の直径はいかほどか。

「算額（その1800）」と図形的には同じで，小球の個数が特定され，大球が与えられて小球の直径を求める点が違う。

大球の半径と中心座標を \(r_1,\ (0,\ 0,\ r_1 - r_2)\)
中球の半径と中心座標を \(r_2,\ (0,\ 0,\ r_2)\)
小球の半径と中心座標を \(r_3,\ (x_3,\ 0,\ r_2)\)
とおき，以下の連立方程式の解を求める。

include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r1::positive, r2::positive, r3::positive, x3::positive
eq1 = x3^2 + r2^2 - (r2 + r3)^2
eq2 = x3^2 + (2r2 - r1 )^2 - (r1 - r3)^2
eq3 = x3*sind(Sym(22.5)) - r3
res = solve([eq1, eq2, eq3], (r3, x3, r2))[1];

# r3: 小球の半径
ans_r3 = res[1]
@show(ans_r3)

ans_r3 = r1*(-4 + 3*sqrt(2))

\(r_{1} \left(-4 + 3 \sqrt{2}\right)\)

大球の直径が 8 寸のとき，小球の直径は 1.94112549695428 寸である。

ans_r3(r1 => 8/2).evalf() * 2

\(1.94112549695428\)

「答」は小球の直径は 2.314 あまりとなっている。
「術」は以下のとおりであり，大球の直径が 8 寸のとき，小球の直径は 2.314575050761981 となる。
上で求めた 1.94112549695428 と異なる。

大球の直径を 2 倍にすると小球の径も 2 倍になるべきであるが，術ではそのようにならない。この術は間違っている。

function 術(大球径)
甲 = √8 + 3
( (2甲 + 大球径)*大球径)/(甲^2*2)
end;

術(8)

2.314575050761981

術(16)

6.513134198986965

術(16)/術(8)

2.813965439073914

上述の ans_r3 ではちゃんと 2 倍になる。

(ans_r3(r1 => 16)/ans_r3(r1 => 8)).evalf()

\(2.0\)

# 小円の中心の x 座標
@syms d
ans_x3 = res[2]
ans_x3 = apart(ans_x3/r1, d) |> simplify |> x -> r1*x
@show(ans_x3)

ans_x3 = 2*r1*(-1 + sqrt(2))*sqrt(2 - sqrt(2))

\(2 r_{1} \left(-1 + \sqrt{2}\right) \sqrt{2 - \sqrt{2}}\)

ans_x3(r1 => 8/2).evalf()

\(2.53620268449772\)

# r2: 中球の半径
ans_r2 = res[3]
@show(ans_r2)

ans_r2 = sqrt(2)*r1/2

\(\displaystyle \frac{\sqrt{2} r_{1}}{2}\)

ans_r2(r1 => 8/2).evalf()

\(2.82842712474619\)

その他のパラメータも用いて図を描くと，小球の半径 = 1.94112549695428 で矛盾はない。
「答」，「術」ともに間違えている。

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

function draw(r1, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, showaxis=true, label="", fontfamily="IPAMincho")
r3 = r1*(-4 + 3√2)
x3 = 2*r1*(√2 - 1)*sqrt(2 - √2)
r2 = r1/√2
p1 = plot(xlabel="x-axis", ylabel="y-axis")
rotate(x3, 0, r3, :blue, angle=45)
circle(0, 0, r1)
circle(0, 0, r2, :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(x3, 0, "x3", :blue, :right, deltax=-delta/2)
ylims!(-5, 5)
end
p2 = plot(xlabel="x-axis", ylabel="z-axis")
circle(0, r1 - r2, r1)
circle(0, 0, r2, :magenta)
circle(0, 2r2, r2, :magenta)
circle2(x3, r2, r3, :blue)
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(0, 0, "中球:r2,(0,0,r2)", :magenta, :center, delta=-delta)
point(0, 2r2, "中球:r2,(0,0,2r2)", :magenta, :center, :bottom, delta=delta)
point(0, r1 - r2, "大球:r1,(0,0,r1-r2)", :red, :left, delta=-delta, deltax=-17delta)
point(x3, r2, "小球", :blue, :center, delta=-delta)
ylims!(-3.5, 8.5)

end
plot(p1, p2)
end;
draw(8/2, true)