裏 RjpWiki

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

算額(その208)

2023年04月19日 | Julia

算額(その208)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2(113)
長野県長野市信更町田野口 清水神社 文政11年(1828)

鉤が 60 寸,股が 80 寸の鉤股弦内に径が 40 寸の全円を入れる。全円と中斜と股に接する甲円,全円と中斜と鉤に接する乙円の径を求めよ。

甲円,乙円の半径を r1, r2 とする。中心座標は (x1, r1), (r2, y2) とする。
中斜と弦は点 A で垂直に交わる。A の座標は (48*3/5, 48*4/5) である。
甲円,乙円が全円と外接する条件式と甲円,乙円の中心から中斜までの距離に関する条件式の連立方程式を解く。

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 r1::positive, x1::positive, r2::positive, y2::positive;

eq1 = (20 - x1)^2 + (20 - r1)^2 - (20 + r1)^2
eq2 = (20 - r2)^2 + (20 - y2)^2 - (20 + r2)^2
eq3 = distance(0, 0, 48*3/5, 48*4/5, x1, r1) - r1^2
eq4 = distance(0, 0, 48*3/5, 48*4/5, r2, y2) - r2^2
res = solve([eq1, eq2, eq3, eq4], (r1, x1, r2, y2))

   4-element Vector{NTuple{4, Sym}}:
    ((20 - 20*sqrt(3))^2/80, 40 - 20*sqrt(3), 20/9, 20/3)
    ((20 - 20*sqrt(3))^2/80, 40 - 20*sqrt(3), 20, 60)
    ((20 + 20*sqrt(3))^2/80, 20*sqrt(3) + 40, 20/9, 20/3)
    ((20 + 20*sqrt(3))^2/80, 20*sqrt(3) + 40, 20, 60)

using Printf
for i in 1:4
   @printf("r1 = %.5f, x1 = %.5f, r2 = %.5f, y2 = %.5f\n", res[i][1], res[i][2], res[i][3], res[i][4])
end

   r1 = 2.67949, x1 = 5.35898, r2 = 2.22222, y2 = 6.66667
   r1 = 2.67949, x1 = 5.35898, r2 = 20.00000, y2 = 60.00000
   r1 = 37.32051, x1 = 74.64102, r2 = 2.22222, y2 = 6.66667
   r1 = 37.32051, x1 = 74.64102, r2 = 20.00000, y2 = 60.00000

最初の組が適切解である。「甲円,乙円の直径は 5寸3分5厘あまりあり,4寸4分4厘あまりあり

(20 - 20*sqrt(3))^2/80*2, 20/9*2

   (5.358983848622452, 4.444444444444445)

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:blue, 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 draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x1, r2, y2) = res[1]
   @printf("r1 = %.5f, x1 = %.5f, r2 = %.5f, y2 = %.5f\n", r1, x1, r2, y2)
   plot([0, 80, 0, 0, 48*3/5], [0, 0, 60, 0, 48*4/5], color=1, lw=0.5)
   circle(20, 20, 20)
   circle(x1, r1, r1, :blue)
   circle(r2, y2, r2, :green)
   if more == true
       point(20, 20, "(20,20) 全円", :red)
       point(x1, r1, "  (x1,r1) 甲円", :blue)
       point(r2, y2, "    (r2,y2) 乙円", :green, :left, :bottom)
       point(48*3/5, 48*4/5, "  A (48*3/5, 48*4/5)", :orange, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

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

算額(その207)

2023年04月19日 | Julia

算額(その207)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2(113)
長野県長野市信更町田野口 清水神社 文政11年(1828)

直線の上に甲円,2 個の正方形とその上に乙円がある。正方形の一辺の長さと乙円の直径は 10 寸である。甲円の直径を求めよ。

甲円の半径を r,正方形の位置(図参照) x として連立方程式を解く。

using SymPy

@syms x::positive, r::positive;

eq1 = x^2 + (r -10)^2 - r^2
eq2 = (x + 10)^2 + (r - 15)^2 - (r + 5)^2;
res = solve([eq1, eq2], (r, x))

   1-element Vector{Tuple{Sym, Sym}}:
    (10*sqrt(2) + 20, 10 + 10*sqrt(2))

using Printf
@printf("r = %.5f;  x = %.5f\n", 10*sqrt(2) + 20, 10 + 10*sqrt(2))

   r = 34.14214;  x = 24.14214

甲円の直径は 68.28428 寸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:blue, 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 draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x) = res[1]
   @printf("r = %.5f;  x = %.5f\n", r, x)
   plot([x+10, x+10, x, x, x+20, x+20, x+10], [0, 10, 10, 0, 0, 10, 10], color=:green, lw=0.5)
   circle(0, r, r)
   circle(x + 10, 15, 5, :blue)
   hline!([0], color=:black, lw=0.5)
   if more == true
       point(0, r, " r", :red)
       point(15, 40, "甲円", :red, mark=false)
       point(40, 20, "乙円", :blue, mark=false) 
       point(x + 2, 6, "正方形", :green, mark=false)
       point(x + 12, 6, "正方形", :green, mark=false)
       point(x, 0, " x", :green, :left, :bottom)
       point(x + 10, 15, " (x+10,15)", :blue)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

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

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