算額（その479）

算額（その479）

宮城県丸森町小斎日向 鹿島神社 明治9年

徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf

外円内に甲，乙，丙，丁円が入っている。丁円の直径を知って，丙円の直径を求めよ。

甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (0, r2), (0, r0 - r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (x4, r4)
とおいて以下の連立方程式の数値解を求める。

include("julia-source.txt")

using SymPy
@syms r1::positive, x1::positive, y1::positive, r2::positive,
     r3::positive, x3::positive, r4::positive,
     x4::positive, r0::positive;
r0 = 4r2
y1 = 2r2
r4 = 1//2
eq1 = x1^2 + y1^2 - (r0 - r1)^2
eq2 = x3^2 + r3^2 - (r0 - r3)^2
eq3 = x1^2 + (r0 - r2 - y1)^2 - (r1 + r2)^2
eq4 = (x3 - x1)^2 + (y1 - r3)^2 - (r1 + r3)^2
eq5 = (x1 - x4)^2 + (y1 - r4)^2 - (r1 + r4)^2
eq6 = x4^2 + (r2 - r4)^2 - (r2 + r4)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6])#, (r1, x1, r2, r3, x3, x4))

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)
   (r1, x1, r2, r3, x3, x4) = u
   return [
       4*r2^2 + x1^2 - (-r1 + 4*r2)^2,  # eq1
       r3^2 + x3^2 - (4*r2 - r3)^2,  # eq2
       r2^2 + x1^2 - (r1 + r2)^2,  # eq3
       -(r1 + r3)^2 + (2*r2 - r3)^2 + (-x1 + x3)^2,  # eq4
       -(r1 + r4)^2 + (2*r2 - r4)^2 + (x1 - x4)^2,  # eq5
       x4^2 + (r2 - r4)^2 - (r2 + r4)^2,  # eq6
   ]
end;

r0 = 4r2
y1 = 2r2
r4 = 1/2

iniv = [big"2.6", 4.3, 2.2, 1.5, 7.1, 3.0]
res = nls(H, ini=iniv)

  (BigFloat[1.311396103067892771960759925894364135356352343919626632929505882095829615307969, 2.141500868793756520305245638521174546550359518581729934628195842934091980001108, 1.092830085889910643300633271578636779463626953266355527441254901746524679423308, 0.7285533905932737622004221810524245196424179688442370182941699344976831196155384, 3.569168114656260867175409397535290910917265864302883224380326404890153300001858, 1.478397839480233171313044189429409031462889497069357846136076631538706602818749], true)

   r1 = 1.3114;  x1 = 2.1415;  r2 = 1.09283;  r3 = 0.728553;  x3 = 3.56917;  x4 = 1.4784
   丙円の直径は丁円の直径の 1.45711 倍

丙円の半径は丁円の半径の 1.4571067811865475 倍である。
この数値は術によれば sqrt(1/2)+3/4 である。

Float64(res[1][4]/r4) |> println
sqrt(1/2)+3/4 |> println

   1.4571067811865475
   1.4571067811865475

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x1, r2, r3, x3, x4) = res[1]
   r0 = 4r2
   y1 = 2r2
   r4 = 1/2
   @printf("r1 = %g;  x1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  x4 = %g\n", r1, x1, r2, r3, x3, x4)
   @printf("丙円の直径は丁円の直径の %g 倍\n", r3/r4)
   plot()
   circle(0, 0, 4r2, :black)
   circle4(x1, y1, r1, :magenta)
   circle(0, r0 - r2, r2, :blue)
   circle(0, r2 - r0, r2, :blue)
   circle(0, r2, r2, :blue)
   circle(0, -r2, r2, :blue)
   circle4(x3, r3, r3, :green)
   circle4(x4, r4, r4, :orange)
   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(x1, y1, "甲円:r1,(x1,y1)", :magenta, :center, :top, delta=-delta/2)
       point(0, r0 - r2, " 乙円:r2,(0,r0-r2)", :blue, :left, :vcenter)
       point(0, r2, " 乙円:r2,(0,r2)", :blue, :left, :vcenter)
       point(x3, r3, "丙円::r3,(x3,r3) ", :green, :right, :vcenter)
       point(x4, r4, " 丁円::r4,(x4,r4) ", :orange, :left, :top, delta=-delta/2)
       point(0, r0, " r0", :black, :left, :bottom, delta=delta/3)
   end
end;

 

算額（その478）

算額（その478）

宮城県丸森町小斎日向 鹿島神社 明治9年

徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf

鈎股弦があり，鈎，股，弦を直径としその中点を中心とする 3 個の円，天円，地円，人円がある。
天円の中に星円があり，星園は人円，地円と外接し，天円に内接する。
地円，人円の直径を知って星円の直径を求めよ。

地円の直径は股，人円の直径は鈎である。
天円の半径と中心座標を r1, (r2, r3)
地円の半径と中心座標を r2, (r2, 0)
人円の半径と中心座標を r3, (0, r3)
星円の半径と中心座標を r4, (x4, y4)
とおいて以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     r1::positive, r2::positive, r3::positive,
     r4::positive, x4::positive, y4::positive;

r1 = sqrt(鈎^2 + 股^2)/2
(r2, r3) = (股, 鈎).//2;

eq1 = (r2 - x4)^2 + (y4 - r3)^2 - (r1 - r4)^2
eq2 = (r2 - x4)^2 + y4^2 - (r2 + r4)^2
eq3 = x4^2 + (y4 - r3)^2 - (r3 + r4)^2
res = solve([eq1, eq2, eq3], (r4, x4, y4))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (股^2*鈎^2*(4*股^3 + 2*股^2*鈎 + 4*股^2*sqrt(股^2 + 鈎^2) + 3*股*鈎^2 + 2*股*鈎*sqrt(股^2 + 鈎^2) + 鈎^3 + 鈎^2*sqrt(股^2 + 鈎^2))/(8*股^6 + 8*股^5*sqrt(股^2 + 鈎^2) + 8*股^4*鈎^2 + 4*股^3*鈎^3 + 4*股^3*鈎^2*sqrt(股^2 + 鈎^2) + 5*股^2*鈎^4 + 4*股^2*鈎^3*sqrt(股^2 + 鈎^2) + 4*股*鈎^5 + 4*股*鈎^4*sqrt(股^2 + 鈎^2) + 2*鈎^6 + 2*鈎^5*sqrt(股^2 + 鈎^2)), 股*鈎^2*(4*股^4 + 6*股^3*鈎 + 4*股^3*sqrt(股^2 + 鈎^2) + 7*股^2*鈎^2 + 6*股^2*鈎*sqrt(股^2 + 鈎^2) + 5*股*鈎^3 + 5*股*鈎^2*sqrt(股^2 + 鈎^2) + 2*鈎^4 + 2*鈎^3*sqrt(股^2 + 鈎^2))/(8*股^6 + 8*股^5*sqrt(股^2 + 鈎^2) + 8*股^4*鈎^2 + 4*股^3*鈎^3 + 4*股^3*鈎^2*sqrt(股^2 + 鈎^2) + 5*股^2*鈎^4 + 4*股^2*鈎^3*sqrt(股^2 + 鈎^2) + 4*股*鈎^5 + 4*股*鈎^4*sqrt(股^2 + 鈎^2) + 2*鈎^6 + 2*鈎^5*sqrt(股^2 + 鈎^2)), 股^2*鈎*(4*股^3 + 2*股^2*鈎 + 4*股^2*sqrt(股^2 + 鈎^2) + 3*股*鈎^2 + 2*股*鈎*sqrt(股^2 + 鈎^2) + 鈎^3 + 鈎^2*sqrt(股^2 + 鈎^2))/(4*股^5 + 4*股^4*sqrt(股^2 + 鈎^2) + 3*股^3*鈎^2 + 2*股^2*鈎^3 + 股^2*鈎^2*sqrt(股^2 + 鈎^2) + 2*股*鈎^4 + 2*股*鈎^3*sqrt(股^2 + 鈎^2) + 2*鈎^5 + 2*鈎^4*sqrt(股^2 + 鈎^2)))

星円の半径のみならず中心座標も大変長い式になる。
鈎 = 3，股 = 4 のとき，星円の直径は 2 になる。

2res[1][1](鈎 => 3, 股 => 4) |> println

   2

   星円の直径 = 2;  r1 = 2.5;  r2 = 2;  r3 = 1.5;  r4 = 1;  x4 = 2;  y4 = 3

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股) = (3, 4)
   r1 = sqrt(鈎^2 + 股^2)/2
   (r2, r3) = (股, 鈎).//2
   (r4, x4, y4) = (股^2*鈎^2*(4*股^3 + 2*股^2*鈎 + 4*股^2*sqrt(股^2 + 鈎^2) + 3*股*鈎^2 + 2*股*鈎*sqrt(股^2 + 鈎^2) + 鈎^3 + 鈎^2*sqrt(股^2 + 鈎^2))/(8*股^6 + 8*股^5*sqrt(股^2 + 鈎^2) + 8*股^4*鈎^2 + 4*股^3*鈎^3 + 4*股^3*鈎^2*sqrt(股^2 + 鈎^2) + 5*股^2*鈎^4 + 4*股^2*鈎^3*sqrt(股^2 + 鈎^2) + 4*股*鈎^5 + 4*股*鈎^4*sqrt(股^2 + 鈎^2) + 2*鈎^6 + 2*鈎^5*sqrt(股^2 + 鈎^2)), 股*鈎^2*(4*股^4 + 6*股^3*鈎 + 4*股^3*sqrt(股^2 + 鈎^2) + 7*股^2*鈎^2 + 6*股^2*鈎*sqrt(股^2 + 鈎^2) + 5*股*鈎^3 + 5*股*鈎^2*sqrt(股^2 + 鈎^2) + 2*鈎^4 + 2*鈎^3*sqrt(股^2 + 鈎^2))/(8*股^6 + 8*股^5*sqrt(股^2 + 鈎^2) + 8*股^4*鈎^2 + 4*股^3*鈎^3 + 4*股^3*鈎^2*sqrt(股^2 + 鈎^2) + 5*股^2*鈎^4 + 4*股^2*鈎^3*sqrt(股^2 + 鈎^2) + 4*股*鈎^5 + 4*股*鈎^4*sqrt(股^2 + 鈎^2) + 2*鈎^6 + 2*鈎^5*sqrt(股^2 + 鈎^2)), 股^2*鈎*(4*股^3 + 2*股^2*鈎 + 4*股^2*sqrt(股^2 + 鈎^2) + 3*股*鈎^2 + 2*股*鈎*sqrt(股^2 + 鈎^2) + 鈎^3 + 鈎^2*sqrt(股^2 + 鈎^2))/(4*股^5 + 4*股^4*sqrt(股^2 + 鈎^2) + 3*股^3*鈎^2 + 2*股^2*鈎^3 + 股^2*鈎^2*sqrt(股^2 + 鈎^2) + 2*股*鈎^4 + 2*股*鈎^3*sqrt(股^2 + 鈎^2) + 2*鈎^5 + 2*鈎^4*sqrt(股^2 + 鈎^2)))
   @printf("星円の直径 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  r4 = %g;  x4 = %g;  y4 = %g\n", 2r4, r1, r2, r3, r4, x4, y4)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   circle(r2, r3, r1)
   circle(r2, 0, r2, :magenta)
   circle(0, r3, r3, :green)
   circle(x4, y4, r4, :orange)
   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(r2, r3, " 天円:r1,(r2,r3)", :red, :left, :top, delta=-delta)
       point(r2, 0, " 地円:r2,(r2,0)", :magenta, :left, :top, delta=-delta)
       point(0, r3, " 人円:r3,(0,r3)", :green, :left, :top, delta=-delta)
       point(x4, y4, " 星円:r4,(x4,y4)", :orange, :left, :top, delta=-delta)
   end
end;

算額（その477）

算額（その477）

宮城県丸森町小斎日向 鹿島神社 明治9年

徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf

円内に等円 2 個，菱形 2 個が入っている。等円の直径が与えられたとき，菱形の一辺の長さを求めよ。

菱形の短い方の対角線の長さを 2a，長い方の対角線の長さを 2b とする。菱形の右の頂点の座標は (b, 2r - a) である。頂点は円周上にあるので，b = sqrt((2r)^2 - (2r - a)^2) である。

以下の方程式で a を求める。

include("julia-source.txt")

using SymPy
@syms a::positive, b::positive, r::positive;

b = sqrt((2r)^2 - (2r - a)^2)
eq = distance(0, 2r - 2a, b, 2r - a, r, 0) - r^2
res = solve(eq, a);

6 個の解が得られるが，2 番目のものが適解である。

すなわち，a = r*(2 - √7/2) である。

res[2] |> println

   r*(2 - sqrt(7)/2)

菱形の一辺の長さは sqrt(a^2 + b^2) = r*(√7 - 1) である。

@syms a::positive, b::positive, r::positive;
a = r*(2 - sqrt(Sym(7))/2)
b = sqrt((2r)^2 - (2r - a)^2)
sqrt(a^2 + b^2) |> simplify |> sympy.sqrtdenest |> println

   r*(-1 + sqrt(7))

等円の直径が 1 寸のとき，菱形の一辺の長さは 0.5(√7 - 1) = 0.8228756555322954 = 8分2厘2毛8糸有奇である。

0.5(√7 - 1) |> println

   0.8228756555322954

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1//2
   a = r*(2 - sqrt(7)/2)
   b = sqrt((2r)^2 - (2r - a)^2)
   length = r*(sqrt(7) - 1)
   @printf("a = %g; b = %g;  菱形の一辺の長さ = %g\n", a, b, length)
   plot([0, b,  0, -b, 0], [2r - 2a, 2r - a, 2r, 2r - a, 2r - 2a], color=:brown, lw=0.5)
   plot!([0, b,  0, -b, 0], [2a - 2r, a - 2r, -2r, a - 2r, 2a - 2r], color=:brown, lw=0.5)
   circle(0, 0, 2r, :blue)
   circle(r, 0, r)
   circle(-r, 0, 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(0, 2r, " 2r", :blue, :left, :bottom, delta=delta/2)
       point(0, 2r - a, " 2r-a", :brown, :left, :bottom, delta=delta/2)
       point(b, 2r - a, "(b,2r-a)", :brown, :left, :bottom, delta=delta/2)
       point(0, 2r - 2a, " 2r-2a", :black, :left, :top)
       point(r, 0, "等円:r,(r,0)", :red, :center, :bottom, delta=delta)
   end
end;