裏 RjpWiki

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

算額(その101)

2023年01月18日 | Julia

算額(その101)

岩手県盛岡市の天満宮 文化 4 年(1807)
一関市博物館 平成29年度,中級問題
https://www.city.ichinoseki.iwate.jp/museum/wasan/h29/images/normal.pdf

正方形の中に半円 2 個,天円 1 個,地円 2 個がある。正方形の 1 辺が 8 寸のとき,天円径,地円径はいくらか。

正方形の一辺の長さの半分(半円の半径)を r として,下図のように記号を定め,方程式を解く。
x2 = r/2, x3 = r - x1 である。

using SymPy

@syms r::positive, x1::positive, x2::positive, x3::positive, r1::positive, r2::positive;
r = 8 // 2
x2 = r / 2
x3 = r - x1
eq1 = (r - x1)^2 + x1^2 - (r - r1)^2
eq2 = (r - x2)^2 + x2^2 - (r - r2)^2
eq3 = 2(x2 - x1)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (x1, r1, r2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-4*sqrt(2)*sqrt(11 - 6*sqrt(2))/7 + 2, -4/7 + 6*sqrt(2)/7, 4 - 2*sqrt(2))
    (4*sqrt(2)*sqrt(11 - 6*sqrt(2))/7 + 2, -4/7 + 6*sqrt(2)/7, 4 - 2*sqrt(2))

地円の径(直径)は (12√2 - 8) / 7 = 1.2815089640681632
天円の径(直径)は 8 - 4√2 = 2.3431457505076194

(12√2 - 8) / 7, 8 - 4√2

   (1.2815089640681632, 2.3431457505076194)

地円の中心座標 (x1, x1) は -4*sqrt(2)*(3 - sqrt(2))/7 + 2 から二重根号を外して (22 - 12√2) / 7 = 0.7184910359318367

sympy.sqrtdenest(-4*sqrt(Sym(2))*(3 - sqrt(Sym(2)))/7 + 2) |> println

   22/7 - 12*sqrt(2)/7

(22 - 12√2) / 7

   0.7184910359318367

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=: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=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, r1, r2) = 4 .* (-sqrt(2)*sqrt(11 - 6*sqrt(2))/7 + 1/2, -1/7 + 3*sqrt(2)/14, 1 - sqrt(2)/2)
   (x1, r1, r2) = (-4*sqrt(2)*sqrt(11 - 6*sqrt(2))/7 + 2, -4/7 + 6*sqrt(2)/7, 4 - 2*sqrt(2))
   r = 8 // 2
   x2 = r / 2
   x3 = r - x1
    plot(r .*[0, 2, 2, 0, 0], r .*[0, 0, 2, 2, 0], color=:black, lw=0.5)
   circle(r, 0, r, :green, beginangle=0, endangle=180)
   circle(0, r, r, :green, beginangle=270, endangle=450)
   circle(x1, x1, r1)
   circle(x3, x3, r1)
   circle(x2, x2, r2, :blue)
    if more
       info = @sprintf("x1 = %.3f, x2 = %.3f, x3 = %.3f, r1 = %.3f, r2 = %.3f\n", x1, x2, x3, r1, r2) *
              @sprintf("地円の直径 = %.3f,  天円の直径 = %.3f\n", 2r1, 2r2)
       point(0.1, 1.2r, info, :black, :left, mark=false)
       point(x1, x1, " 地円 (x1,x1); 半径:r1", :red, :left)
       point(x2, x2, " 天円 (x2,x2); 半径:r2", :green, :left)
    end
end;

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

算額(その100)

2023年01月18日 | Julia

算額(その100)

東京都杉並区和田 立法寺 天明8年(1788)正月
http://www.wasan.jp/tokyo/ryuhoji.html

直角三角形の中に甲,乙 2 つの円と,正方形が入っている。正方形の一辺の長さと乙円の直径が等しい。

乙円の半径を x として,左右反転させて,下図のように記号を定め,方程式を解く。

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))
end;

@syms x::positive, x2::positive, r::positive, x0::positive, y0::positive;
x = 1
eq1 = (x2 - x)^2 + (r - x)^2 - (r + x)^2
eq2 = (2x - x2)^2 + (2x - r)^2 - r^2
eq3 = x0 * (y0 - 4x) - 2x * y0
eq4 = distance(0, y0, x0, 0, x2, r) - r
res = solve([eq1, eq2, eq3, eq4], (r, x2, x0, y0))

   1-element Vector{NTuple{4, Sym}}:
    (25/16, 7/2, 25*sqrt(23)/28 + 173/28, -1676/49 + 400*sqrt(23)/49)

(25/16, 7/2, 25*sqrt(23)/28 + 173/28, -1676/49 + 400*sqrt(23)/49)

   (1.5625, 3.5, 10.460563860100642, 4.945563455614028)

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=: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=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x = 1
   (r , x2, x0, y0) = (25/16, 7/2, 25*sqrt(23)/28 + 173/28, -1676/49 + 400*sqrt(23)/49)
    plot([0, x0, 0, 0], [0, 0, y0, 0], color=:blue, lw=0.5)
    circle(x, x, x, :green)
   plot!([2x, 0, 0, 2x, 2x], [2x, 2x, 4x, 4x, 2x], color=:magenta, lw=0.5)
   circle(x2, r, r)
    if more
       info = @sprintf("x = %.0f, r = %.4f, x2 = %.1f, x0 = %.3f, y0 = %.3f\n", x, r, x2, x0, y0)
       point(1.7, 4.5, info, :black, :left, mark=false)
       point(x2, r, "甲円 (x2,r)", :red, :center)
       point(2x, 2x, "(2x,2x)", :magenta)
       point(x, x, "乙円 (x,x)", :green, :center)
       point(0, y0, " y0", :black, :left, :bottom)
       point(x0, 0, " x0", :black, :left, :bottom)
    end
end;

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

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

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