裏 RjpWiki

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

算額(その1123)

2024年07月05日 | Julia

算額(その1123)

二十八 一関市萩荘 赤萩観音寺 天保2年(1831)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html
キーワード:円6個,長方形

長方形の中に甲円 1 個,乙円 1 個,丙円 2 個,丁円 2 個を容れる。長方形の短辺の長さが 1 寸のとき,長辺の長さはいかほどか。

「岩手の算額」も,解けるものは解いてきたが,解く前から「こりゃー,やってみても解けるかどうかわからない」というのが残っている。この算額の問題も,図を見ただけで,「なんだこりゃ」と後回しにしてきた。こんなの算額の図じゃない。案の定,以下のような図になるものであった。しかも,これもこの本に載っている問題の常で,答えが怪しい。算額の元々がおかしいのか,山村の解説(?)も首をひねるものが多い。

このシリーズでは,得られた解で図を描いて提示している。図が正しく描かれているということは,解も正しいという一つの証拠ではあろう。
多くの算額問題では,条件を与えて結果としてパラメータの一つを求めるというのが常道である。しかし,一つのパラメータ以外のパラメータも含めて図を描いて検証するということは殆ど行われてこなかったのではないかと思う次第である。

長方形の長辺,短辺を a, 2r1
甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (r2, r1)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (r4, r4)
とおき,以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms a::positive, r1::positive, r2::positive,
     r3::positive, x3::positive,
     r4::positive
eq1 = (a - r1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq2 = (r2 - x3)^2 + (r1 - r3)^2 - (r2 + r3)^2
eq3 = (r2 - r4)^2 + (r1 -r4)^2 - (r2 + r4)^2
eq4 = (x3 - r4)^2 + (r3 - r4)^2 - (r3 + r4)^2
eq5 = 2r1 + 2r2 - a
res = solve([eq1, eq2, eq3, eq4, eq5], (a, r2, r3, x3, r4))[1]

   (4*r1*(-6*sqrt(239 + 169*sqrt(2)) + 4 + 3*sqrt(2) + 4*sqrt(478 + 338*sqrt(2))), r1*(-12*sqrt(239 + 169*sqrt(2)) + 7 + 6*sqrt(2) + 8*sqrt(478 + 338*sqrt(2))), 2*r1*(-145*sqrt(239 + 169*sqrt(2)) - 100*sqrt(478 + 338*sqrt(2)) + 3131 + 2214*sqrt(2))/(-44*sqrt(239 + 169*sqrt(2)) - 31*sqrt(2)*sqrt(239 + 169*sqrt(2)) + 960 + 679*sqrt(2)), r1*(-2085*sqrt(478 + 338*sqrt(2)) - 2944*sqrt(239 + 169*sqrt(2)) + 64416 + 45549*sqrt(2))/(-44*sqrt(239 + 169*sqrt(2)) - 31*sqrt(2)*sqrt(239 + 169*sqrt(2)) + 960 + 679*sqrt(2)), r1*(-sqrt(1912 + 1352*sqrt(2)) + 31 + 22*sqrt(2)))

res[1] |> println
res[1](r1 => 1/2).evalf() |> println

   4*r1*(-6*sqrt(239 + 169*sqrt(2)) + 4 + 3*sqrt(2) + 4*sqrt(478 + 338*sqrt(2)))
   1.48071257037594

長方形の一辺の長さは甲円の半径の 4(-6*sqrt(239 + 169√2) + 4 + 3√2 + 4sqrt(478 + 338√2)) 倍である。
これはこれ以上簡約化できない。

長方形の短辺が 1 寸のとき,長辺は 1.48071 寸である。

「術」では,1.481136081497563 倍としているが,数値は微妙に違うし,山村は 術の 2 を 2.21 に変更して(理論が不明であるし),「答」に合わせているが,それは本末転倒である。

天 = √2 + 3
(2.21 - √天)*天 + 1

   1.481136081497563

その他のパラメータは以下の通りである。

  r1 = 0.5;  a = 1.48071;  r2 = 0.240356;  r3 = 0.153333;  x3 = 0.426939;  r4 = 0.137032

function draw(r1, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r2, r3, x3, r4) = (4*r1*(-6*sqrt(239 + 169*sqrt(2)) + 4 + 3*sqrt(2) + 4*sqrt(478 + 338*sqrt(2))), r1*(-12*sqrt(239 + 169*sqrt(2)) + 7 + 6*sqrt(2) + 8*sqrt(478 + 338*sqrt(2))), 2*r1*(-145*sqrt(239 + 169*sqrt(2)) - 100*sqrt(478 + 338*sqrt(2)) + 3131 + 2214*sqrt(2))/(-44*sqrt(239 + 169*sqrt(2)) - 31*sqrt(2)*sqrt(239 + 169*sqrt(2)) + 960 + 679*sqrt(2)), r1*(-2085*sqrt(478 + 338*sqrt(2)) - 2944*sqrt(239 + 169*sqrt(2)) + 64416 + 45549*sqrt(2))/(-44*sqrt(239 + 169*sqrt(2)) - 31*sqrt(2)*sqrt(239 + 169*sqrt(2)) + 960 + 679*sqrt(2)), r1*(-sqrt(1912 + 1352*sqrt(2)) + 31 + 22*sqrt(2)))
   @printf("長方形の短辺が %g のとき,長辺は %g である。\n", 2r1, a)
   @printf("r1 = %g;  a = %g;  r2 = %g;  r3 = %g;  x3 = %g;  r4 = %g\n", r1, a, r2, r3, x3, r4 )
   plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   circle(a - r1, r1, r1)
   circle(r2, r1, r2, :green)
   circle(x3, r3, r3, :orange)
   circle(x3, 2r1 - r3, r3, :orange)
   circle(r4, r4, r4, :magenta)
   circle(r4, 2r1 - r4, r4, :magenta)
   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/2)
       point(a - r1, r1, "甲円:r1,(a-r1,r1)", :red, :center, delta=-delta/2)
       point(r2, r1, "乙円:r2,(r2,r1)", :green, :center, delta=-delta/2)
       point(x3, r3, "丙円:r3\n(x3,r3)", :orange, :center, delta=-delta/2)
       point(r4, r4, "丁円:r4\n(r4,r4)", :magenta, :center, delta=-delta/2)
       point(a, 2r1, "(a,2r1)", :black, :right, :bottom, delta=delta/2)
   end
end;

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

算額(その1122)

2024年07月05日 | Julia

算額(その1122)

二十九 一関市山ノ目 配志和神社 推定文政8年頃
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html
キーワード:円8個,外円,二等辺三角形

外円の中に二等辺三角形と甲円 4 個,乙円 1 個,丙円 2 個を容れる。丙円の直径が 1 寸のとき,乙円の直径はいかほどか。

二等辺三角形の底辺と y 軸の交点座標を (0, y)
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (r1, y + r1), (r1, y - r1)
乙円の半径と中心座標を r2, (0, y2)
丙円の半径と中心座標を r3, (x3, y3)
とおき,以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms R::positive, y::negative, r1::positive,
     r2::positive, y2::positive,
     r3::positive, x3::positive, y3::positive
eq1 = r1^2 + (y - r1)^2 - (R - r1)^2
eq2 = x3^2 + y3^2 - (R - r3)^2
eq3 = r1^2 + (y2 - y - r1)^2 - (r1 + r2)^2
eq4 = dist2(0, R, sqrt(R^2 - y^2), y, r1, y + r1, r1)
eq5 = dist2(0, R, sqrt(R^2 - y^2), y, 0, y2, r2)
eq6 = dist2(0, R, sqrt(R^2 - y^2), y, x3, y3, r3)
eq7 = y3/x3 - sqrt(R^2 - y^2)/(R - y);

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, r1, r2, y2, x3, y3) = u
   return [
       r1^2 - (R - r1)^2 + (-r1 + y)^2,  # eq1
       x3^2 + y3^2 - (R - r3)^2,  # eq2
       r1^2 - (r1 + r2)^2 + (-r1 - y + y2)^2,  # eq3
       (-R*(-R^2 + 2*R*r1 + R*y + 2*r1*sqrt(R^2 - y^2) + y^2)/2 + r1^2*sqrt(R^2 - y^2) + r1*y^2 + r1*y*sqrt(R^2 - y^2) + y^3/2)/R,  # eq4
       R^3 + R^2*y - 2*R^2*y2 - 2*R*r2^2 - 2*R*y*y2 + R*y2^2 + y*y2^2,  # eq5
       (R*(R^2 + R*y - 2*R*y3 - 2*r3^2 + x3^2 - 2*x3*sqrt(R^2 - y^2) - 2*y*y3 + y3^2) - x3^2*y + 2*x3*y3*sqrt(R^2 - y^2) + y*y3^2)/(2*R),  # eq6
       -sqrt(R^2 - y^2)/(R - y) + y3/x3,  # eq7
   ]
end;
r3 = 1/2
iniv = BigFloat[2.5, -0.7, 0.8, 11/16, 42/31, 1.6, 1.2]
res = nls(H, ini=iniv)

   ([2.5, -0.7, 0.8, 0.6875388202501893, 1.3541019662496845, 1.6, 1.2], true)

丙円の直径が 1 寸のとき,乙円の直径は 2*0.6875388202501893 = 1.3750776405003786 寸である。

function draw(R, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1/2
   (R, y, r1, r2, y2, x3, y3) = res[1]
   @printf("丙円の直径が %g のとき,乙円の直径は %g である。\n", 2r3, 2r2)
   x = sqrt(R^2 - y^2)
   plot([x, 0, -x, x], [y, R, y, y], color=:black, lw=0.5)
   circle(0, 0, R)
   circle2(r1, y + r1, r1, :green)
   circle2(r1, y - r1, r1, :green)
   circle(0, y2, r2, :magenta)
   circle2(x3, y3, r3, :orange)
   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, 0, " R", :magenta, :left, :bottom, delta=delta/2)
       point(r1, y - r1, "甲円:r1,(r1,y-r1)", :green, :center, delta=-delta/2)
       point(r1, y + r1, "甲円:r1,(r1,y+r1)", :green, :center, delta=-delta/2)
       point(0, y2, "乙円:r2,(0,y2)", :magenta, :center, delta=-delta/2)
       point(x3, y3, "丙円:r3\n(x3,y3)", :blue, :center, delta=-delta/2)
   end
end;

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

算額(その1120)

2024年07月05日 | Julia

算額(その1120)

十四 前沢町赤生津長根 前沢月山神社 明治11年(1878)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html
キーワード:球3個,四角柱,3次元

四角柱(底面が正方形)の中に等球を 3 個容れる。底面の正方形の一辺の長さが 20 寸,高さが 30 寸のとき,等球の直径を求めよ。

正方形の一辺の長さを a,高さを h
等球の半径と中心座標を下から順に r, (r, r, r), (a - r, a - r, z), (r, r, h - r); 上下対称なので z = h/2
とおき,以下の方程式を解く。

include("julia-source.txt")

using SymPy
@syms a::positive, h::positive, r::positive, z::positive
z = h/2
eq = 2(a - 2r)^2  + (z - r)^2 - 4r^2;

ans_r = solve(eq, r)[1];

ans_r(a => 20, h => 30) |> println
ans_r(a => 20, h => 30).evalf() |> println

   19 - 2*sqrt(39)
   6.51000400320320

底面の正方形の一辺の長さが 20 寸,高さが 30 寸のとき,等球の半径は 19 - 2√39 で,ほぼ 6.51 寸である。直径は 13.02 寸である。

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

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

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