算額（その529）

算額（その529）

和算図形問題あれこれ - 令和4年4月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html

正方形内に大円 1 個と小円 12 個が入っている，小円の直径が 1 寸のとき，大円の直径はいかほどか。

大円の半径と中心座標を r1, (0, 0)
小円の半径と中心座標を r2, (r2, 3r2)
とおき，以下の連立方程式を r2 を変数にしたまま解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive;

eq  = r2^2 + (4r2 - r2)^2 - (r1 + r2)^2
res = solve(eq, r1)[1]
res |> println

   r2*(-1 + sqrt(10))

大円の半径(r1)は，小円の半径(r2)の √10 - 1 倍である。

したがって，小円の直径が 1 寸のとき，大円の直径は √10 - 1 ≒ 2.1622776601683795 寸である。

√10 - 1

   2.1622776601683795

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   r1 = (√10 - 1)r2
   a = 4r2
   @printf("r1 = %g;  r2 = %g\n", r1, r2)
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:black, lw=0.5)
   circle(0, 0, r1, :blue)
   circle4(r2, a - r2, r2)
   circle4(a - r2, a - r2, r2)
   circle4(a - r2, r2, r2)
    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, 4r2, " 4r2", :red, :left, :bottom, delta=delta/2)
       point(0, r1, " r1", :red, :left, :top, delta=-delta/2)
       point(r2, 3r2, "(r2,3r2)", :blue, :center, :top, delta=-delta/2)
    end
end;

算額（その528）

算額（その528）

和算図形問題あれこれ - 令和4年5月の問題-No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html

外円内に弦と小円 4 個を入れる。外円の直径が 10 寸のとき，小円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
小円の半径と中心座標を r, (x, 3r - R)
とおき，以下の連立方程式を R を変数にしたまま解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r::positive, x::positive;

eq1  = x^2 + (3r - R)^2 - (R - r)^2
eq2 = x^2 + (2R - 4r)^2 - (r + r)^2
res = solve([eq1, eq2], (r, x))

   1-element Vector{Tuple{Sym, Sym}}:
    (-(-8*R^2 + 4*R^2*(-11/2 + 5*sqrt(5)/2))/(20*R), 2*R*sqrt(-11/2 + 5*sqrt(5)/2))

小円の半径は，外円の半径の (3 - sqrt(5))/2 倍である。

res[1][1] |> simplify |> println

   R*(3 - sqrt(5))/2

外円の直径が 10 寸のとき，小円の直径は 3.819660112501051 である。

10(3 - sqrt(5))/2

   3.819660112501051

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 10//2
   (r, x) = (R*(3 - sqrt(5))/2, 2*R*sqrt(-11/2 + 5*sqrt(5)/2))
   @printf("R = %g;  r = %g;  小円の直径 = %g\n", R, r, 2r)
   plot()
   circle(0, 0, R)
   circle(x, 3r - R, r, :blue)
   circle(-x, 3r - R, r, :blue)
   circle(0, R - r, r, :blue)
   circle(0, r - R, r, :magenta)
   y = 2r - R
   segment(-sqrt(R^2 - y^2), y, sqrt(R^2 - y^2), y, :green)
    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, R - r, " R-r", :blue, :left, :vcenter)
       point(x, 3r - R, " (x,3r-R)", :blue, :left, :vcenter)
       point(0, r - R, " r-R", :blue, :left, :vcenter)
       point(0, 2r - R, " 2r-R", :blue, :left, :bottom, delta=delta/2)
    end
end;

算額（その527）

算額（その527）

和算図形問題あれこれ - 令和4年8月の問題-No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html

(7) 滋賀県 吉田博氏宅　p.116
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.

キーワード：円10個，累円

外円内に甲円，乙円，丙円，丁円，戊円が入っている。外円が 40 寸のとき，丙円を最も大きくしたときの戊円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
乙円の半径と中心座標を r2, (r2, y2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (x4, y4)
戊円の半径と中心座標を r5, (x5, y5)
とおき，以下の連立方程式を解く r1 について解く。その後 r3についての式をr1で微分し導関数が0になるときの r1 を求める。そのときの r1 で表される r5 が求める値である。
しかし，連立方程式の記号解が求まらないので，数値微分を用いて同じことをする。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, y1::positive,
   r2::positive, y2,
   r3::positive, x3::positive, y3,
   r4::positive, x4::positive, y4,
   r5::positive, x5::positive, y5;
r0 = 40//2
y1 = r0 - r1
x2 = r2
eq1  = r2^2 + y2^2 - (r0 - r2)^2
eq3  = x3^2 + y3^2 - (r0 - r3)^2
eq6  = x4^2 + y4^2 - (r0 - r4)^2
eq9  = x5^2 + y5^2 - (r0 - r5)^2

eq2  = r2^2 + (y1 - y2)^2 - (r1 + r2)^2
eq4  = x3^2 + (y1 - y3)^2 - (r1 + r3)^2
eq7  = x4^2 + (y1 - y4)^2 - (r1 + r4)^2
eq10 = x5^2 + (y1 - y5)^2 - (r1 + r5)^2

eq5  = (x3 - x2)^2 + (y3 - y2)^2 - (r3 + r2)^2
eq8  = (x4 - x3)^2 + (y4 - y3)^2 - (r4 + r3)^2
eq11 = (x5 - x4)^2 + (y5 - y4)^2 - (r5 + r4)^2;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11], (r2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5))

using NLsolve

function nls(r1, func, params...; ini = [0.0])
    if typeof(ini) <: Number
        r = nlsolve((vout, vin) -> vout[1] = func(r1, vin[1], params..., [ini]), ftol=big"1e-40")
        v = r.zero[1]
    else
        r = nlsolve((vout, vin)->vout .= func(r1, vin, params...), ini, ftol=big"1e-40")
        v = r.zero
    end
    return v, r.f_converged
end;

function H(r1, u)
    (r2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5) = u
    return [
        r2^2 + y2^2 - (20 - r2)^2,  # eq1
        r2^2 - (r1 + r2)^2 + (-r1 - y2 + 20)^2,  # eq2
        x3^2 + y3^2 - (20 - r3)^2,  # eq3
        x3^2 - (r1 + r3)^2 + (-r1 - y3 + 20)^2,  # eq4
        (-r2 + x3)^2 - (r2 + r3)^2 + (-y2 + y3)^2,  # eq5
        x4^2 + y4^2 - (20 - r4)^2,  # eq6
        x4^2 - (r1 + r4)^2 + (-r1 - y4 + 20)^2,  # eq7
        -(r3 + r4)^2 + (-x3 + x4)^2 + (-y3 + y4)^2,  # eq8
        x5^2 + y5^2 - (20 - r5)^2,  # eq9
        x5^2 - (r1 + r5)^2 + (-r1 - y5 + 20)^2,  # eq10
        -(r4 + r5)^2 + (-x4 + x5)^2 + (-y4 + y5)^2,  # eq11
   ]
end;

numerical_derivative(f, x, h=1e-10) = (f(x + h) - f(x - h)) / (2 * h)
numerical_derivative2(f, x, h=1e-10) = (f(x + h) - 2 * f(x) + f(x - h)) / h^2

r1 = big"10"
while true
    iniv = BigFloat[6.1, -13.3, 4.8, 15.3, -5.7, 3.4, 16.0, 4.6, 2.3, 15.2, 9.5]
    h = 1e-10
    r3 = nls(r1, H, ini=iniv)[1][3]
    r3l = nls(r1 - h, H, ini=iniv)[1][3]
    r3r = nls(r1 + h, H, ini=iniv)[1][3]
    numerical_dervative = (r3r - r3l)/2h
    numerical_dervative2 = (nls(r1 + h, H, ini=iniv)[1][3] -2nls(r1, H, ini=iniv)[1][3] + nls(r1 - h, H, ini=iniv)[1][3])/h^2
    abs(numerical_dervative/numerical_dervative2) < big"1e-20" && break
    r1 -= numerical_dervative/numerical_dervative2
    println("r1 = $r1")
end
println("r1 = $r1")

   r1 = 12.48538011695906438403271202135121844125836007027810600838640181720678979984704
   r1 = 12.01201496172095097990312251445735848809684416949930658448600519799368331897538
   r1 = 12.00000898032969429019237242846760671180689989079357424987938861517199168254016
   r1 = 12.00000000000504037748540527008744540545918004544761142950056638766185544102556
   r1 = 11.99999999999999999999979325438981227049694304973714727017681791708893764281108
   r1 = 11.99999999999999999999979325438981227049694304973714727017681791708893764281108

甲円の半径が　12 のとき，乙円が最大で，戊円の半径を求める。

r1 = 12
iniv = BigFloat[6.1, -13.3, 4.8, 15.3, -5.7, 3.4, 16.0, 4.6, 2.3, 15.2, 9.5]
res = nls(r1, H, ini=iniv)

 (BigFloat[7.500000000000000000000000000000000000000000000000000000000000000420264272338678, -10.00000000000000000000000000000000000000000000000000000000000000168105708935471, 4.999999999999999999999999999999999999999999999995121239468409802261948427817372, 15.00000000000000000000000000000000000000000000004390884457514065793291839983194, 1.951504212636079095220628871952364045461575089304458517289222454862992166675505e-47, 2.999999999999999999999999999999999999999999998473295509659893374550847421982661, 15.00000000000000000000000000000000000000000000329535771623506103989092567879276, 8.000000000000000000000000000000000000000000006106817961360426501796610311927722, 1.874999999999999999999999999999999999999999998197629693673087468497217798802652, 13.12499999999999999999999999999999999999999999841365663452270670651293440891112, 12.5000000000000000000000000000000000000000000072094812253076501260111288044371], true)

戊円の半径は 1.875（直径は 3.75）である。

その他のパラメータは以下の通り。

r2 = 7.5;  y2 = -10;  r3 = 5;  x3 = 15;  y3 = 1.9515e-47;  r4 = 3;  x4 = 15;  y4 = 8;  r5 = 1.875;  x5 = 13.125;  y5 = 12.5
甲円の直径 = 24;  丙円の直径 = 10;  戊円の直径 = 3.75

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5) = res[1]
    r0 = 40//2
    x1 = 0
    y1 = r0 - r1
    x2 = r2
    @printf("r2 = %g;  y2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x4 = %g;  y4 = %g;  r5 = %g;  x5 = %g;  y5 = %g\n", r2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5)
    @printf("甲円の直径 = %g;  丙円の直径 = %g;  戊円の直径 = %g\n", 2r1, 2r3, 2r5)
    plot()
    circle(-r2, y2, r2, :blue)
    circle(x3, y3, r3, :magenta)
    circle(-x3, y3, r3, :magenta)
    circle(x4, y4, r4, :brown)
    circle(-x4, y4, r4, :brown)
    circle(x5, y5, r5, :orange)
    circle(-x5, y5, r5, :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, " 甲円", mark=false)
        point(x2, y2, "乙円", mark=false)
        point(x3, y3, "丙円", mark=false)
        point(x4, y4, "丁円", mark=false)
        point(x5, y5, "戊円", mark=false)
    end
end;