30 岩手県一関市赤荻字宿 正慶山観音寺 天保2年(1831)
安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円6個，半円，四分円，直角三角形，正方形
#Julia #SymPy #算額 #和算 #数学
正方形の中に四分円1個，半円2個，元円1個，亨円2個，利円1個，貞円2個，直角三角形1個，小さな正方形1個を容れる。元円か亨円の直径の和が与えられたとき，丁円の直径と小さな正方形の一辺の長さの和はいかほどか。
注：「只謂知元円径歟亨円径歟」の「歟」は「元円か亨円の直径を知って」ということであろう。実は，元円の直径と亨円の直径は等しい（点対称であるため）。
大正方形の一辺の長さを \(a\)
四分円の半径と中心座標を \(r_6,\ (0,\ 0);\ r_6 = a\)
半円の半径と中心座標を \(r_5,\ (0,\ r_5),\ (a,\ r_5);\ r_5 = a/2\)
元円の半径と中心座標を \(r_1,\ (a - r_1,\ y_1)\)
亨円の半径と中心座標を \(r_2,\ (r_2,\ y_2),\ (r_2,\ y_2 + 2r_2)\)
利円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
貞円の半径と中心座標を \(r_4,\ (x_4,\ y_4)\)
小正方形の一辺の長さを \(2b\)
とおき，以下の連立方程式を解く。
大正方形の一辺の長さ \(a\) が決まれば，残りのパラメータは \(a\) の倍数で決まる。
正方形の一辺の長さを \(a = 1\) としても，一般性は失われない。
まず，元円の半径と中心の \(y\) 座標を求める。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, r1::positive, r2::positive, r3::positive,
r4::positive, r5::positive, r6::positive,
x::positive, y1::positive, y2::positive
a = 1
r5 = a//2
r6 = a
eq1 = (a - r1)^2 + y1^2 - (a + r1)^2 # 元円と四分円が外接
eq2 = r1^2 + (y1 - r5)^2 - (r5 - r1)^2 # 元円が右の半円に内接
res1 = solve([eq1, eq2], (r1, y1))[1];
# r1
res1[1]
\( \displaystyle \frac{4}{25}\)
# y1
res1[2]
\( \displaystyle \frac{4}{5}\)
次に，下にある亨円の半径と中心の \(y\) 座標を求める。
y = r5 - sqrt(r5^2 - x^2) # (x, y) が左の半円の周上にある
eq3 = r2^2 + (r5 - y2)^2 - (r5 - r2)^2 # 下の亨円が左の半円に内接
eq4 = dist(0, r6, x, y, r2, y + r2) - r2^2 # 上の亨円が直角三角形に内接する
eq5 = (y2 + 2r2) - (y + r2); # 長さの関係式
res2 = solve([eq3, eq5], (y2, x))[2]; # 2 of 2
# y2
res2[1]
\( \displaystyle \frac{1}{2} - \frac{\sqrt{1 - 4 r_{2}}}{2}\)
# x
res2[2]
\( \displaystyle \frac{\sqrt{1 - \left(- 2 r_{2} + \sqrt{1 - 4 r_{2}}\right)^{2}}}{2}\)
eq4 の \(y_2,\ x\) に代入し，方程式 eq14 を得る。
eq14 = eq4(y2 => res2[1], x => res2[2]) |> simplify |> numerator |> expand
\( \displaystyle 128 r_{2}^{5} - 256 r_{2}^{4} \sqrt{1 - 4 r_{2}} + 64 r_{2}^{4} \sqrt{- 4 r_{2}^{2} + 4 r_{2} \sqrt{1 - 4 r_{2}} + 4 r_{2}} - 32 r_{2}^{4} \sqrt{4 r_{2}^{2} - 4 r_{2} \sqrt{1 - 4 r_{2}} - 4 r_{2} + 1} - 864 r_{2}^{4} - 64 r_{2}^{3} \sqrt{1 - 4 r_{2}} \sqrt{- 4 r_{2}^{2} + 4 r_{2} \sqrt{1 - 4 r_{2}} + 4 r_{2}} + 64 r_{2}^{3} \sqrt{1 - 4 r_{2}} \sqrt{4 r_{2}^{2} - 4 r_{2} \sqrt{1 - 4 r_{2}} - 4 r_{2} + 1} + 448 r_{2}^{3} \sqrt{1 - 4 r_{2}} - 32 r_{2}^{3} \sqrt{- 4 r_{2}^{2} + 4 r_{2} \sqrt{1 - 4 r_{2}} + 4 r_{2}} \sqrt{4 r_{2}^{2} - 4 r_{2} \sqrt{1 - 4 r_{2}} - 4 r_{2} + 1} - 160 r_{2}^{3} \sqrt{- 4 r_{2}^{2} + 4 r_{2} \sqrt{1 - 4 r_{2}} + 4 r_{2}} + 256 r_{2}^{3} \sqrt{4 r_{2}^{2} - 4 r_{2} \sqrt{1 - 4 r_{2}} - 4 r_{2} + 1} + 896 r_{2}^{3} + 32 r_{2}^{2} \sqrt{1 - 4 r_{2}} \sqrt{- 4 r_{2}^{2} + 4 r_{2} \sqrt{1 - 4 r_{2}} + 4 r_{2}} \sqrt{4 r_{2}^{2} - 4 r_{2} \sqrt{1 - 4 r_{2}} - 4 r_{2} + 1} + 96 r_{2}^{2} \sqrt{1 - 4 r_{2}} \sqrt{- 4 r_{2}^{2} + 4 r_{2} \sqrt{1 - 4 r_{2}} + 4 r_{2}} - 128 r_{2}^{2} \sqrt{1 - 4 r_{2}} \sqrt{4 r_{2}^{2} - 4 r_{2} \sqrt{1 - 4 r_{2}} - 4 r_{2} + 1} - 256 r_{2}^{2} \sqrt{1 - 4 r_{2}} + 64 r_{2}^{2} \sqrt{- 4 r_{2}^{2} + 4 r_{2} \sqrt{1 - 4 r_{2}} + 4 r_{2}} \sqrt{4 r_{2}^{2} - 4 r_{2} \sqrt{1 - 4 r_{2}} - 4 r_{2} + 1} + 128 r_{2}^{2} \sqrt{- 4 r_{2}^{2} + 4 r_{2} \sqrt{1 - 4 r_{2}} + 4 r_{2}} - 160 r_{2}^{2} \sqrt{4 r_{2}^{2} - 4 r_{2} \sqrt{1 - 4 r_{2}} - 4 r_{2} + 1} - 288 r_{2}^{2} + 32 r_{2} \sqrt{1 - 4 r_{2}} \sqrt{4 r_{2}^{2} - 4 r_{2} \sqrt{1 - 4 r_{2}} - 4 r_{2} + 1} + 32 r_{2} \sqrt{1 - 4 r_{2}} - 32 r_{2} \sqrt{- 4 r_{2}^{2} + 4 r_{2} \sqrt{1 - 4 r_{2}} + 4 r_{2}} \sqrt{4 r_{2}^{2} - 4 r_{2} \sqrt{1 - 4 r_{2}} - 4 r_{2} + 1} - 32 r_{2} \sqrt{- 4 r_{2}^{2} + 4 r_{2} \sqrt{1 - 4 r_{2}} + 4 r_{2}} + 32 r_{2} \sqrt{4 r_{2}^{2} - 4 r_{2} \sqrt{1 - 4 r_{2}} - 4 r_{2} + 1} + 32 r_{2}\)
eq14 は下図のような見かけで，0.15 前後に解がある。
pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(eq14, xlims=(0,0.25))
hline!([0])
Roots.find_zero で解を求めると, \(r_2 = 0.16\) となる。
using Roots
ans_r2 = find_zero(eq14, 0.1, 0.2)
0.15999999999999992
直角三角形の斜辺の長さは 0.8 である。
sqrt(x^2 + (r6 - y)^2)(x => res2[2])(r2 => 0.16) |> println
0.800000000000000
また，\(x\) の数値解は 0.48 である。
ans_x = res2[2](r2 => 0.16) |> println
0.480000000000000
さらに，利円，貞円のパラメータ \(r_3,\ x_3,\ y_3,\ r_4,\ x_4,\ y_4\) を求める（数値解）。
@syms r3, x3, y3, r4, x4, y4
eq5 = x3^2 + (y3 - r5)^2 - (r5 - r3)^2 # 利円が半円に内接
eq6 = x4^2 + (y4 - r5)^2 - (r5 - r4)^2 # 貞円が半円に内接
eq7 = (x4 - x3)^2 + (y4 - y3)^2 - (r4 + r3)^2 # 利円と貞円が外接
eq8 = 0.8^2 - 4*(2r5)*(2r4) # 算法助術-公式29
eq9 = dist2(0, r6, x, y, x3, y3, r3) |> numerator # 利円が直角三角形の斜辺に接する
eq10 = dist2(0, r6, x, y, x4, y4, r4) |> numerator # 貞円が直角三角形の斜辺に接する
function H(u)
(r3, x3, y3, r4, x4, y4) = u
return [
x3^2 - (1/2 - r3)^2 + (y3 - 1/2)^2, # eq5
x4^2 - (1/2 - r4)^2 + (y4 - 1/2)^2, # eq6
-(r3 + r4)^2 + (-x3 + x4)^2 + (-y3 + y4)^2, # eq7
0.64 - 8.0*r4, # eq8
-2*r3^2*x^4 + 2*r3^2*x^2*sqrt(1 - 4*x^2) + 4*r3^2*x^2 - r3^2*sqrt(1 - 4*x^2) - r3^2 + 4*x^5*x3*y3 - 4*x^5*x3 + x^4*x3^2*sqrt(1 - 4*x^2) + 5*x^4*x3^2 - x^4*y3^2*sqrt(1 - 4*x^2) - 3*x^4*y3^2 + 2*x^4*y3*sqrt(1 - 4*x^2) + 6*x^4*y3 - x^4*sqrt(1 - 4*x^2) - 3*x^4 - 4*x^3*x3*y3*sqrt(1 - 4*x^2) - 8*x^3*x3*y3 + 4*x^3*x3*sqrt(1 - 4*x^2) + 8*x^3*x3 - 3*x^2*x3^2*sqrt(1 - 4*x^2) - 5*x^2*x3^2 + x^2*y3^2*sqrt(1 - 4*x^2) + x^2*y3^2 - 2*x^2*y3*sqrt(1 - 4*x^2) - 2*x^2*y3 + x^2*sqrt(1 - 4*x^2) + x^2 + 2*x*x3*y3*sqrt(1 - 4*x^2) + 2*x*x3*y3 - 2*x*x3*sqrt(1 - 4*x^2) - 2*x*x3 + x3^2*sqrt(1 - 4*x^2) + x3^2, # eq9
-2*r4^2*x^4 + 2*r4^2*x^2*sqrt(1 - 4*x^2) + 4*r4^2*x^2 - r4^2*sqrt(1 - 4*x^2) - r4^2 + 4*x^5*x4*y4 - 4*x^5*x4 + x^4*x4^2*sqrt(1 - 4*x^2) + 5*x^4*x4^2 - x^4*y4^2*sqrt(1 - 4*x^2) - 3*x^4*y4^2 + 2*x^4*y4*sqrt(1 - 4*x^2) + 6*x^4*y4 - x^4*sqrt(1 - 4*x^2) - 3*x^4 - 4*x^3*x4*y4*sqrt(1 - 4*x^2) - 8*x^3*x4*y4 + 4*x^3*x4*sqrt(1 - 4*x^2) + 8*x^3*x4 - 3*x^2*x4^2*sqrt(1 - 4*x^2) - 5*x^2*x4^2 + x^2*y4^2*sqrt(1 - 4*x^2) + x^2*y4^2 - 2*x^2*y4*sqrt(1 - 4*x^2) - 2*x^2*y4 + x^2*sqrt(1 - 4*x^2) + x^2 + 2*x*x4*y4*sqrt(1 - 4*x^2) + 2*x*x4*y4 - 2*x*x4*sqrt(1 - 4*x^2) - 2*x*x4 + x4^2*sqrt(1 - 4*x^2) + x4^2, # eq10
]
end;
x = 0.48
iniv = BigFloat[0.11, 0.3, 0.7, 0.079, 0.4, 0.6]
res = nls(H, ini=iniv)
([0.1, 0.3200000000000002, 0.7399999999999998, 0.08, 0.41133126291998995, 0.5848916494400132], true)
最後に小さい正方形の一辺の長さは以下の方程式を解いて求める。
@syms b
eq11 = (r5 - b)^2 + (r5 - 2b)^2 - r5^2
ans_b = solve(eq11, b)[1]
ans_b |> println
1/10
まとめると，大正方形の一辺の長さが \(a = 1\) のとき，
\(r_1 = 0.16;\ r_2 = 0.16;\ r_3 = 0.1;\ r_4 = 0.08;\ r_5 = 0.5;\ r_6 = 1;\ b = 0.1\)
元円の直径は 0.32
亨円の直径は 0.32
利円の直径は 0.20
貞円の直径は 0.16
正方形の一辺の長さは 0.2
元円と右下の亨円は点対称なので，計算しなくても等しいことはわかる。
元円と亨円の直径はともに 0.32
貞円の直径と正方形の一辺の長さの和は 0.36
0.36/0.32 = 9/8 = 1.125
「貞円の直径と正方形の一辺の長さの和」は「元円（または亨円）の直径」の 9/8 = 1.125 倍である。
描画関数プログラムのソースを見る
function draw(a, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r6 = a
r5 = a/2
(r1, y1) = (4*a/25, 4*a/5)
(r2, y2, x) = [1.6, 2.0, 4.8]./10
y = r5 - sqrt(r5^2 - x^2)
(r3, x3, y3, r4, x4, y4) = [0.1, 0.3200000000000002, 0.7399999999999998, 0.08, 0.41133126291998995, 0.5848916494400132]
b = 0.1
@printf("r1 = %g; r2 = %g; r3 = %g; r4 = %g; r5 = %g; r6 = %g; b = %g\n", r1, r2, r3, r4, r5, r6, b)
θ1 = atand(y3 - r5, x3)
θ2 = atand(y4 - r5, x4)
length = sqrt(x4^2 + (y4- r5)^2)
x42 = length*cosd(2θ1 - θ2)
y42 = length*sind(2θ1 - θ2) + r5
plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
plot!(r5 .+ [-b, b, b, -b, -b], [0, 0, 2b, 2b, 0], color=:orange, lw=0.5)
circle(0, 0, r6, beginangle=0, endangle=90)
circle(a, r5, r5, :blue, beginangle=90, endangle=270)
circle(0, r5, r5, :blue, beginangle=-90, endangle=90)
circle(a - r1, y1, r1, :magenta)
circle(r2, y2, r2, :purple)
circle(r2, y2 + 2r2, r2, :purple)
segment(0, y, x, y)
segment(0, a, x, y)
circle(x3, y3, r3, :green)
circle(x4, y4, r4, :tomato)
circle(x42, y42, r4, :tomato)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
vline!([0], color=:gray80, lw=0.5)
hline!([0], color=:gray80, lw=0.5)
point(a - r1, y1, "元円:r1,(a-r1,y1)", :magenta, :center, delta=-delta)
point(r2, y2, "亨円:r2,(r2,y2)", :purple, :center, delta=-delta)
point(r2, y2 + 2r2, "亨円:r2\n(r2,y2+2r2)", :purple, :center, delta=-delta)
point(x3, y3, "利円:r3\n(x3,y3)", :green, :center, delta=-delta/2)
point(x4, y4, " 貞円:r4,(x4,y4)", :black, :left, :vcenter)
point(x42, y42, " 貞円:r4,(x42,y42)", :black, :left, :vcenter)
#point(r2, y + r2, "(r2,y+r2)", :magenta, :left, :bottom, delta=2delta)
point(x, y, "(x,y) ", :black, :right, delta=-delta)
point(0, r6, "a=r6 ", :red, :right, :vcenter)
point(0, r5, "2/a=r5 ", :blue, :right, :vcenter)
point(r5 - b, 2b, "(r5-b,2b)", :orange, :left, :bottom, delta=delta/2, deltax=1.75delta)
point(r5 - b, 0, " r5-b", :orange, :left, :bottom, delta=delta/2)
point(r5 + b, 0, " r5+b", :orange, :left, :bottom, delta=delta/2)
xlims!(-12delta, r6 + 3delta)
println(4r2, " ", 2r4 + 2b)
end
end;
draw(1, true)
以下のアイコンをクリックして応援してください
