折れた楓
算額（その1420）
福島県二本松市亀谷 坂上観音堂 寛政5年(1793)令和3年復元奉納
http://www.wasan.jp/fukusima/sakauekanondo.html
キーワード：円5個，円弧
#Julia, #SymPy, #算額, #和算
円弧（弓形）の中に 5 個の円を入れる。そのうちの甲円，乙円，丙円の直径が 16 寸，25 寸，9 寸のとき，外円の直径はいかほどか。
甲円，乙円，丙円が右から順に一つおきになっているが，残りの円も右から丁円，丙円と名付ける。
外円（円弧，弓形を構成する円）の半径と中心座標を R, (0, 0)
円弧を構成する水平な弦と y 軸の交点座標を (0, y)
甲円の半径と中心座標を r1, (x1, y + r1)
乙円の半径と中心座標を r2, (x2, y + r2)
丙円の半径と中心座標を r3, (x3, y + r3)
丁円の半径と中心座標を r4, (x4, y + r4)
戊円の半径と中心座標を r5, (x5, y + r5)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R, y, r1, x1, r2, x2, r3, x3, r4, x4, r5, x5
eq1 = x1^2 + (y + r1)^2 - (R - r1)^2
eq2 = x2^2 + (y + r2)^2 - (R - r2)^2
eq3 = x3^2 + (y + r3)^2 - (R - r3)^2
eq4 = x4^2 + (y + r4)^2 - (R - r4)^2
eq5 = x5^2 + (y + r5)^2 - (R - r5)^2
eq6 = (x1 - x4)^2 +(r1 - r4)^2 - (r1 + r4)^2
eq7 = (x4 - x2)^2 +(r4 - r2)^2 - (r4 + r2)^2
eq8 = (x2 - x5)^2 +(r2 - r5)^2 - (r2 + r5)^2
eq9 = (x5 - x3)^2 +(r5 - r3)^2 - (r5 + r3)^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (R, y, x1, x2, x3, r4, x4, r5, x5))
using NLsolve
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 Float64.(v), r.f_converged
end;
function H(u)
(R, y, x1, x2, x3, r4, x4, r5, x5) = u
return [
x1^2 - (R - r1)^2 + (r1 + y)^2, # eq1
x2^2 - (R - r2)^2 + (r2 + y)^2, # eq2
x3^2 - (R - r3)^2 + (r3 + y)^2, # eq3
x4^2 - (R - r4)^2 + (r4 + y)^2, # eq4
x5^2 - (R - r5)^2 + (r5 + y)^2, # eq5
(r1 - r4)^2 - (r1 + r4)^2 + (x1 - x4)^2, # eq6
(-r2 + r4)^2 - (r2 + r4)^2 + (-x2 + x4)^2, # eq7
(r2 - r5)^2 - (r2 + r5)^2 + (x2 - x5)^2, # eq8
(-r3 + r5)^2 - (r3 + r5)^2 + (-x3 + x5)^2, # eq9
]
end;
(r1, r2, r3) = (16, 25, 9) ./ 2
iniv = BigFloat[70, 45, 33, -12, -43, 12, 13, 8, -31]
res = nls(H, ini=iniv)
([69.73157051282051, 43.72996794871795, 33.686751312703706, -10.660364339463198, -43.92070107858837, 12.139917695473251, 13.976922133962859, 8.642578125, -31.448074801416432], true)
# 直径の分数表示
rationalize(2*69.73157051282051)
87025//624
# 帯分数で表示するための計算 商と余りの計算
divrem(87025,624)
(139, 289)
外円の直径は 2*69.73157051282051 = 139.46314102564102 であり，分数で表すと 87025//624 である。
帯分数表示というのは小学校以来使うことはないが，139 + 289/624 であり，「答」の「外円径一百三十九寸六百二十四分之二百八十九寸とピタリ一致する（長精度演算のおかげ）。
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3) = (16, 25, 9) ./ 2
(R, y, x1, x2, x3, r4, x4, r5, x5) = [70, 45, 33, -12, -43, 12, 13, 8, -31]
(R, y, x1, x2, x3, r4, x4, r5, x5) = [69.73157051282051, 43.72996794871795, 33.686751312703706, -10.660364339463198, -43.92070107858837, 12.139917695473251, 13.976922133962859, 8.642578125, -31.448074801416432]
x =sqrt( R^2 - y^2)
θ = atand(y, x)
plot()
circle(0, 0, R, beginangle=θ, endangle=180 - θ, n=500)
circle(x1, y + r1, r1, :blue)
circle(x2, y + r2, r2, :magenta)
circle(x3, y + r3, r3, :green)
circle(x4, y + r4, r4, :orange)
circle(x5, y + r5, r5, :tomato)
plot!([-x, 0, x, -x], [y, 0, y, y], color=:gray70, 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(x1, y + r1, "甲円:r1,(x1,y+r1)", :blue, :center, delta=-8.5delta)
point(x2, y + r2, "乙円:r2\n(x2,y+r2)", :magenta, :center, delta=-delta)
point(x3, y + r3, "丙円:r3,(x3,y+r3)", :green, :left, delta=-5delta, deltax=-5delta)
point(x4, y + r4, "丁円:r4\n(x4,y+r4)", :orange, :center, delta=-delta)
point(x5, y + r5, "戊円:r5\n(x5,y+r5)", :tomato, :center, :bottom, delta=delta)
point(0, y, " y", :black, :left, delta=-delta)
point(0, R, "R", :red, :center, :bottom, delta=delta)
end
end;
draw(true)
綾川町小野 セルフ こだわり麺や 綾南店
県内に十数店舗ある
いつもの冷かけととり天
算額（その1419）
本庄市都島 正観寺 享保11年(1126)
https://www.honjo-kanko.jp/sightseeing/shokanji.html
https://tamalotus2.exblog.jp/24658561/
https://yamabukiwasan.sakura.ne.jp/page3.html#shoukan
山口正義：やまぶき 第26号
https://yamabukiwasan.sakura.ne.jp/ymbk26.pdf
山口正義：やまぶき 第40号
https://yamabukiwasan.sakura.ne.jp/ymbk40.pdf
参考古三面
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：立方体，体積，等差数列
#Julia, #SymPy, #算額, #和算
「現在（当時）現存第三古」ということであったが，2024年11月現在では「埼玉県内最古，全国で5番目に古い古額」だそうである。「埼玉の算額」にも出ていない。
甲，乙，丙，丁，戊 の立方体がある。甲，乙の立方体の体積の和は 189，丙，丁，戊の立方体の体積の和は 36 である。各立方体の一辺の長さは等差数列である。それぞれの一辺の長さはいかほどか。
すでに「算額（その388）」で解いたが，補足的な内容も含め新たに記事を立てる。
https://blog.goo.ne.jp/r-de-r/e/0ab004aeda635279911a670ca5ef16ef
各立方体の一辺の長さを，甲, 乙, 丙, 丁 とする。
甲，乙の立方体の体積の和を K1，丙，丁，戊の立方体の体積の和を K2，公差 を d とおき，以下の連立方程式を解く。
K1, K2 を変数のまま解くのは SymPy では無理のようである。
K1, K2 に実際の値を与えれば，簡単に解が求まる。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms 甲::positive, 乙::positive, 丙::positive, 丁::positive, 戊::positive, K1::positive, K2::positive, d::positive
(K1, K2) = (189, 36)
丁 = 戊 + d
丙 = 戊 + 2d
乙 = 戊 + 3d
甲 = 戊 + 4d
eq1 = 甲^3 + 乙^3 - K1
eq2 = 丙^3 + 丁^3 + 戊^3 - K2;
solve([eq1, eq2], (戊, d))[1] # 1 of 3
(1, 1)
戊 = 1，公差 = 1 なので，甲 = 5, 乙 = 4, 丙 = 3, 丁 = 2 である。
---
これではあまりにもあっけないので，SymPy が実際はどのようにして解を求めているのはわからないが，手作業で（SymPy の助けを借りながら）解いてみよう。
まず，eq1, eq2 がどのように表されているのかを見てみる。
eq1 |> expand |> println
eq2 |> expand |> println
91*d^3 + 75*d^2*戊 + 21*d*戊^2 + 2*戊^3 - 189
9*d^3 + 15*d^2*戊 + 9*d*戊^2 + 3*戊^3 - 36
d, 戊はともに 3 次式で，これが方程式を解くのを困難にしている一因であろう。そこで，定数倍した式を引き算することで，2 次式に字数を落とす。
eq3 = 9eq1 - 91eq2 として，これを解き d を求める。
eq3 = 9eq1 - 91eq2 |> expand |> factor
eq3 |> println
-15*(46*d^2*戊 + 42*d*戊^2 + 17*戊^3 - 105)
2 つの解が得られるが，2 つ目のものが適解である。
ans_d = solve(eq3, d)[2] # 2 of 2
ans_d |> println
-21*戊/46 + sqrt(4830 - 341*戊^3)/(46*sqrt(戊))
eq1 に d を代入し，eq4 とする。
eq4 = eq1(d => ans_d) |> simplify |> numerator
eq4 |> println
359359*戊^(9/2) - 14711697*戊^(3/2) - 5551*戊^3*sqrt(4830 - 341*戊^3) + 219765*sqrt(4830 - 341*戊^3)
これを解き 戊 を求める。
res = solve(eq4, 戊)[1]
res |> println
1
戊は 1 である。
遡って，ans_d の戊に 1 を代入して d = 1 を得る。
ans_d(戊 => 1) |> println
1
よって，甲，乙，丙，丁，戊は 5, 4, 3, 2, 1 である。
PVアクセスランキング にほんブログ村
