裏 RjpWiki

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

算額(その186)

2023年04月06日 | Julia

算額(その186)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1(86)
長野県中野市田上 田上観音堂 文化6年(1809)

第3問 台形内を対角線で区切り,大円,小円を入れる。大円の径が36寸,小円の径が20寸,長が84寸のとき,大頭はいかほどか。

「長」を x,「大頭」を y とする。また,大円,小円の半径と中心座標を r1, (x1, y1), r2, (x2, y2) とする。
台形の右上の頂点座標を (x0, y0) とする。

x1 = r1 = 36/2, r2 = 20/2, x = x0 = 84, x2 = x - r2 のように記号を定め方程式を解く。
4 本の方程式はすべて,大円,小円と対角線の距離に関するものである。なお,それぞれの方程式を simplify() しないと解が求まらない。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す!!
end;

@syms x1::positive, y1::positive, r1::positive,
     x2::positive, y2::positive, r2::positive,
     x::positive, y::positive, x0::positive, y0::positive;

x1 = r1 = 36//2
r2 = 20//2
x0 = x = 84
x2 = x - r2
eq1 = distance(0, y, x, 0, x1, y1) - r1^2 |> simplify
eq2 = distance(0, 0, x0, y0, x1, y1) - r1^2 |> simplify
eq3 = distance(0, y, x, 0, x2, y2) - r2^2 |> simplify
eq4 = distance(0, 0, x0, y0, x2, y2) - r2^2 |> simplify

res = solve([eq1, eq2, eq3, eq4], (y, y0, y1, y2))

   3-element Vector{NTuple{4, Sym}}:
    (63, 35, 27, 20)
    (9*(-171*sqrt(1358) + sqrt(53615198))*sqrt(3*sqrt(39481) + 1025)/49664, 5*(85*sqrt(1358) + sqrt(53615198))*sqrt(3*sqrt(39481) + 1025)/86912, -387*sqrt(1358)*sqrt(3*sqrt(39481) + 1025)/86912 + 9*sqrt(53615198)*sqrt(3*sqrt(39481) + 1025)/86912, 5*sqrt(6*sqrt(39481)/679 + 2050/679))
    (3*(sqrt(1738846) + 281*sqrt(46))*sqrt(sqrt(37801) + 405)/4508, 5*(-sqrt(1738846) + 307*sqrt(46))*sqrt(sqrt(37801) + 405)/10304, -39*sqrt(46)*sqrt(sqrt(37801) + 405)/4508 + 3*sqrt(1738846)*sqrt(sqrt(37801) + 405)/4508, 5*sqrt(sqrt(37801)/46 + 405/46))

最初の組の解が適切解である。

using Printf
for i in 1:3
   @printf("y = %6.3f; y0 = %6.3f;  y1 = %6.3f; y2 = %6.3f\n", res[i][1].evalf(), res[i][2].evalf(), res[i][3].evalf(), res[i][4].evalf())
end

   y = 63.000; y0 = 35.000;  y1 = 27.000; y2 = 20.000
   y =  7.447; y0 = 24.216;  y1 = 23.922; y2 = 10.926
   y = 52.537; y0 =  9.071;  y1 = 20.048; y2 = 18.049

大頭の長さは y で,元の単位では 63 寸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=:black, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (y, y0, y1, y2) = res[1]
   x1 = r1 = 36//2
   r2 = 20//2
   x0 = x = 84
   x2 = x - r2
   @printf("y = %6.3f; y0 = %6.3f;  y1 = %6.3f; y2 = %6.3f\n", y, y0, y1, y2)
   plot([0, x, x, 0, 0], [0, 0, y0, y, 0], color=:black, lw=0.5)
   circle(x1, y1, r1, :red)
   circle(x2, y2, r2, :blue)
   segment(0, 0, x0, y0, :black, linewidth=0.25)
   segment(0, y, x, 0, :black, linewidth=0.25)
   if more == true
       point(x1, y1, "大円:(x1,y1,r1)", :green, :bottom)
       point(x2, y2, "小円:(x2,y2,r2)", :green, :bottom)
       point(x0, y0, " (x0,y0)", :green, :left, :bottom)
       point(0, y, " y", :green, :left, :bottom)
       point(x, 0, " x", :green, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

   y = 63.000; y0 = 35.000;  y1 = 27.000; y2 = 20.000

 

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

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

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