算額（その305）

算額（その305）

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 菅原神社 嘉永7年(1854)

問題文4　

外大円の中に長方形を起き，菱形を切り取った残りの領域に円を置く。円の直径が 3 尺，長方形の残りの長さ BE が 4.5 尺である。このとき，菱形の一辺（菱面=BD），長軸（同立=AD），短軸（同横= BC），外大円の円周，矢（r0 - y）を求めよ。

外大円の半径を r0 とおく。r0 = 同立/2
長方形の右上の頂点の座標を (x, y) とする。
以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x::positive, y::positive,
   菱面::positive, 同立::positive, 同横::positive,
   外大円円周::positive, 矢::positive;

r1 = 3/2
eq1 = 4.5 + 2y - sqrt(4.5^2 + (2y)^2) - 2r1
eq2 = 4.5^2 + (2y)^2 - 菱面^2  # ⊿BDE において
eq3 = (菱面 + 4.5)^2 + (2y)^2 - 同立^2  # ⊿ADE において
eq4 = (同立/2)^2 + (同横/2)^2 - 菱面^2  # ⊿OBE において
eq5 = 外大円円周 - 同立*PI
eq6 = 矢 - (同立/2 - y);

eq1 〜 eq6 までの6元連立方程式にすると解けなくなってしまうので，eq1 〜 eq5 の5元連立方程式を解く。

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (y, 菱面, 同立, 同横, 外大円円周, 矢))

res = solve([eq1, eq2, eq3, eq4, eq5], (y, 菱面, 同立, 同横, 外大円円周))

   1-element Vector{NTuple{5, Sym}}:
    (3.00000000000000, 7.50000000000000, 13.4164078649987, 6.70820393249937, 42.1488883862444)

res[1][3]/2 - res[1][1] |> println  # 矢

   3.70820393249937

なお，solve(eq1) で y だけを求めれば，ドミノ倒しで残りの変数は確定する。

res = solve(eq1)[1] |> println  # y

   3.00000000000000

sqrt(6^2 + 4.5^2)  # 菱面 = sqrt((2y)^2 + 4.5^2)

   7.5

sqrt((7.5 + 4.5)^2 + (2*3)^2)  # 同立 = sqrt((菱面 + 4.5) + (2y)^2)

   13.416407864998739

sqrt(7.5^2 - (13.416407864998739/2)^2)*2  # 同横 = sqrt(菱面^2-(同立/2)^2)*2

   6.708203932499367

13.4164078649987*pi  # 外大円円周 = 同立*π

   42.14888838624424

(13.416407864998739/2) - 3  # 矢 = 同立/2 - y

   3.7082039324993694

   y = 3.000000
   菱面 = 7.500000
   同立 = 13.416408
   同横 = = 6.708204
   外大円円周 = = 42.148888
   矢 = 3.708204

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (y, 菱面, 同立, 同横, 外大円円周) = (3.00000000000000, 7.50000000000000, 13.4164078649987, 6.70820393249937, 42.1488883862444)
   矢 = 3.70820393249937
   @printf("y = %.6f\n", y)
   @printf("菱面 = %.6f\n", 菱面)
   @printf("同立 = %.6f\n", 同立)
   @printf("同横 = = %.6f\n", 同横)
   @printf("外大円円周 = = %.6f\n", 外大円円周)
   @printf("矢 = %.6f\n", 矢)
   plot([6, 6, -6, -6, 6], [-y, y, y, -y, -y], color=:black, lw=0.5)
   plot!([6, 1.5, -6, -1.5, 6], [-y, y, y, -y, -y], color=:red, lw=0.5)
   circle(4.5, 1.5, 1.5, :green)
   segment(1.5, y, -1.5, -y, :blue)
   segment(-6, y, 6, -y, :magenta)
   if more
       point(4.5, 1.5, "(4.5,1.5)", :green, :center)
       point(1.5, y, "B:(1.5,y) ", :red, :left, :bottom)
       point(6, y, "E:(6,y)", :black, :left, :bottom)
       point(-6, y, " A", :green, :left, :bottom)
       point(-1.5, -y, " C", :green, :left, :bottom)
       point(6, -y, " D", :green, :left, :bottom)
       point(0, 0, " O", :green, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;