算額(その317)
p02-03 特集「算額」(杉浦).indd
https://www.s-coop.net/lifestage/backnumber/2011/pdf/1106_02-03.pdf
御香宮の算額
甲,乙,丙の 3 つの正方形は,以下の条件を満たす。
1. 各面積の和はある一定の値を取る
2. 甲の一辺の 3 乗根,乙の一辺の 5 乗根,丙の一辺の 7 乗根の和はある一定の値を取る
3. 甲と乙の一辺の長さの差と,乙と丙の一辺の長さの差は等しい
このとき,甲,乙,丙それぞれの辺の長さを求めよ。
御香宮神社の算額
http://tadahikostar.blog21.fc2.com/blog-entry-3847.html
には,鮮明な写真が掲載されている。
読んだところでは,
1. 甲乙丙の各平積の和が既知,
2. 甲の開立方,乙の開四乗方,丙の開六乗方の和が既知
3. 甲,乙,丙の一辺の長さは等差数列(甲>乙>丙)
で,条件 2. が違う。
また,
御香宮(京都・伏見)の算額の問題 - 高精度計算サイト - Keisan
https://keisan.casio.jp/exec/user/1465639202
では,
2. 甲の1辺の3乗と乙の1辺の5乗と丙の1辺の7乗がある値をとり
と書かれている。これは明らかに誤りである。
「開立」は立方根を表す。三乗は「再自乗」という。
以下では,「甲の開立方,乙の開四乗方,丙の開六乗方の和」を採用する。
既知である値を記号(変数)として解くのは solve() でも無理である。
また,既知の値を具体的に与えても solve() では解けないので,nlsolve() で数値解を求める。
具体的に,既知の値を求めておく。甲,乙,丙の辺の長さをそれぞれ 6, 4, 2 として,条件 1, 2 で表される数値を計算する。
(甲, 乙, 丙) = (6, 4, 2)
甲^2 + 乙^2 + 丙^2 |> println # a とする
甲^(1/3) + 乙^(1/4) + 丙^(1/6) |> println # b とする
56
4.353796203514608
丙^2 + 乙^2 + 甲^2 - 56, # eq1
丙^0.166666666666667 + 乙^0.25 + 甲^0.333333333333333 - 4.35379620351461, # eq2
丙 - 2*乙 + 甲, # eq3
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=1e-14)
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
v = r.zero
end
return v, r.f_converged
end;
function H(u)
(甲, 乙, 丙) = u
return [
丙^2 + 乙^2 + 甲^2 - 56, # eq1
丙^0.166666666666667 + 乙^0.25 + 甲^0.333333333333333 - 4.35379620351461, # eq2
丙 - 2*乙 + 甲, # eq3
]
end;
初期値により,二通りの解が示される。両方ともに,条件を満たしている。
iniv = [big"7.0", 4.1, 1.1]
res = nls(H, ini=iniv);
println(res);
(BigFloat[5.999999999999978421814505360998991236330686683731818877458530801081713663484236, 4.000000000000010789092747335389143058164701865210416514719281049866373768067708, 2.000000000000043156370989309779294879998717046689014151980031298601986612507766], true)
(甲, 乙, 丙) = res[1] # 6, 4, 2
丙^2 + 乙^2 + 甲^2 - 56, # eq1
丙^0.166666666666667 + 乙^0.25 + 甲^0.333333333333333 - 4.35379620351461, # eq2
丙 - 2*乙 + 甲 # eq3
(2.566627137897993161354744300647618625279842621939623293996006529315262445625456e-25, -3.867992033135278266328028707026632666657885527869404966426262217141461630526278e-27, -4.904726014337897321140092172924707876804539997277609351100695166069929253914544e-65)
iniv = [big"3.0", 2.1, 1.2]
res = nls(H, ini=iniv);
println(res);
(BigFloat[4.417532066579428278706241618180857546089192311716748510275041555425975958939494, 4.319756156028245623996284797567623698490811039184724945910874439881535311042085, 4.221980245477062969286327976954389850892429766652701381546707324388269228906724], true)
(甲, 乙, 丙) = res[1] # 4.417532066579429, 4.319756156028245, 4.221980245477063
丙^2 + 乙^2 + 甲^2 - 56, # eq1
丙^0.166666666666667 + 乙^0.25 + 甲^0.333333333333333 - 4.35379620351461, # eq2
丙 - 2*乙 + 甲 # eq3
(7.928474500472361481718927810447226817906506437873289558404530451020621884955364e-25, -5.071109598580144750966960369570845052151000656128305908528664840899463931562674e-27, 5.117456576204827560281226639245275222567489206961363157351246079614929159341669e-65)