算額（その564）

算額（その564）

一 群馬県高崎市石原町 清水寺 享保20年(1735)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

群馬の算額 1 清水寺
http://takasakiwasan.web.fc2.com/yattemimasennka1.html

キーワード：円6個，外円

外円内に大円 2 個，中円 1 個，小円 2 個が入っている。大円，中円，小円の直径が与えられたとき，外円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (r1, y1)
中円の半径と中心座標を r2, (0, y2)
小円の半径と中心座標を r3, (r3, y3)
とおき，以下の連立方程式を解く。

流石に，未知数が 4 個とはいえ，他の 3 変数を変数のまま連立方程式を解くのは無理なようで，r1, r2, r3 を数値で与えれば解ける。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, y1::negative, r2::positive, y2::positive, r3::positive, y3::positive
(r1, r2, r3) = (3, 2, 1)
eq1 = r1^2 + y1^2 - (R - r1)^2
eq2 = r3^2 + y3^2 - (R - r3)^2
eq3 = r1^2 + (y2 - y1)^2 - (r1 + r2)^2
eq4 = r3^2 + (y3 - y2)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3, eq4], (R, y1, y2, y3))

   1-element Vector{NTuple{4, Sym}}:
    (22/7 + 16*sqrt(2)/7, -8/7 - 2*sqrt(2)/7, 20/7 - 2*sqrt(2)/7, 12*sqrt(2)/7 + 20/7)

大円，中円，小円の直径がそれぞれ 6, 4, 2 のとき，外円の直径は (22 + 16√2)/7 である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, y1, y2, y3) = (22/7 + 16*sqrt(2)/7, -8/7 - 2*sqrt(2)/7, 20/7 - 2*sqrt(2)/7, 12*sqrt(2)/7 + 20/7)
   plot()
   circle(0, 0, R)
   circle(r1, y1, r1, :blue)
   circle(-r1, y1, r1, :blue)
   circle(0, y2, r2, :green)
   circle(r3, y3, r3, :magenta)
   circle(-r3, y3, r3, :magenta)
   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(r1, y1, "大円:r1,(r1,y1)", :blue, :center, delta=-delta/2)
       point(0, y2, " 中円:r2\n (0,y2)", :green, :left, :vcenter)
       point(r3, y3, " 小円:r3,(r3,y3)", :magenta, :left, :vcenter)
       point(R, 0, "R ", :red, :right, :bottom, delta=delta/2)
   end
end;

 

算額（その563）

算額（その563）

群馬の算額 20 旧福寿庵
http://takasakiwasan.web.fc2.com/gunnsann/g020.html

長方形内に 2 本の斜線，上円 1 個，下円　2 個が入っている。長方形の縦×横が 336 平方寸，上円と下円の直径の積が 63 平方寸のとき，上円の直径はいかほどか。

長方形の短辺，長辺の長さを a, b とおく。
上円の半径と中心座標を r1, (a/2, b - r1)
下円の半径と中心座標を r2, (r2, r2)
とおき，以下の連立方程式を解く。
簡単な連立方程式に見えるが SymPy では解けないので数値解を求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, a::positive, b::positive

弦 = sqrt(b^2 + a^2/4)
eq1 = r1*弦 - a*(b - r1)/2
eq2 = b + a/2 - 弦 - 2r2
eq3 = a*b - 336
eq4 = 2r1*2r2 - 63;
# res = solve([eq1, eq2, eq3, eq4], (r1, r2, a, b))

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, r2, a, b) = u
   return [
       -a*(b - r1)/2 + r1*sqrt(a^2/4 + b^2),  # eq1
       a/2 + b - 2*r2 - sqrt(a^2/4 + b^2),  # eq2
       a*b - 336,  # eq3
       4*r1*r2 - 63,  # eq4
   ]
end;

iniv = BigFloat[5, 3.1, 15, 25]
res = nls(H, ini=iniv)

   (BigFloat[5.250000000000000000000000000000000000000000009939436168066665434919743445334633, 2.999999999999999999999999999999999999999999989678292454450465893181229233646876, 13.99999999999999999999999999999999999999999983811297245359217105296244127297891, 24.00000000000000000000000000000000000000000008495905738224684384339759581898053], true)

上円の半径は 5.25 寸，直径は 10.5 寸である。

その他のパラメータは以下の通り。
上円の直径 = 10.5;  r1 = 5.25;  r2 = 3;  a = 14;  b = 24

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, a, b) = res[1]
   @printf("上円の直径 = %g;  r1 = %g;  r2 = %g;  a = %g;  b = %g\n", 2r1, r1, r2, a, b)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:blue, lw=0.5)
   circle(a/2, b - r1, r1)
   circle(r2, r2, r2, :green)
   circle(a - r2, r2, r2, :green)
   plot!([0, a/2, a], [b, 0, b], color=:magenta, lw=0.5)
   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(a/2, b - r1, "上円:r1,(a/2,b-r1)", :red, :center, delta=-delta)
       point(r2, r2, "下円:r2,(r2,r2)", :green, :center, delta=-delta)
       point(a, 0, "a ", :black, :right, :bottom, delta=delta/2)
       point(a/2, 0, "a/2 ", :black, :right, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :bottom, delta=delta/2)
   end
end;