裏 RjpWiki

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

算額(その613)

2024年01月05日 | Julia

算額(その613)

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション

https://dl.ndl.go.jp/pid/3508431/1/6

円弧内に斜線と等円 2 個を入れる。等円の直径が 4 寸で,弦の長さから「子 = 14 寸」を引いた丑の長さ」(図での xa - x)を求めよ。

円弧をなす円の半径と中心座標を R, (0, 0)
弦と y 軸の交点座標を (0, y) とする。
等円の半径と中心座標を r, (x, y + r), (x2, y2); x = 14 - sqrt(R^2 - y^2)

include("julia-source.txt");

using SymPy
@syms r::positive, R::positive, x::positive, y::positive,
     x2::negative, y2::positive, y3::positive, xa::positive, xb::positive
r = 4//2
xa = sqrt(R^2 - y^2)
xb = sqrt(R^2 - y3^2)
eq1 = (14 - sqrt(R^2 - y^2))^2 + (y + r)^2 - (R - r)^2
eq2 = x2^2 + y2^2 - (R - r)^2
eq3 = (y3 - y)/(xa + xb) * y2/x2 + 1
eq4 = distance(-xa, y, xb, y3, 14 - sqrt(R^2 - y^2), y + r) - r^2
eq5 = (xb - xa)^2/4 + ((y3 - y)/2+y)^2  - (R - 2r)^2;

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 v, r.f_converged
end;

function H(u)
   (R, y, x2, y2, y3) = u
   return [
       (14 - sqrt(R^2 - y^2))^2 - (R - 2)^2 + (y + 2)^2,  # eq1
       x2^2 + y2^2 - (R - 2)^2,  # eq2
       1 + y2*(-y + y3)/(x2*(sqrt(R^2 - y^2) + sqrt(R^2 - y3^2))),  # eq3
       (y - ((y + 2)*(R^2 - y*y3 + sqrt(R^2 - y^2)*sqrt(R^2 - y3^2)) - (sqrt(R^2 - y^2) + sqrt(R^2 - y3^2))*(7*y - 7*y3 + sqrt(R^2 - y^2) + sqrt(R^2 - y3^2)))/(R^2 - y*y3 + sqrt(R^2 - y^2)*sqrt(R^2 - y3^2)) + 2)^2 + (-sqrt(R^2 - y^2) - ((14 - sqrt(R^2 - y^2))*(R^2 - y*y3 + sqrt(R^2 - y^2)*sqrt(R^2 - y3^2)) - (y - y3)*(7*y - 7*y3 + sqrt(R^2 - y^2) + sqrt(R^2 - y3^2)))/(R^2 - y*y3 + sqrt(R^2 - y^2)*sqrt(R^2 - y3^2)) + 14)^2 - 4,  # eq4
       -(R - 4)^2 + (y/2 + y3/2)^2 + (-sqrt(R^2 - y^2) + sqrt(R^2 - y3^2))^2/4,  # eq5
   ]
end;

iniv = BigFloat[8.8, 2.5, -1.67, 6.6, 6.7]
res = nls(H, ini=iniv)

   (BigFloat[9.165816326530612244897959183673469387755102040816326530612244897959183673469382, 2.839183673469387755102040816326530612244897959183673469387755102040816326530621, -2.00642857142857142857142857142857142857142857142857137843826725205845906344671, 6.879183673469387755102040816326530612244897959183673492933533457683883923317655, 7.079183673469387755102040816326530612244897959183673469387755102040816326530546], true)

丑の長さは xa - x = 8.715 - 5.285 = 3.43 寸である。

   R = 9.16582;  x = 5.285;  y = 2.83918;  x2 = -2.00643;  y2 = 6.87918;  y3 = 7.07918;  xa = 8.715;  xb = 5.82214

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 2
   (R, y, x2, y2, y3) = res[1]
   xa = sqrt(R^2 - y^2)
   xb = sqrt(R^2 - y3^2)
   x = 14 - xa
   @printf("R = %g;  x = %g;  y = %g;  x2 = %g;  y2 = %g;  y3 = %g;  xa = %g;  xb = %g\n", R, x, y, x2, y2, y3, xa, xb)
   @printf("丑 = %g\n", xa - x)
   plot()
   circle(0, 0, R, beginangle=0, endangle=180)
   circle(x, y + r, r)
   circle(x2, y2, r)
   plot!([sqrt(R^2 - y3^2), -sqrt(R^2 - y^2), sqrt(R^2 - y^2)], [y3, y, y], color=:blue, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, y, " y", :blue, :left, delta=-delta/2)
       point(-xa, y, "(-xa,y)", :blue)
       point(xa, y, "(xa,y)", :blue, :right, delta=-delta/2)
       point(x, y, "(x,y)", :blue, :center, delta=-delta/2)
       point(sqrt(R^2 - y3^2), y3, " (xb,y3)", :blue, :left, :vcenter)
       point(x, y + r, "(x,y+r)", :red, :center, delta=-delta/2)
       point(x2, y2, "(x2,y2)", :red, :center, delta=-delta/2)
       point(R, 0, "R ", :red, :right, :bottom, delta=delta/2)
       segment(x, y, xa, y, lw=2)
       point((x+xa)/2, y, "丑", :black, :center, :bottom, delta=delta, mark=false)
   end
end;

 

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

算額(その612)

2024年01月05日 | Julia

算額(その612)

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション

https://dl.ndl.go.jp/pid/3508431/1/9

鈎股弦(直角三角形)内に長方形と中鈎(直角の頂点とから斜辺におろした垂線)と 2 個の等円を入れる。長方形の長辺と短辺が 40 寸,20 寸であるとき,等円の直径を求めよ。

長方形の長辺と短辺を「長」,「平」
直角三角形の直角を挟む二辺を「鈎」,「股」
中鈎と斜辺の交点座標を (x2, y2)
中鉤と長方形の長辺の交点座標を (x, 20)
等円の半径と中心座標を r, (r, a), (b, 20 + r)
とおき,以下の連立方程式を解き,r, x, 鈎, 股, x2, y2 を求める。
なお,等円の半径 r を求めるだけであれば,eq1, eq2 のみの二元連立方程式を解くだけでよい。

include("julia-source.txt");

using SymPy
@syms 長::positive, 平::positive,
     鈎::positive, 股::positive, x::positive, y::positive, r::positive,
     x2::positive, y2::positive, a::positive, b::positive
eq1 = sqrt(x^2 + 平^2) - (長 - x)
eq2 = 平 + x - sqrt(x^2 + 平^2) - 2r
eq3 = x/平 - 鈎/股
eq4 = 長*平 + 長*(鈎 - 平)/2 + 平*(股 - 長)/2 - 鈎*股/2
eq5 = x2^2 + y2^2 + (鈎 -y2)^2 + x2^2 - 鈎^2
eq6 = x2^2 + y2^2 + (股 - x2)^2 + y2^2 - 股^2;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (r, x, 鈎, 股, x2, y2))

   1-element Vector{NTuple{6, Sym}}:
    ((-平^2 + 長*(2*平 + 長) - sqrt(4*平^2*長^2 + (平^2 - 長^2)^2))/(4*長), (-平^2 + 長^2)/(2*長), (平^2 + 長^2)/(2*平), -長*(平^2 + 長^2)/(平^2 - 長^2), 長*(-平^2 + 長^2)/(平^2 + 長^2), 2*平*長^2/(平^2 + 長^2))

長,平が 40, 20 のとき,等円の半径は 5 である。

半径 = (-平^2 + 長*(2*平 + 長) - sqrt(4*平^2*長^2 + (平^2 - 長^2)^2))/(4*長)
半径(長 => 40, 平 => 20)

   5

半径の式中の sqrt の中は平方になるので sqrt を剥ぎ取ることができる。

sqrt(4*平^2*長^2 + (平^2 - 長^2)^2) |> factor |> println

   平^2 + 長^2

2倍して直径を求める式を得る。

2(-平^2 + 長*(2*平 + 長) - (平^2 + 長^2))/(4*長) |> simplify |> println

   平*(-平 + 長)/長

カッコの中を「長」で割ると,平*(1 - 平/長) となり,術の「平以長除之得數以減一個余剰平得等径合問」に一致する。
長,平が 40, 20 のとき,等円の直径は 10 である。

直径 = 平*(1 - 平/長)
直径(長 => 40, 平 => 20)

   10

「平」と「長」を未知数のままにしておいても構わないが,煩雑なので数値解を求めると,
(r, x, 鈎, 股, x2, y2) = (5, 15, 50, 200/3, 24, 32)
となる。このパラメータを既知として,更に図を描くために,等円の中心座標 (r, a), (b, 20 + r) の a, b を求める。

@syms 鈎::positive, 股::positive, x::positive, y::positive, r::positive,
     x2::positive, y2::positive, a::positive, b::positive
(r, x, 鈎, 股, x2, y2) = (5, 15, 50, 200/3, 24, 32)
eq7 = distance(0, 0, x2, y2, r, a) - r^2
eq8 = distance(0, 0, x2, y2, b, 20 + r) - r^2
solve([eq7, eq8], (a, b))

   2-element Vector{Tuple{Sym, Sym}}:
    (15, 25/2)
    (15, 25)

2 番目のものが適解である。すなわち等円の中心座標は (5, 15), (25, 25) である

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x, 鈎, 股, x2, y2) = (5, 15, 50, 200/3, 24, 32)
   (a, b) = (15, 25)
   @printf("等円の直径 = %g;  r = %g;  x = %g;  鈎 = %g;  股 = %g;  x2 = %g;  y2 = %g;  a = %g;  b = %g\n",
       2r, r, x, 鈎, 股, x2, y2, a, b)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   rect(0, 0, 40, 20, :blue)
   segment(0, 0, x2, y2, :green)
   circle(r, a, r)
   circle(b, 20 + r, r)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x, 20, "(x,20)", :black)
       point(x2, y2, " (x2,y2)", :black, :left, :vcenter)
       point(股, 0, "股", :black, :left, :bottom, delta=delta/2)
       point(0, 鈎, " 鈎", :black, :left, :bottom, delta=delta/2)
       point(r, a, "(r,a)", :red, :center, delta=-delta/2)
       point(b, 20 + r, "(b,20+r)", :red, :center, delta=-delta/2)
   end
end;

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

万博の縮小延期ではなく,中止せよ

2024年01月05日 | 雑感

吉村知事、万博の縮小延期を否定 能登地震との「二者択一でない」 2024/01/04 15:56共同通信https://news.goo.ne.jp/article/kyodo_nor/politics/kyodo_nor-2024010401000831.html

吉村知事の発言に対して、万博の縮小延期を否定する姿勢に疑問が残る。

まず、吉村知事が「能登半島地震の影響で縮小や延期をする可能性を否定した」と発言しているが、この主張には十分な説明や根拠が欠如している。能登半島地震がもたらす影響を具体的に検証し、それに基づいてなぜ縮小や延期が不要であると結論づけているのかについて、より詳細な情報が求められる。

また、「二者択一の関係ではない。万博があるから(復興の)費用が削減されるものではない」との発言も疑問を呼び起こす。吉村知事はなぜ、万博の開催と災害復興の費用には直接の因果関係がないと断言できるのか、その理由や具体的な根拠についての説明が必要だ。災害復興に使われるべき財源が万博に流用される可能性はないのか、十分な検証が求められる。

さらに、吉村知事は「被災地の救助活動や復興支援については国と自治体が全力でやるべきである」と主張しているが、それが万博の開催にどのように影響するかについてもっと具体的に説明されるべきだ。万博の成功に向けた勝負の年とされる中、ほんとうにそれほどの期待が寄せられているのか、その背景や具体的な計画についての情報が不足している。

総じて、吉村知事の発言には論理的な説明や具体的な根拠が欠如しており、災害と万博の関連性やリスクについて十分な説明が必要だ。

猛省と方針変換を望む。

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

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

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