算額（その320）

算額（その320）
解析解は 算額（その834）を参照のこと
長野県長野市戸隠 戸隠山中院権現堂 安永5年(1776)
中村信弥「改訂増補 長野県の算額」(26)
http://www.wasan.jp/zoho/zoho.html

外円の中に中円1個，小円2個が入っている。外円の面積から中円・小円の面積を引くと 120 歩である。また，中円と小円の直径の差は 5 寸である。小円の直径を求めよ。

 

以下の連立方程式を解く。
eq4 の 3 は，この算額で用いる円周率の近似値である(120 はこれに依存している）。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive, y::negative;
@syms R, r1, r2, y;

eq1 = r1 - r2 - 5//2
eq2 = r2^2 + (R - r1 - y)^2 - (r1 + r2)^2
eq3 = r2^2 + y^2 - (R - r2)^2
eq4 = 3*(R^2 - (r1^2 + 2r2^2)) - 120
# eq4 = r1 - 2r2
res = solve([eq1, eq2, eq3, eq4], (R, r1, r2, y))

   2-element Vector{NTuple{4, Sym}}:
    (-sqrt(185)/2, 5/2, 0, -sqrt(185)/2)
    (sqrt(185)/2, 5/2, 0, sqrt(185)/2)

この連立方程式は正しいにも関わらず，小円の半径 = 0 という不適切な解を与える。

nlsolve() で数値解を求めると解が求まる。

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")

   r1 - r2 - 5/2,  # eq1
   r2^2 - (r1 + r2)^2 + (R - r1 - y)^2,  # eq2
   r2^2 + y^2 - (R - r2)^2,  # eq3
   3*R^2 - 3*r1^2 - 6*r2^2 - 120,  # eq4

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (R, r1, r2, y) = u
   return [
       2*r1 - 2*r2 - 5,  # eq1
       r2^2 - (r1 + r2)^2 + (R - r1 - y)^2,  # eq2
       r2^2 + y^2 - (R - r2)^2,  # eq3
       3*R^2 - 3*r1^2 - 6*r2^2 - 120,  # eq4
   ]
end;

iniv = [big"10.0", 6, 4, -5]
res = nls(H, ini=iniv);
println(res);
   (BigFloat[10.56270582162731798080457930499016663309319075011358623115618193017898504536659, 6.40671194097832904571536744155326822289140318621162475377480880617077430209802, 3.90671194097832904571536744155326822289140318621162475377480880881194616675602, -5.388864105677014011100110408852971086316399261553975564921545695903341835297999], true)

   R = 10.562706;  r1 = 6.406712;  r2 = 3.906712;  y = -5.388864
   小円の直径 = 7.813423881956658

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1, r2, y) = res[1]
   @printf("R = %.6f;  r1 = %.6f;  r2 = %.6f;  y = %.6f\n", R, r1, r2, y)
   println("小円の直径 = $(Float64(2r2))")
   plot()
   circle(0, 0, R, :black)
   circle(0, R - r1, r1)
   circle(r2, y, r2, :blue)
   circle(-r2, y, r2, :blue)
   if more
       point(0, R - r1, " R-r1")
       point(R, 0, " R")
       point(r2, y, "(r2,y)")
       annotate!(0.8, -4, text("小円の直径 = " * @sprintf("%.5f", 2r2), :left, 9))
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その319）

算額（その319）

長野県長野市戸隠 戸隠山中院権現堂 安永5年(1776)
中村信弥「改訂増補 長野県の算額」(27)
http://www.wasan.jp/zoho/zoho.html

鈎股内に中鈎を引き，区切られた領域に中円，小円を入れる。中円の直径と弦の和が 48 寸，股と小円の直径の和が 52 寸であるとき，中鈎はいかほどか。

以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive, 中鈎::positive,
     長弦::positive, 短弦::positive, 中円径::positive, 小円径::positive,
     x::positive, y::positive;

eq1 = 中円径 + 長弦 - 48
eq2 = 股 + 小円径 - 52
eq3 = 鈎*股/2 - 中鈎*弦/2
eq4 = (長弦 + 股 + 中鈎)*中円径/4 + (鈎 + 短弦 + 中鈎)*小円径/4 - 中鈎*弦/2
eq5 = 短弦^2 + 中鈎^2 - 鈎^2
eq6 = 中鈎^2 + 長弦^2 - 股^2
eq7 = 長弦 + 中鈎 - 股 - 中円径
eq8 = 短弦 + 中鈎 - 鈎 - 小円径
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (鈎, 股, 弦, 中鈎, 長弦, 短弦, 中円径, 小円径))

   1-element Vector{NTuple{8, Sym}}:
    (30, 40, 50, 24, 32, 18, 16, 12)

中鈎は 24 寸

図を描くためには中円，小円の中心座標が必要であるが，上の連立方程式とまとめて解こうとすると solve() では無理である。そこで，上で得られた値に基づいて連立方程式を解く。

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive, 中鈎::positive,
     長弦::positive, 短弦::positive, 中円径::positive, 小円径::positive,
     x::positive,  # 中円の中心の x 座標
     y::positive;  # 小円の中心の y 座標

(鈎, 股, 弦, 中鈎, 長弦, 短弦, 中円径, 小円径) = (30, 40, 50, 24, 32, 18, 16, 12)
eq9 = distance(0, 鈎, 股, 0, x, 中円径//2) - (中円径//2)^2
eq10 = distance(0, 鈎, 股, 0, 小円径//2, y) - (小円径//2)^2
res2 = solve([eq9, eq10], (x, y))

   4-element Vector{Tuple{Sym, Sym}}:
    (16, 18)
    (16, 33)
    (128/3, 18)
    (128/3, 33)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, 弦, 中鈎, 長弦, 短弦, 中円径, 小円径) = (30, 40, 50, 24, 32, 18, 16, 12)
   (x, y) = (16, 18)
   @printf("中円の直径 = %.6f\n", 中円径)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(x, 中円径/2, 中円径/2)
   circle(小円径/2, y, 小円径/2, :blue)
   (x1, y1) = (中鈎*3/5, 中鈎*4/5)
   segment(0, 0, x1, y1)
   if more
       point(0, 0, " O", :green, :left, :bottom)
       point(股, 0, "  A", :green, :left, :bottom)
       point(0, 鈎, " B", :green, :left, :bottom)
       point(x1, y1, " C", :green, :left, :bottom)
       annotate!(15, 25, text("鈎:OB, 股:OA, 弦:AB, 中鈎:OC\n長弦:AC, 短弦:BC", 10, :left))
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;