裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

算額(その317)

2023年07月07日 | Julia

算額(その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)

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その316)

2023年07月07日 | Julia

算額(その316)

(4) 京都府長岡京市天神 長岡天満宮 寛政2年(1970)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

森田健(2020): 日本文化としての数学:和算と算額, 日本語・日本文化,2020,47,p81-107.
https://ir.library.osaka-u.ac.jp/repo/ouka/all/75881/JLC_47_081.pdf

p02-03 特集「算額」(杉浦).indd
https://www.s-coop.net/lifestage/backnumber/2011/pdf/1106_02-03.pdf

キーワード:円2個,正三角形,直線上

大円と小円の間に正三角形が挟まっている。大円と小円の直径はそれぞれ 18cm,8cm である。正三角形の辺の長さを求めよ。

大円の半径を r1,中心座標を (0, r1)
小円の半径を r2,中心座標を (x2, r2)
正三角形の一辺の長さを 2a,底辺の中点の x 座標を x1 とする。
以下の連立方程式を解く。
なお,eq1 は直線の上にある外接する 2 つの円の中心間の距離を表すもので,距離 = x2 = 2sqrt(r1, r2) である。eq2, eq3 はそれれぞれの円の中心から正三角形の斜辺までの距離についての方程式である。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::positive, x1::positive, a::positive;

(r1, r2) = (9, 4)
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = distance(x1 - a, 0, x1, sqrt(Sym(3))a, 0, r1) - r1^2
eq3 = distance(x1 + a, 0, x1, sqrt(Sym(3))a, x2, r2) - r2^2;
res = solve([eq1, eq2, eq3], (x1, a, x2))

   3-element Vector{Tuple{Sym, Sym, Sym}}:
    (6 - 5*sqrt(3)/2, 6 + 13*sqrt(3)/2, 12)
    (6 + 7*sqrt(3)/2, sqrt(3)/2 + 6, 12)
    (5*sqrt(3)/6 + 6, 6 - 13*sqrt(3)/6, 12)

三番目の組が適解である。

(5*sqrt(3)/6 + 6, 6 - 13*sqrt(3)/6, 12)

   (7.4433756729740645, 2.247223250267433, 12)

   正三角形の一辺の長さ = 4.494446500534866 cm

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (9, 4)
   x1 = 5*sqrt(3)/6 + 6
   a = 6 - 13*sqrt(3)/6
   x2 = 12
   println("正三角形の一辺の長さ = $(2a) cm")
   plot()
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   plot!([x1 - a, x1 + a, x1, x1 - a], [0, 0, √3a, 0], color=:green, lw=0.5)
   if more
       point(0, r1, " r1", :red)
       point(x2, r2, "(x2,r2)", :blue)
       point(x1, 0, "x1", :green, :left, :bottom)
       point(x1 - a, 0, "x1-a", :green, :right, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村