結論：DWARF 3 の圧勝
撮影は 2026/01/06
DWARF 3
暗部と半暗部がくっきりと写っている
かなり拡大してもそれなりに見える
Seestar S50
ピント調整してもぼんやりとしている
拡大しようという気になれない
算額（その1997）
六十二 群馬県高崎市石原町 清水寺 天保10年(1839)
キーワード：楕円，長方形，最小値
#Julia #SymPy #算額 #和算 #数学
楕円の中に長方形を容れる。長方形の長辺の長さが与えられたとき，楕円の面積が最小になるときの楕円の長径を求める術を述べよ。
楕円の長径，短径，面積を \(2a,\ 2b,\ S; S = \pi a b\)
長方形の長辺，短辺を \(2x,\ 2y\)
とおき，以下の連立方程式を解いて \(S,\ b\) を求める。
以下の計算に用いるのは，長半径 \(a\)，短半径 \(b\)，長辺の半分 \(x\)，短辺の半分 \(y\) なので混乱しないように注意すること。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a, b, x::positive, y, S
eq1 = x^2/a^2 + y^2/b^2 ⩵ 1 # 楕円の方程式
eq2 = S ⩵ PI*a*b; # 楕円の面積
res = solve( (eq1, eq2), (S, b))[2] # 2 of 2
(pi*a^2*y*sqrt(1/( (a - x)*(a + x))), a*y*sqrt(1/( (a - x)*(a + x))))
# S: 楕円の面積
f = res[1]
\(\displaystyle \pi a^{2} y \sqrt{\frac{1}{\left(a - x\right) \left(a + x\right)}}\)
たとえば，\(x = 2,\ y = 1\) のときの \(S,\ a\) のグラフを描いてみると，以下のようになる。
using LaTeXStrings
using Plots
pyplot(size=(300, 150), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(f(x => 2, y => 1), xlims=(2.4, 3.7), xlabel=L"a", ylims=(12, 14.5), ylabel=L"S=f(a)")
\(x = 2,\ y = 1\) のとき，楕円の面積は \(a =2.8\) 近辺で最小値を取ることがわかる。
正確にいえば，まず \(f(a)\) の導関数 \(g(a)\) を求め，\(g(a_0) = 0\) となる \(a_0\) を求める。
そして，\(a = a_0\) のとき \(f(a_0)\) は最小値を取る。
# 導関数
g = diff(f, a) |> simplify
\(\displaystyle \frac{\pi a y \left(- a^{2} + 2 x^{2}\right) \sqrt{\frac{1}{a^{2} - x^{2}}}}{- a^{2} + x^{2}}\)
# a0: 楕円の面積が最小になるときの楕円の長半径
a0 = solve(g, a)[3]
\(\sqrt{2} x\)
# a = √2x のときの，楕円の面積の最小値
f(a => √2x)
\(2.0 \pi x y\)
\(a_0\) が長方形の一辺の長さの半分 \(x\) の\(\sqrt{2}\) 倍のとき，楕円の面積は最小値 \(2\pi x y\) を取る。
たとえば，\(x = 2,\ y = 1\) のとき，\(a_0 = \sqrt{2}x = 2.82842712474619\) のとき，楕円の面積の最小値は \(2\pi x y = 4\pi = 12.5663706143592\) になる。
なお，そのときの楕円の短半径は \(b_0 = a_0 y/\sqrt{a_0^2 - x^2} = 1.41421356237309\) である。楕円の面積は \(\pi a b = \pi\cdot 2.82842712474619 \cdot 1.41421356237309 = 4\pi = 12.566370614359128\) となり，当然ながら一致する。
a0(x => 2).evalf()
\(2.82842712474619\)
f(a => √2x)(x => 2, y => 1).evalf()
\(12.5663706143592\)
# b0: 楕円の面積が最小になるときの楕円の短半径
res[2](a => √2x, y => 1) |> println
1.41421356237309
using LaTeXStrings
描画関数プログラムのソースを見る
function draw(x, y, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x, y) = (x, y) ./ 2
a = √2x
b = a*y/sqrt(a^2 - x^2)
@printf("長方形の長辺 = %g, 短辺 = %g のとき，楕円の長径 = %g, 短径 = %g のとき面積が最小値 %g になる。\n", 2x, 2y, 2a, 2b, 2pi*x*y)
plot()
rect(-x, -y, x, y, :blue)
ellipse(0, 0, a, b, color=:red, lw=0.5)
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(x, y, L"(x,y)", :red, :left, :bottom, delta=delta)
point(√2x, 0, L"\sqrt{2}x", :red, :left, :bottom, delta=delta)
point(0, a*y/sqrt(a^2 - x^2), L"\frac{a y}{\sqrt{a^2 - x^2}}", :red, :left, :bottom, delta=delta)
end
end;
draw(4, 2, true)
長方形の長辺 = 4, 短辺 = 2 のとき，楕円の長径 = 5.65685, 短径 = 2.82843 のとき面積が最小値 12.5664 になる。
以下のアイコンをクリックして応援してください
算額（その1996）
五十 群馬県高崎市八幡町 八幡宮 文政11年(1828)
五十一 群馬県安中市柳瀬 稲荷神社 文政11年(1828)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円8個，外円
#Julia #SymPy #算額 #和算 #数学
外円の中に外円に内接し，互いに外接する 4 個の円があり，4 円に囲まれて隙間に，隣同士また外の 4 円に外接する 3 個の円がある。外の直径が 20 寸，外側の 4 円のうち左右の円の直径がそれぞれ 5 寸，4 寸のとき，内部にある 3 円のうち中央にあるものの直径はいかほどか。
この算額問題は，算額（その1995）を一段階複雑にした（全体を内接する外円を加えた）ものである。
「群馬の算額」では，
「高崎市八幡町 八幡宮 文政十一年戊子十月 関流岩井重遠門人 上毛碓氷郡新井邑 岩井藤右衛門藤原貴重，上原多蔵，橘秀重」
「安中市柳瀬 稲荷神社 文政十一年戊子十月 上毛高崎藩 小野栄重門人 富田又五郎源篤忠」
が所載されている。前者の出典は算法雑爼，後者は第2問目で，「第1問は算法雑爼，第2,第3問は中曽根家蔵の資料によった」という注がついている。
両者は円の名前が違うだけで，与える条件が 20寸，5寸，4寸，答が 1 寸で，術も同じである。
算法雑爼には岩井らの名前で同じ算額題がある。これも，円の名前はまた別であるが，与える条件も同じ（20寸，5寸，4寸）である。しかし，答えは 1 寸ではなく，分数で 20/21 寸とされている。
https://kokusho.nijl.ac.jp/biblio/100245066/22?ln=ja
外円の中心を座標原点に置き，等円の中心が \(y\) 軸上になるように配置する
外円の半径と中心座標を \(R,\ (0,\ 0)\)
天円の半径と中心座標を \(r_1,\ (x_1,\ y_1)\)
地円の半径と中心座標を \(r_2,\ (x_2,\ y_2)\)
人円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
東円の半径と中心座標を \(r_4,\ (0,\ y_4);\ y_4 = R - r_4\)
西円の半径と中心座標を \(r_5,\ (x_5,\ y_5)\)
南円の半径と中心座標を \(r_6,\ (x_6,\ y_6)\)
北円の半径と中心座標を \(r_7,\ (x_7,\ y_7)\)
とおく。
\(R = 20/2,\ r_1 = 5/2,\ r_2 = 4/2\) が与えられ，
\(x_1,\ y_1,\ x_2,\ y_2,\ r_3,\ x_3,\ y_3,\ r_4,\ y_4,\ r_5,\ x_5,\ y_5,\ r_6,\ x_6,\ y_6,\ r_7,\ x_7,\ y_7\) の 18 変数を求める。
外円を含め 8 個の円から外接・内接する 2 円の半径の間に成り立つピタゴラスの定理を表す方程式が以下の 18 本（eq1〜eq18）なので，18 元連立方程式を解けばよいことになる。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R, r1, x1, y1, r2, x2, y2, r3, x3, y3,
r4, x4, y4, r5, x5, y5, r6, x6, y6, r7, x7, y7
eq1 = x1^2 + y1^2 - (R - r1)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = y4 - (R - r4)
eq4 = x5^2 + y5^2 - (R - r5)^2
eq5 = (x1 - x5)^2 + (y1 - y5)^2 - (r1 + r5)^2
eq6 = (x5 - x6)^2 + (y5 - y6)^2 - (r5 + r6)^2
eq7 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq8 = (x5 - x7)^2 + (y5 - y7)^2 - (r5 + r7)^2
eq9 = (x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2
eq10 = x1^2 + (y1 - y4)^2 - (r1 + r4)^2
eq11 = x6^2 + (y4 - y6)^2 - (y4 + y6)^2
eq12 = x3^2 + (y3 - y4)^2 - (r3 + r4)^2
eq13 = x7^2 + (y4 - y7)^2 - (r4 + r7)^2
eq14 = x2^2 + (y2 - y4)^2 - (r2 + r4)^2
eq15 = (x1 - x6)^2 + (y1 - y6) - (r1 + r6)^2
eq16 = (x3 - x6)^2 + (y3 - y6)^2 - (r3 + r6)^2
eq17 = (x3 - x7)^2 + (y3 - y7)^2 - (r3 + r7)^2
eq18 = (x2 - x7)^2 + (y2 - y7)^2 - (r2 + r7)^2;
SymPy の能力的に solve() で解析解を求めるのは無理なので，NLSsolve() で数値解を求める。
初期値は GeoGebra で略図を描いて座標値などを読み取って用いる。
Julia では，多重ディスパッチ機能のお陰で，nls(), nlsolve() の引数に BigFloat を使用することにより，内部で自動的に BigFloat で長精度計算が行われ，計算誤差の影響を極力抑えることが可能になる。
通常の使用では結果を Float64（倍精度実数）に丸めて返すようにしているが，今回は BigFloat のまま返すようにした。
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve( (vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
v = r.zero[1]
else
r = nlsolve( (vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
v = r.zero
end
return v, r.f_converged
# return Float64.(v), r.f_converged
end;
function H(u)
(x1, y1, x2, y2, r3, x3, y3, r4, y4, r5, x5, y5, r6, x6, y6, r7, x7, y7) = u
return [
x1^2 + y1^2 - (R - r1)^2, # eq1
x2^2 + y2^2 - (R - r2)^2, # eq2
y4 - (R - r4), # eq3
x5^2 + y5^2 - (R - r5)^2, # eq4
(x1 - x5)^2 + (y1 - y5)^2 - (r1 + r5)^2, # eq5
(x5 - x6)^2 + (y5 - y6)^2 - (r5 + r6)^2, # eq6
(x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2, # eq7
(x5 - x7)^2 + (y5 - y7)^2 - (r5 + r7)^2, # eq8
(x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2, # eq9
x1^2 + (y1 - y4)^2 - (r1 + r4)^2, # eq10
x6^2 + (y4 - y6)^2 - (r4 + r6)^2, # eq11
x3^2 + (y3 - y4)^2 - (r3 + r4)^2, # eq12
x7^2 + (y4 - y7)^2 - (r4 + r7)^2, # eq13
x2^2 + (y2 - y4)^2 - (r2 + r4)^2, # eq14
(x1 - x6)^2 + (y1 - y6) - (r1 + r6)^2, # eq15
(x3 - x6)^2 + (y3 - y6)^2 - (r3 + r6)^2, # eq16
(x3 - x7)^2 + (y3 - y7)^2 - (r3 + r7)^2, # eq17
(x2 - x7)^2 + (y2 - y7)^2 - (r2 + r7)^2 # eq18
]
end;
(R, r1, r2) = (20, 5, 4)./2
iniv = BigFloat[3.85, 6.49, # x1, y1
-3.76, 7.06, # x2, y2
0.55, -0.26, 5.82, # r3, x3, y3
1.84, 8.16, # r4, y4
7.57, -1.20, -2.13, # 5
0.6, 0.86, 5.77, # 6
0.61, -1.37, 6.04 # 7
]
res = nls(H, ini=iniv);
数値解を求めるプロセスで，収束誤差はほとんどすべてが \(10^{-48}\) 以下である。
(R, r1, r2) = (20, 5, 4)./2
(x1, y1, x2, y2, r3, x3, y3, r4, y4, r5, x5, y5, r6, x6, y6, r7, x7, y7) = res[1];
err = Float64[x1^2 + y1^2 - (R - r1)^2, # eq1
x2^2 + y2^2 - (R - r2)^2, # eq2
y4 - (R - r4), # eq3
x5^2 + y5^2 - (R - r5)^2, # eq4
(x1 - x5)^2 + (y1 - y5)^2 - (r1 + r5)^2, # eq5
(x5 - x6)^2 + (y5 - y6)^2 - (r5 + r6)^2, # eq6
(x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2, # eq7
(x5 - x7)^2 + (y5 - y7)^2 - (r5 + r7)^2, # eq8
(x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2, # eq9
x1^2 + (y1 - y4)^2 - (r1 + r4)^2, # eq10
x6^2 + (y4 - y6)^2 - (r4 + r6)^2, # eq11
x3^2 + (y3 - y4)^2 - (r3 + r4)^2, # eq12
x7^2 + (y4 - y7)^2 - (r4 + r7)^2, # eq13
x2^2 + (y2 - y4)^2 - (r2 + r4)^2, # eq14
(x1 - x6)^2 + (y1 - y6) - (r1 + r6)^2, # eq15
(x3 - x6)^2 + (y3 - y6)^2 - (r3 + r6)^2, # eq16
(x3 - x7)^2 + (y3 - y7)^2 - (r3 + r7)^2, # eq17
(x2 - x7)^2 + (y2 - y7)^2 - (r2 + r7)^2 # eq18
]
println(err)
[3.2420074929264436e-49, 1.3827635414815871e-49, 1.3817869688151111e-76, 2.0609691246985557e-50, 1.8139310823774985e-49, 8.066445768180866e-49, 1.1921631520035702e-48, 3.8013880507087696e-49, 2.2118860295102422e-49, 2.054238379708111e-49, 1.0515566542938712e-48, 1.0080238769804959e-48, 2.305837382746626e-49, 4.325740633857128e-50, 4.1581213270227974e-50, 3.5266186897459675e-48, 2.3706486600615317e-50, 5.762909830943632e-51]
長精度の結果を表示する。
names = ("x1", "y1", "x2", "y2", "r3", "x3", "y3", "r4", "y4", "r5", "x5", "y5", "r6", "x6", "y6", "r7", "x7", "y7")
for (s, x) in zip(names, res[1])
println(s, " = ", x)
end
x1 = 3.9931122639050275990681304130169810535135524160526295770802464444491295058537
y1 = 6.348626185865038409569537109664589380369938389152814991944504128854924760564188
x2 = -3.726816517969065155816050678907394993167551308865867366418687120620827240380214
y2 = 7.078900948692030727655629687731671504295950711322251993396821096967697475773503
r3 = 0.4793196182266118254204998254094961948183291828789941139244845589472333276665717
x3 = -0.2437605578704490094389560995573206639087383474369953455585263463989130088209249
y3 = 5.78988238923057952377261607518815808989767523942241194519851443913406798643635
r4 = 1.871734427014128400223955397974279154052027993086429144970295838662820285477141
y4 = 8.128265572985871599776044602025720845947972006913570855029704161337179714522997
r5 = 7.499999380223074495604088525695951163168559786215925457373750192763773134186097
x5 = -1.333470226742621839462912773345721310497530662323208852522431401702108970067703
y5 = -2.114677340228525819596416959143938476049803774783763905503314643067434022188707
r6 = 0.6771303413884256148653332818715354338627062958725435074628840676459152561960812
x6 = 0.9119351353477939512342204702263355272434614285612389589501720468792751203818402
y6 = 5.748121477238304041624410108913693121481204581555038603319292170270580479370239
r7 = 0.6264023391342447139809467498908996331599696045356664081550706537598574273406719
x7 = -1.327000291833266999322080431895226158570523074071858113270293313467688612663804
y7 = 6.01172180356916623766874010517771385505362853160558648181889752552837471195104
人円の直径は約 0.9586 寸であり，群馬の算額での「答：人円径は 1 寸」ではない。
# 人円の直径
2 * res[1][5]
0.9586392364532236508409996508189923896366583657579882278489691178944666553331433
rationalize(2*res[1][5], tol=1e-6), 1043/1088
(1043//1088, 0.9586397058823529)
rationalize(2*res[1][5], tol=1e-3), 23/24
(23//24, 0.9583333333333334)
rationalize(0.958)
479//500
算法雑爼では人円の直径を 20/21 = 0.9523809523809523 としているが，正確に 20/21 ではないのではないかと思われる。
最初に示した図は，得られた解に基づくものである。
模式図的に描かれている算額での図とはかなり異なることが明らかになった。
問答では「外円径 = 20 寸，天円径 5 寸で，地円径は 4 寸，人円径は 1 寸」天円径は外円径の 1/4 = 0.25 倍，地円径は外円径の 4/20 = 0.2，人円径は外円径の 1/20 = 0.05 である。
群馬の算額に描かれている図では（無謀であるが）ノギスで測定すると，外円径 = 37.2mm，天円径 = 11.7mm, 地円径 = 14.3mm，人円径 2.4mm なので，11.7/37.2 = 0.314516129032258，14.3/37.2 = 0.3844086021505376，2.4/37.2 = 0.06451612903225806 である。まるっきり参考にできない。
描画関数プログラムのソースを見る
function draw(R, r1, r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x1, y1, x2, y2, r3, x3, y3, r4, y4, r5, x5, y5, r6, x6, y6, r7, x7, y7) =
[
3.85, 6.49, # x1, y1
-3.76, 7.06, # x2, y2
0.55, -0.26, 5.82, # r3, x3, y3
1.84, 8.16, # r4, y4
7.57, -1.20, -2.13, # 5
0.6, 0.86, 5.77, # 6
0.61, -1.37, 6.04 # 7
]
(x1, y1, x2, y2, r3, x3, y3, r4, y4, r5, x5, y5, r6, x6, y6, r7, x7, y7) = res[1]
@printf("外円 = %.15g, 天円 = %.15g, 地円 = %.15g のとき，人円 = %.15g\n", 2R, 2r1, 2r2, 2r3)
plot()
circle(0, 0, R, :purple)
circle(x1, y1, r1)
circle(x2, y2, r2, :blue)
circle(x3, y3, r3, :green)
circle(0, y4, r4, :magenta)
circle(x5, y5, r5, :darkorange)
circle(x6, y6, r6, :tomato)
circle(x7, y7, r7, :skyblue)
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, "外:R,(0,0)", :purple, :center, :bottom, delta=delta)
point(x1, y1, "天:r1\n(x1,y1)", :red, :center, :bottom, delta=delta)
point(x2, y2, "地:r2\n(x2,y2)", :blue, :center, :bottom, delta=delta)
point(x3, y3, "人:r3,(x3,y3)", :green, :center, delta=-3delta)
point(0, y4, "東:r4\n(0,R-r4)", :magenta, :center, :bottom, delta=delta)
point(x5, y5, "西:r5,(x5,y5)", :darkorange, :center, :bottom, delta=delta)
point(x6, y6, "南:r6\n(x6,y6)", :black, :left, :vcenter, deltax=3delta)
point(x7, y7, "北:r7\n(x7,y7)", :black, :right, :vcenter, deltax=-2delta)
point(0, R, "R", :purple, :center, :bottom, delta=delta)
end
end;
draw(20/2, 5/2, 4/2, true)
外円 = 20, 天円 = 5, 地円 = 4 のとき，人円 = 0.958639236453224
以下のアイコンをクリックして応援してください
算額（その1995）
『雑題』三十（大西佐兵衛）
愛媛和算研究会：『雑題』三十巻（大西佐兵衛著）を読むにあたって（その2），令和5年7月。
https://ehimewasan.com/wp-content/uploads/2024/09/72a03830eef2b8f8be223d0e0ff00b5e-1.pdf
キーワード：円7個
#Julia #SymPy #算額 #和算 #数学
雑題【11の4】日円，月円，火円，土円が互いに接している。その隙間に，水円，木円，金円を容れる。火円，水円，木円，金円の直径が 17.16 寸，3.9 寸，2.86 寸，4.29 寸のとき，土円の直径はいかほどか。
月円の中心を座標原点に置き，木円の中心が \(y\) 軸上になるように配置する。
日円の半径と中心座標を \(r_1,\ (x_1,\ y_1)\)
月円の半径と中心座標を \(r_2,\ (0,\ 0)\)
火円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
水円の半径と中心座標を \(r_4,\ (x_4,\ y_4)\)
木円の半径と中心座標を \(r_5,\ (0,\ r_2 + r_5)\)
金円の半径と中心座標を \(r_6,\ (x_6,\ y_6)\)
土円の半径と中心座標を \(r_7,\ (x_7,\ y_7)\)
とおき，以下の連立方程式の数値解を求める。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r1, x1, y1, r2, r3, x3, y3, r4, x4, y4,
r5, y5, r6, x6, y6, r7, x7, y7
y5 = r2 + r5
x5 = 0
eq1 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq2 = x3^2 + y3^2 - (r2 + r3)^2
eq3 = x7^2 + y7^2 - (r2 + r7)^2
eq4 = (x1 - x7)^2 + (y1 - y7)^2 - (r1 + r7)^2
eq5 = (x1 - x4)^2 + (y1 - y4)^2 - (r1 + r4)^2
eq6 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq7 = x4^2 + (y4 - y5)^2 - (r4 + r5)^2
eq8 = x6^2 + (y5 - y6)^2 - (r5 + r6)^2
eq9 = (x6 - x7)^2 + (y6 - y7)^2 - (r6 + r7)^2
eq10 = (x1 - x6)^2 + (y1 - y6)^2 - (r1 + r6)^2
eq11 = x1^2 + (y1 - y5)^2 - (r1 + r5)^2
eq12 = x4^2 + y4^2 - (r4 + r2)^2
eq13 = x6^2 + y6^2 - (r2 + r6)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13], (r1, x1, y1, r2, x3, y3, x4, y4, x6, y6, r7, x7, y7))
function H(u)
(r1, x1, y1, r2, x3, y3, x4, y4, x6, y6, r7, x7, y7) = u
return [
(x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2, # eq1
x3^2 + y3^2 - (r2 + r3)^2, # eq2
x7^2 + y7^2 - (r2 + r7)^2, # eq3
(x1 - x7)^2 + (y1 - y7)^2 - (r1 + r7)^2, # eq4
(x1 - x4)^2 + (y1 - y4)^2 - (r1 + r4)^2, # eq5
(x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2, # eq6
x4^2 + (y4 - (r2 + r5))^2 - (r4 + r5)^2, # eq7
x6^2 + (r2 + r5 - y6)^2 - (r5 + r6)^2, # eq8
(x6 - x7)^2 + (y6 - y7)^2 - (r6 + r7)^2, # eq9
(x1 - x6)^2 + (y1 - y6)^2 - (r1 + r6)^2, # eq10
x1^2 + (y1 - (r2 + r5))^2 - (r1 + r5)^2, # eq11
x4^2 + y4^2 - (r4 + r2)^2, # eq12
x6^2 + y6^2 - (r2 + r6)^2, # eq13
]
end;
(r3, r4, r5, r6) = (17.16, 3.9, 2.86, 4.29) ./ 2
iniv = BigFloat[14, 0.6, 22, 5, 13, 4, 3, 7, -4, 7, 13, -18, 3]
res = nls(H, ini=iniv)
([14.211137652235761, 0.6147384017903117, 22.463598083149616, 5.404545528578323, 13.496531066824165, 3.662125448852095, 3.3670139679392537, 6.538543941248187, -3.570333542168671, 6.651943782515173, 13.200000000000012, -18.404061703588194, 2.7238992520764294], true)
# r7: 土円の直径
2 * res[1][11]
26.400000000000023
土円の直径は 26.4 寸である。
術は以下のような式を示している。
(火径, 水径, 木径, 金径) = (17.16, 3.9, 2.86, 4.29)
天 = 火径*(水径 - 木径)
地 = 天*水径
人 = 天 - (金径 - 木径)*(火径 - 金径)*水径^2/金径^2
土 = 地/人
26.400000000000034
他の値を試しても数値解と同じ答えが得られる。
やはり，例示された図は与えられた数値（直径）を反映したものではないようだ。
描画関数プログラムのソースを見る
function draw(r3, r4, r5, r6, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, x1, y1, r2, x3, y3, x4, y4, x6, y6, r7, x7, y7) =res[1]
# [14.211137652235761, 0.6147384017903117, 22.463598083149616, 5.404545528578323, 13.496531066824165, 3.662125448852095, 3.3670139679392537, 6.538543941248187, -3.570333542168671, 6.651943782515173, 13.200000000000012, -18.404061703588194, 2.7238992520764294]
x5 = 0
y5 = r2 + r5
@printf("火円 = %.15g, 水円 = %.15g, 木円 = %.15g, 金円 = %.15g のとき，土円 = %.15g\n", 2r3, 2r4, 2r5, 2r6, 2r7)
plot()
circle(x1, y1, r1)
circle(0, 0, r2, :blue)
circle(x3, y3, r3, :green)
circle(x3, y3, r3, :magenta)
circle(x4, y4, r4, :darkorange)
circle(x5, y5, r5, :purple)
circle(x6, y6, r6, :tomato)
circle(x7, y7, r7, :blueviolet)
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(x1, y1, "日円:r1,(x1,y1)", :red, :center, delta=-delta)
point(0, 0, "月円:r2,(0,0)", :blue, :center, delta=-delta)
point(x3, y3, "火円:r3,(x3,y3)", :magenta, :center, delta=-delta)
point(x4, y4, "水円:r4,(x4,y4)", :darkorange, :left, :vcenter, delta=delta, deltax=5delta)
point(x5, y5, "木円:r5,(x5,y5)", :purple, :center, delta=-4delta)
point(x6, y6, "金円:r6,(x6,y6)", :tomato, :right, :vcenter, deltax=-4delta)
point(x7, y7, "土円:r7,(x7,y7)", :blueviolet, :center, delta=-delta)
end
end;
draw(17.16/2, 3.9/2, 2.86/2, 4.29/2, true)
火円 = 17.16, 水円 = 3.9, 木円 = 2.86, 金円 = 4.29 のとき，土円 = 26.4
以下のアイコンをクリックして応援してください
算額（その1994）
続神壁算法 付録〇三 東都芝増上寺境内内天満宮 寛政11年(1799)
豫州松山 丸山良玄門人 大西佐兵衛義全
平田浩一：大西佐兵衛著『雑題』の研究--第1巻より--，松山大学論集，第33巻第1号，2021年4月.
https://ehimewasan.com/wp-content/uploads/2024/12/6dba4f7ea9ead10da8d07ad8ddcf8d4a.pdf
キーワード：円2個，外円，正方形2個
#Julia #SymPy #算額 #和算 #数学
円弧の中に，大小の正方形と円を容れる。大正方形の一辺の長さが 11.31 寸，円の直径が 7.54 寸のとき，小正方形の一辺の長さはいかほどか。
この図形は算法助術の公式33 に取り上げられたものと基本的に同じであるが，平田の論文が最も丁寧でわかりやすい。
大正方形の辺が \(x-y\) 軸に並行になるように時計廻りに回転させて考える。このとき，必ず小正方形の対角線が\(x-y\) 軸に並行になり，大正方形の辺の延長線と円の交点が小正方形の頂点になる。
大正方形，小正方形の一辺の長さを \(a,\ b\)
外円の半径と中心座標を \(R,\ (0,\ 0)\)
円の半径と中心座標を \(r,\ (x,\ y)\)
とおくと，以下の関係式が成り立つ。
\( (\sqrt{2} - 1)r^2 - 2(a + b)r + \sqrt{2}a b = 0\)
この式は \(a,\ b,\ r\) を含むので，どれか 2 つの変数に実値を与えれば残りの変数の値が決定できる。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r, a, b
eq = (√Sym(2) - 1)*r^2 - 2(a + b)r + √Sym(2)a*b ⩵ 0
\(\sqrt{2} a b + r^{2} \left(-1 + \sqrt{2}\right) - r \left(2 a + 2 b\right) = 0\)
1. \(a,\ b\) を与えて \(r\) を求める
ans_r = solve(eq, r)[1] |> factor # 1 of 2
@show(ans_r)
ans_r = (1 + sqrt(2))*(a + b - sqrt(a^2 + sqrt(2)*a*b + b^2))
\(\left(1 + \sqrt{2}\right) \left(a + b - \sqrt{a^{2} + \sqrt{2} a b + b^{2}}\right)\)
ans_r(a => 19.6, b => 24.5).evalf()
\(8.00018413855404\)
2. \(b,\ r\) を与えて \(a\) を求める
ans_a = solve(eq, a)[1]
@show(ans_a)
ans_a = r*(2*b - sqrt(2)*r + r)/(sqrt(2)*b - 2*r)
\(\displaystyle \frac{r \left(2 b - \sqrt{2} r + r\right)}{\sqrt{2} b - 2 r}\)
ans_a(b => 47.95, r => 38.36/2).evalf()
\(57.2800011343656\)
3. \(a,\ r\) を与えて \(b\) を求める
大西佐兵衛が掲出した算額の場合である。
ans_b = solve(eq, b)[1]
@show(ans_b)
ans_b = r*(2*a - sqrt(2)*r + r)/(sqrt(2)*a - 2*r)
\(\displaystyle \frac{r \left(2 a - \sqrt{2} r + r\right)}{\sqrt{2} a - 2 r}\)
ans_b(a => 11.31, r => 7.54/2).evalf()
\(9.3900083909132\)
ans_r(a => 0.405417522627599, b => 0.359715193221452).evalf() |> println
0.140086998431592
4. 外円の半径
外円の半径 \(R\) は以下の関係式において，\(a,\ b\) を与えれば求まる。
@syms R
eq2 = 2R^2 ⩵ (a^2 + √Sym(2)a*b + b^2)
\(2 R^{2} = a^{2} + \sqrt{2} a b + b^{2}\)
ans_R = solve(eq2, R)[2] # 2 of 2
@show(ans_R)
ans_R = sqrt(2)*sqrt(a^2 + sqrt(2)*a*b + b^2)/2
\(\displaystyle \frac{\sqrt{2} \sqrt{a^{2} + \sqrt{2} a b + b^{2}}}{2}\)
ans_R(a => 0.405417522627599, b => 0.359715193221452).evalf()
\(0.5\)
描画関数プログラムのソースを見る
function draw(a, b, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
plot()
R = sqrt(2)*sqrt(a^2 + sqrt(2)*a*b + b^2)/2
circle(0, 0, R, :magenta)
y1 = a/2
x1 = sqrt(R^2 -y1^2)
(x2, y2) = (x1 - a, y1)
plot!([x1, x2, x2, x1, x1], [y1, y1, -y1, -y1, y1], color=:red, lw=0.5)
(x3, y3) = (x2 - b*cos(pi/4), y2 + b*sin(pi/4))
(x4, y4) = (x3 - b*sin(pi/4), y3 - b*cos(pi/4))
(x5, y5) = (x4 + b*cos(pi/4), y4 - b*sin(pi/4))
plot!([x2, x3, x4, x5, x2], [y2, y3, y4, y5, y2], color=:blue, lw=0.5)
r = (1 + sqrt(2))*(a + b - sqrt(a^2 + sqrt(2)*a*b + b^2))
y = a/2 + r
x = sqrt( (R - r)^2 -y^2)
circle(x, y, r, :green)
length = 0.2a
circle(x2, y2, length, :black, beginangle=225, endangle=270, lw=2)
segment(x4, y4, x2, y2, :red, lw=0.2)
plot!([0, x, x, 0], [0, y, 0, 0], color=:green, lw=0.2)
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(x1, y1, "(sqrt(R^2-a^2/4),a/2)", :black, :right, :bottom, delta=delta, deltax=5delta)
point(x2 - 1.1length*sin(pi/8), y2 - 1.1length*cos(pi/8), "45°", :black, :center, delta=-delta, mark=false)
point(x, y, "(sqrt( (R-r)^2-(a/2+r)^2),a/2+r", :black, :center, :bottom, delta=delta)
point(0, R, "R", :magenta, :center, :bottom, delta=delta)
end
end;
draw(0.405417522627599, 0.359715193221452, true)
以下のアイコンをクリックして応援してください
宝くじ--買う人は◯◯
はじめしゃちょー、年末ジャンボ宝くじを5700万円購入 衝撃の結末に机叩く
という記事がありましたけど。
たとえ，発行される宝くじ「全部を買ったとしても，戻ってくるのは5割以下」
買う前から結果は分かっているのよ。
算額（その216）解析解
長野県長野市元善町 善光寺 天保3年(1832)
中村信弥「改訂増補 長野県の算額」県内の算額2（134）
http://www.wasan.jp/zoho/zoho.html
長野県長野市鬼無里 松厳寺 天保10年(1839)
和算家 北明 寺島宗伴
http://www.wasan.jp/terasima/terasima.html
http://www.wasan.jp/terasima/terasima3.pdf
キーワード：円9個，外円，正三角形，斜線2本
#Julia #SymPy #算額 #和算 #数学
円内に正三角形と乙円 6 個，正三角形に 2 本の斜線を引き正三角形内に甲円 2 個を入れる。乙円の径を知って，甲円の径を求めよ。
外円の半径 \(R\) と甲円の半径 \(r_1\) の関係は，算法助術の公式58 が適用できる。
正三角形の外接円と内接円の半径を \(R,\ r_0\)，正三角形の一辺の長さを \(a\) とすれば，以下の式が成り立つ。
\(r_0 = R/2\)
\(a = \sqrt{3}R\)
\(r_0(4r_1^2 + a^2) = a^2\cdot 2r_1\)
using SymPy
@syms R::positive, r0, a, r1, r2::positive
r0 = R/2
a = √Sym(3)*R/2*2
eq1 = r0*(4r1^2 + a^2) - a^2*2r1
\(\displaystyle - 6 R^{2} r_{1} + \frac{R \left(3 R^{2} + 4 r_{1}^{2}\right)}{2}\)
eq1 を \(r_1\) について解く（最終的に \(r_1\) を求めることになる）。
ans_r1 = solve(eq1, r1)[1]
\(\displaystyle \frac{R \left(3 - \sqrt{6}\right)}{2}\)
外円の半径 \(R\) と乙円の半径 \(r_2\) の関係は，算法助術の公式48 が適用できる。
2つの円の反対側にある矢を \(y\) とすれば，以下の式が成り立つ。
\(y = 3R/2\)
\(2r_2 y + 2r_2^2 = 2\sqrt{2r_2^2 R y}\)
eq2 = 2r2*3R/2 + 2r2^2 - 2sqrt(2r2^2*R*3R/2)
\(- 2 \sqrt{3} R r_{2} + 3 R r_{2} + 2 r_{2}^{2}\)
eq2 を \(R\) について解く（\(r_2\) が与えられるので，計算により \(R\) が求まる）。
ans_R = solve(eq2, R)[1] |> simplify
\(\displaystyle \frac{2 r_{2} \left(3 + 2 \sqrt{3}\right)}{3}\)
求められた \(R\) を，先に示した \(ans_{r1}\) の式に代入すれば，甲円の半径 \(r_1\) が求まる。
ans_r1(R => 2r2*(3 + 2√Sym(3))/3)
\(\displaystyle \frac{r_{2} \left(3 + 2 \sqrt{3}\right) \left(3 - \sqrt{6}\right)}{3}\)
乙円の半径 \(r_2 = 1\) のとき，甲円の半径は \(r_1 = 1.18618474760839\) である。
ans_r1(R => 2r2*(3 + 2√Sym(3))/3)(r2 => 1).evalf()
\(1.18618474760839\)
図形を描くためには他のパラメータも決定しなければならない。
算額（その216）を参照のこと
以下のアイコンをクリックして応援してください
算額（その1992）
埼玉県加須市 秋葉山本坊(秋葉宮) 明治27年(1894)
和算の館
http://www.wasan.jp/saitama/akibahonbo2.html
キーワード：円11個，外円
#Julia, #SymPy, #算額, #和算
全円の中に，大円 2 個，中円 2 個，小円 6 個を容れる。小円の直径が 1 寸のとき，全円の直径はいかほどか。
全円の半径と中心座標を \(R,\ (0,\ 0)\)
大円の半径と中心座標を \(r_1,\ (0,\ R - r_1)\)
中円の半径と中心座標を \(r_2,\ (0,\ R - r_2)\)
小円の半径と中心座標を \(r_3,\ (R - r_3,\ 0),\ (x_3,\ y_3)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, r1::positive, r2::positive,
r3::positive, x3::positive, y3::positive
eq1 = R - (r1 + r2)
eq2 = x3^2 + (R - r2 - y3)^2 - (r2 + r3)^2
eq3 = x3^2 + (R - r1 - y3)^2 - (r1 - r3)^2
eq4 = x3^2 + (y3 + R - r1)^2 - (r1 + r3)^2
eq5 = (R - r3)^2 + (R - r1)^2 - (r1 + r3)^2;
res = solve([eq1, eq2, eq3, eq4, eq5], (R, r1, r2, x3, y3))[1]
(5*r3, 10*r3/3, 5*r3/3, 4*sqrt(3)*r3/3, 2*r3)
# R: 全円の半径
ans_R = res[1]
\(5 r_{3}\)
全円の半径は小円の半径の 5 倍である。
小円の直径が 1 寸のとき，全円の半径は 5 寸である。
# r1: 大円の半径
ans_r1 = res[2]
\(\displaystyle \frac{10 r_{3}}{3}\)
# r2: 中円の半径
ans_r2 = res[3]
\(\displaystyle \frac{5 r_{3}}{3}\)
# x3: 小円の中心 x 座標
ans_x3 = res[4]
\(\displaystyle \frac{4 \sqrt{3} r_{3}}{3}\)
# y3: 小円の中心 y 座標
ans_y3 = res[5]
\(2 r_{3}\)
描画関数プログラムのソースを見る
function draw(r3, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r1, r2, x3, y3) = (5*r3, 10*r3/3, 5*r3/3, 4*sqrt(3)*r3/3, 2*r3)
plot()
circle(0, 0, R, :magenta)
circle22(0, R - r1, r1)
circle22(0, R - r2, r2, :blue)
circle2(R - r3, 0, r3, :green)
circle4(x3, y3, r3, :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(0, R - r1, "大円:r1,(0,R-r1)", :red, :center, delta=-1.5delta)
point(0, R - r2, "中円:r2,(0,R-r2)", :blue, :center, delta=-delta)
point(R - r3, 0, "小円:r3\n(R-r3,0)", :green, :center, delta=-delta)
point(x3, y3, "小円:r3\n(x3,y3)", :green, :center, delta=-delta)
point(0, R, "R", :magenta, :center, :bottom, delta=delta)
end
end;
draw(1/2, true)
以下のアイコンをクリックして応援してください
算額（その1991）
埼玉県加須市 秋葉山本坊(秋葉宮) 明治27年(1894)
和算の館
http://www.wasan.jp/saitama/akibahonbo2.html
キーワード：円12個，外円，1/3円2個，正三角形
#Julia, #SymPy, #算額, #和算
全円の中に正三角形を描き，甲円 3 個，乙円 3 個，大円 1 個，中円 2 個，小円 2 個を容れる。小円の直径が 1 寸のとき，甲円の直径はいかほどか。
全円の半径と中心座標を \(R,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (0,\ r_1 - R)\)
乙円の半径と中心座標を \(r_2,\ (0,\ 2r_2 + (R - 2r_4 - r_5)\sin(\pi/6))\)
大円の半径と中心座標を \(r_3,\ (0,\ 0); r_3 = 2r_4\)
中円の半径と中心座標を \(r_4,\ (0,\ r_4); r_4 = r_1\)
小円の半径と中心座標を \(r_5,\ (0, r_3 - r_5)\)
円弧の半径と中心座標を \(2r_4,\ (0,\ 2r_4)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, r1::positive, r2::positive,
r3::positive, x3::positive, r4::positive, r5::positive
r4 = r1 # = R/4
r3 = 2r4
eq0 = r1 - R/4
eq1 = r2 - (R - 2r4 - r2)*sin(PI/6)
eq2 = (2r4 - r5)^2 + 4r4^2 - (2r4 + r5)^2
res = solve([eq0, eq1, eq2], (r1, r2, R))[1]
(2*r5, 4*r5/3, 8*r5)
# r1: 甲円の半径
ans_r1 = res[1]
@show(ans_r1)
ans_r1 = 2*r5
\(2 r_{5}\)
甲円の半径 \(r_1\) は，小円の半径 \(r_5\) の 2 倍である。
小円の直径が 1 寸のとき，甲円の直径は 2 寸である。
# r2: 乙円の半径
ans_r2 = res[2]
@show(ans_r2)
ans_r2 = 4*r5/3
\(\displaystyle \frac{4 r_{5}}{3}\)
# R: 全円の半径
ans_R = res[3]
@show(ans_R)
ans_R = 8*r5
\(8 r_{5}\)
描画関数プログラムのソースを見る
function draw(r5, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = 8r5
r2 = 4*r5/3
r1 = r4 = R/4
r3 = 2r4
@printf("小円の直径が %.15g のとき，甲円の直径は %.15g である\n", 2r5, 2r1)
plot([R*cosd(30), 0, -R*cosd(30), R*cosd(30)], [-R*sind(30), R, -R*sind(30), -R*sind(30)], color=:green, lw=0.5)
circle(0, 0, R, :magenta)
rotate(0, r1 - R, r1)
circle22(0, r4, r4)
circle(0, 0, r3, :darkorange)
rotate(0, 2r1 + r2, r2, :blue)
circle(0, 2r4, 2r4, :darkcyan, beginangle=210, endangle=330)
circle(0, -2r4, 2r4, :darkcyan, beginangle=30, endangle=150)
circle2(r3 - r5, 0, r5, :purple)
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, r1 - R, "甲円:r1\n(0,r1-R)", :red, :center, delta=-delta)
point(0, -r4, "中円:r4\n(0,-r4)", :red, :center, delta=-delta)
point(0, r4, "中円:r4\n(0,r4)", :red, :center, delta=-delta)
point(0, 0, "大円:2r4\n(0,0)", :darkorange, :center, delta=-delta)
point(0, 2r1 + r2, "乙円:r2\n(0,2r1+r2)", :blue, :center, :bottom, delta=delta)
point(r3 - r5, 0, "小円:r5,(2r4-r5,0)", :purple, :left, delta=-delta, deltax=-2delta)
point(0, 2r4, "円弧:2r4\n(0,2r4)", :darkcyan, :center, delta=-delta)
point(0, R, "R", :magenta, :center, :bottom, delta=delta)
end
end;
draw(1/2, true)
小円の直径が 1 のとき，甲円の直径は 2 である
以下のアイコンをクリックして応援してください
新年初の太陽黒点 2026/01/01
あけましておめでとうございます
本年初の太陽画像では，かなり活発な太陽黒点が観察されました
北の地域ではオーロラなども見られるのかな?
以下の写真は，香川県のソウルフード，白味噌仕立てのあんもち雑煮です
毎年，必ずと言っていいほど，どこかのテレビ局（キー局）で取り上げられるので，ご存知の方も多いかと思います
あんこの入った餅（焼かずに）と，大根，金時人参（香川特産），青のり（今日は，青のりがなかったので水菜を散らしました）を入れます
昔，韓国からの留学生に勧めたら，あと一息で食べてもらえるところを，一緒に来ていた留学生の奥さんに「残したら悪いから，やめといたほうがいいのでは？」と止められて，食べてもらえなかったという逸話｛？）があります
ネット通販でも取り寄せできるので，興味のある方はぜひどうぞ
算額（その1990）
埼玉県加須市 秋葉山本坊(秋葉宮) 明治27年(1894)
和算の館
http://www.wasan.jp/saitama/akibahonbo2.html
キーワード：円14個，外円
#Julia, #SymPy, #算額, #和算
全円の中に，甲円 2 個，乙円 3 個，丙円 2 個，丁円 6 個 を容れる。丁円の直径が 1 寸のとき，乙円，丙円の直径はいかほどか。
全円の半径と中心座標を \(R,\ (0,\ 0),\ (x_1,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (0,\ r_2)\)
乙円の半径と中心座標を \(r_2,\ (0,\ 0),\ (0,\ 2r_2)\)
丙円の半径と中心座標を \(r_3,\ (3r_2 - r_3,\ 0)\)
丁円の半径と中心座標を \(r_4,\ (R - r_4,\ 0)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, X::positive, r0::positive, r1::positive, r2::positive,
r3::positive, x3::positive, r4::positive
r1 = 2r2
r0 = 3r2
R = r0 + 2r4
eq1 = (r0 - r3)^2 + r2^2 - (r1 + r3)^2
eq2 = (R + 3r2) - 2R*cos(PI/6); # x1
res = solve([eq1, eq2], (r2, r3))[1]
(2*r4/3 + 2*sqrt(3)*r4/3, 2*r4/5 + 2*sqrt(3)*r4/5)
# r2: 乙円の半径
ans_r2 = res[1] |> simplify
@show(ans_r2)
ans_r2 = 2*r4*(1 + sqrt(3))/3
\(\displaystyle \frac{2 r_{4} \left(1 + \sqrt{3}\right)}{3}\)
乙円の半径 \(r_2\) は，丁円の半径の \(2(1 + \sqrt{3})/3\) 倍である。
丁円の直径が 1 寸のとき，乙円の直径は 1.82136720504592 寸である。
# 乙円の直径
2 * ans_r2(r4 => 1/2).evalf() |> println
1.82136720504592
# r3: 丙円の半径
ans_r3 = res[2] |> simplify
@show(ans_r3)
ans_r3 = 2*r4*(1 + sqrt(3))/5
\(\displaystyle \frac{2 r_{4} \left(1 + \sqrt{3}\right)}{5}\)
丙円の半径 \(r_2\) は，丁円の半径の \(2(1 + \sqrt{3})/5\) 倍である。
丁円の直径が 1 寸のとき，乙円の直径は 1.09282032302755 寸である。
# 丙円の直径
2 * ans_r3(r4 => 1/2).evalf()
\(1.09282032302755\)
描画関数プログラムのソースを見る
function draw(r4, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, r3) = (2*r4/3 + 2*sqrt(3)*r4/3, 2*r4/5 + 2*sqrt(3)*r4/5)
r1 = 2r2
r0 = 3r2
R = r0 + 2r4
x1 = R + 3r2
plot()
circle(0, 0, r2)
circle22(0, 2r2, r2)
circle22(0, r2, 2r2, :blue)
circle(0, 0, r0, :green)
circle2(r0 - 3r2/5, 0, 3r2/5, :magenta)
circle(0, 0, R, :brown)
rotate(x1, 0, R, :brown, angle=60)
rotate(r0 + r4, 0, r4, :darkorange, angle=60)
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(x1, 0, "x1", :brown, :center, delta=-delta)
point(R*cos(pi/6), R*sin(pi/6), "R*cos(pi/6),R*sin(pi/6)", :brown, :left, delta=-delta/2, deltax=delta)
point(0, 0, "全", :brown)
point(0, r2, "甲", :blue)
point(0, 2r2, "乙", :red)
point(r0 - r3, 0, "丙", :magenta)
point(R - r4, 0, "丁", :darkorange)
point(0, R, "R", :brown, :center, :bottom, delta=delta)
plot!(xlims=(-1.1R, 1.1(x1 + R)), ylims=(-1.1R, 1.1R))
end
end;
draw(1, true)
以下のアイコンをクリックして応援してください
算額（その1989）
埼玉県加須市 秋葉山本坊(秋葉宮) 明治27年(1894)
和算の館
http://www.wasan.jp/saitama/akibahonbo2.html
キーワード：円13個，外円，正方形，対角線
#Julia, #SymPy, #算額, #和算
全円の中に，正方形とその対角線を容れ，甲円 2 個，乙円 6 個，丙円 4 個を容れる。乙円の直径が 1 寸のとき，甲円，丙円の直径はいかほどか。
上のように，図形を時計方向に45°回転させると考えやすい。
全円の半径と中心座標を \(R,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (x_1,\ -r_1)\)
乙円の半径と中心座標を \(r_2,\ (0,\ r_2),\ (2r_2,\ r_2),\ (4r_2,\ r_2)\)
丙円の半径と中心座標を \(r_3,\ ( (R - r_3)/\sqrt{2},\ (R - r_3)/\sqrt{2})\)
とおき，以下の連立方程式を解く。
eq1, eq2 は独立で，かつ，簡単なので，連立方程式に組み入れなくてもよい。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, r1::positive, x1::positive, r2::positive, r3::positive
eq1 = 2r3 - (R - R/√Sym(2))
\(\displaystyle - R + \frac{\sqrt{2} R}{2} + 2 r_{3}\)
eq2 = r2*(8 + 2tand(Sym(675)//10)) - 2R
\(- 2 R + r_{2} \left(2 \sqrt{2} + 10\right)\)
eq3 = x1 - r1 + r2
eq4 = (R - x1)*tand(Sym(225)//10) - r1
\(- r_{1} + \left(-1 + \sqrt{2}\right) \left(R - x_{1}\right)\)
res = solve([eq1, eq2, eq3, eq4], (R, r1, x1, r3))
\begin{equation*}\begin{cases}r_{1} & \text{=>} &- 2 \sqrt{2} r_{2} + 5 r_{2}\\r_{3} & \text{=>} &- \frac{3 \sqrt{2} r_{2}}{4} + 2 r_{2}\\R & \text{=>} &\frac{3 \sqrt{2} r_{2}}{-2 + \sqrt{2}} - \frac{8 r_{2}}{-2 + \sqrt{2}}\\x_{1} & \text{=>} &- 2 \sqrt{2} r_{2} + 4 r_{2}\\\end{cases}\end{equation*}
# r1: 甲円の半径
ans_r1 = res[r1] |> simplify |> factor
@show(ans_r1)
ans_r1 = -r2*(-5 + 2*sqrt(2))
\(- r_{2} \left(-5 + 2 \sqrt{2}\right)\)
甲円の半径 \(r_1\) は，乙円の半径 \(r_2\) の \(5 - 2\sqrt{2}\) 倍である。
乙円の直径が 1 寸のとき，甲円の直径は \(5 - 2\sqrt{2} = 2.17157287525381\) 寸である。
# 甲円の直径
2 * ans_r1(r2 => 1/2).evalf()
\(2.17157287525381\)
# r3: 丙円の半径
ans_r3 = res[r3] |> simplify |> factor
@show(ans_r3)
ans_r3 = -r2*(-8 + 3*sqrt(2))/4
\(\displaystyle - \frac{r_{2} \left(-8 + 3 \sqrt{2}\right)}{4}\)
丙円の半径 \(r_3\) は，乙円の半径 \(r_2\) の \( (8 - 3\sqrt{2})/4\) 倍である。
乙円の直径が 1 寸のとき，丙円の直径は \( (8 - 3\sqrt{2})/4 = 0.939339828220179\) 寸である。
# 丙円の直径
2 * ans_r3(r2 => 1/2).evalf()
\(0.939339828220179\)
# R: 全円の半径
ans_R = res[R] |> simplify |> factor
@show(ans_R)
ans_R = r2*(sqrt(2) + 5)
\(r_{2} \left(\sqrt{2} + 5\right)\)
全円の半径 \(R\) は，乙円の半径 \(r_2\) の \(\sqrt{2} + 5\) 倍である。
乙円の直径が 1 寸のとき，全円の直径は \(\sqrt{2} + 5 = 6.41421356237309\) 寸である。
# 全円の直径
2 * ans_R(r2 => 1/2).evalf()
\(6.41421356237309\)
# x1: 甲円の中心座標
ans_x1 = res[x1] |> simplify |> factor
@show(ans_x1)
ans_x1 = -2*r2*(-2 + sqrt(2))
\(- 2 r_{2} \left(-2 + \sqrt{2}\right)\)
ans_x1(r2 => 1/2).evalf()
\(0.585786437626905\)
描画関数プログラムのソースを見る
function draw(r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = r2*(5 - 2√2)
r3 = r2*(8 - 3√2)/4
R = r2*(√2 + 5)
x1 = r2*(4 - 2√2)
plot([R, 0, -R, 0, R], [0, R, 0, -R, 0], color=:brown, lw=0.5)
circle(0, 0, R, :magenta)
circle2(x1, -r1, r1)
circle(0, -r1, r2, :blue)
circle(0, r2, r2, :blue)
circle2(2r2, r2, r2, :blue)
circle2(4r2, r2, r2, :blue)
rotate( (R - r3)/√2, (R - r3)/√2, r3, :green, angle=90)
segment(-R, 0, R, 0, :black, lw=1)
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( (R - r3)/√2, (R - r3)/√2, "丙円:r3\n( (R-r3)/√2,(R-r3)/√2)", :green, :right, :vcenter, deltax=-10delta)
point(4r2, r2, "乙円:r2,(4r2,r2)", :blue, :right, delta=8delta)
point(x1, -r1, "甲円:r1\n(x1,-r1)", :red, :left, :vcenter, deltax=delta)
point(0, R, "R", :magenta, :center, :bottom, delta=delta)
end
end;
draw(1/2, true)
以下のアイコンをクリックして応援してください
算額（その1988）
埼玉県加須市 秋葉山本坊(秋葉宮) 明治27年(1894)
和算の館
http://www.wasan.jp/saitama/akibahonbo2.html
キーワード：円9個，外円
#Julia, #SymPy, #算額, #和算
全円の中に，大円 2 個，小円 6 個を容れる。小円の直径が 1 寸のとき，大円の直径はいかほどか。
図形的には算額（その1579）と同じで，求めるものが違うだけである。
全円の半径と中心座標を \(R,\ (0,\ 0)\)
大円の半径と中心座標を \(r_1,\ (0,\ R - r_1)\)
小円の半径と中心座標を \(r_2,\ (x_{21},\ 0),\ (x_{22},\ r_2)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, r1::positive, r2::positive, x21::positive, x22::positive
R = 2r1
eq1 = x22^2 + r2^2 - (R - r2)^2
eq2 = x21^2 + (R - r1)^2 - (r1 + r2)^2
eq3 = (x22 - x21)^2 + r2^2 - 4r2^2
res = solve([eq1, eq2, eq3], (r1, x21, x22))[1];
# r1: 大円の半径
ans_r1 = res[1]
@show(ans_r1)
ans_r1 = r2*(5 + sqrt(33))/4
\(\displaystyle \frac{r_{2} \left(5 + \sqrt{33}\right)}{4}\)
大円の半径 \(r1\) は，小円の半径の \(\displaystyle \frac{5 + \sqrt{33}}{4}\) 倍である。
小円の直径が 1 寸のとき，大円の直径は 2.68614066163451 である。
2 * ans_r1(r2 => 1/2).evalf()
\(2.68614066163451\)
描画関数プログラムのソースを見る
function draw(r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, x21, x22) = (r2*(5 + sqrt(33))/4, r2*(5*sqrt(2) + sqrt(66))/sqrt(3*sqrt(33) + 19), r2*sqrt(3*sqrt(33)/2 + 19/2))
R = 2r1
plot()
circle(0, 0, R, :green)
circle22(0, R - r1, r1)
circle2(x21, 0, r2, :blue)
circle4(x22, r2, r2, :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, R, "R", :green, :center, :bottom, delta=delta/2)
point(0, R - r1, "大円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
point(x21, 0, "小円:r2\n(x21,0)", :blue, :center, :bottom, delta=delta/2)
point(x22, r2, "小円:r2\n(x22,r2)", :blue, :center, :bottom, delta=delta/2)
end
end;
draw(1/2, true)
以下のアイコンをクリックして応援してください
算額（その1987）
九十五 群馬県尾島町世良田 八坂神社 安政5年(1858)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：楕円錐台，体積
#Julia, #SymPy, #算額, #和算
上面が円，下面が楕円の物体がある。上面の半径，下面の長半径と短半径，高さが与えられたとき，体積を求めよ。
この物体を上方から見て等高線図を描くと，一番上は円であるがその下はすべて同心楕円である。
上面の半径を \(r\)，下面の長半径と短半径を \(a,\ b\)，高さを \(h\) とすると，高さが \(x\) で水平に切断したときの切断面は長半径 \(a_x\)，短半径 \(b_x\) はそれぞれ以下のようなる。
\(a+x = a - (a - r)x/h\)，\(b_x = b - (b - r)x/h\) で，面積 \(S\) は \(\pi a_x b_x\) である。
体積 \(V\) は \(S\) を 0 から h まで積分することによって得られる。
\(V = \int_{0}^{h} S dx\)
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a, b, r, h, x, ax, bx
ax = a - (a - r)*x/h
bx = b - (b - r)*x/h
V = integrate(PI*ax*bx, (x, 0, h)) |> simplify
\(\displaystyle \frac{\pi h \left(2 a b + a r + b r + 2 r^{2}\right)}{6}\)
長半径，短半径，半径，高さが 5，2，1，10 のときは 151.843644923507 である。
V(a => 5, b => 2, r => 1, h => 10).evalf()
\(151.843644923507\)
術は，( (a + b + 2r)*2r + (2a * 2b)) * h * (pi/4)/3 に相当するので，同じものである。
@syms 長径, 短径, 円径, 高さ
V2 = ( *1* 高さ * (pi/4)/3;
V2(長径 => 10, 短径 => 4, 円径 => 2, 高さ => 10)
\(151.843644923507\)
描画関数プログラムのソースを見る
function draw(a, b, r, h, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
plot()
δ =(a - r)/5
c = 0
for x in 0:10
x2 = a - x*(a - r)/h
y2 = b - x*(b - r)/h
c += 1
ellipse(0, 0, x2, y2, color= c == 5 ? :red : c)
end
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", :black, :left, :bottom, delta=delta)
point(r, 0, "r ", :black, :right, :bottom, delta=delta)
point(0, b, "b", :black, :center, :bottom, delta=delta)
point(0, r, "r", :black, :center, delta=-delta)
end
end;
draw(5, 2, 1, 10, true)
以下のアイコンをクリックして応援してください
*1:長径 + 短径)/2 + 円径)*円径 + (長径 * 短径
