算額（その336）

算額（その336）

愛知県名古屋市熱田区 熱田神宮 文化3年(1806)
http://www.wasan.jp/aichi/atuta1.html

二等辺三角形の中に直交する斜線甲と斜線乙があり，面積の等しい円が 3 個内接している。斜線乙が 2 寸のとき，円の直径はいくらか。

図のように記号を定め，nlsolve() で連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, c::positive, d::positive,
     h::positive, w::positive, r::positive,
     x1::negative, y1::positive, x2::positive, y2::positive,
     X1::negative, X2::positive, Y2::positive, Y3::positive,
     dummy;

a = sqrt((x1 + w)^2 + y1^2)
b = sqrt((x2 - x1)^2 + (y2 - y1)^2)
c = sqrt((w - x2)^2 + y2^2)
d = sqrt(h^2 + w^2)
eq1 = a + 2 - 2w - 2r
eq2 = b + 2 - c - 2r
eq3 = 2a - (2w + a + 2)r
eq4 = 2b - (b + c + 2)r
eq5 = distance(0, h, w, 0, 0, Y3) - r^2
eq6 = distance(0, h, w, 0, X2, Y2) - r^2
eq7 = distance(-w, 0, x2, y2, X2, Y2) - r^2
eq8 = distance(-w, 0, x2, y2, 0, Y3) - r^2
eq9 = distance(-w, 0, x2, y2, X1, r) - r^2
eq10 = distance(x1, y1, w, 0, X1, r) - r^2
eq11 = distance(x1, y1, w, 0, X2, Y2) - r^2;

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")
println(eq6, ",  # eq6")
println(eq7, ",  # eq7")
println(eq8, ",  # eq8")
println(eq9, ",  # eq9")
println(eq10, ",  # eq10")
println(eq11, ",  # eq11")

   -2*r - 2*w + sqrt(y1^2 + (w + x1)^2) + 2,  # eq1
   -2*r - sqrt(y2^2 + (w - x2)^2) + sqrt((-x1 + x2)^2 + (-y1 + y2)^2) + 2,  # eq2
   -r*(2*w + sqrt(y1^2 + (w + x1)^2) + 2) + 2*sqrt(y1^2 + (w + x1)^2),  # eq3
   -r*(sqrt(y2^2 + (w - x2)^2) + sqrt((-x1 + x2)^2 + (-y1 + y2)^2) + 2) + 2*sqrt((-x1 + x2)^2 + (-y1 + y2)^2),  # eq4
   h^2*w^2*(-Y3 + h)^2/(h^2 + w^2)^2 - r^2 + (Y3 - h*(Y3*h + w^2)/(h^2 + w^2))^2,  # eq5
   -r^2 + (X2 - w*(X2*w - Y2*h + h^2)/(h^2 + w^2))^2 + (Y2 - h*(-X2*w + Y2*h + w^2)/(h^2 + w^2))^2,  # eq6
   -r^2 + (X2 - (X2*(w^2 + 2*w*x2 + x2^2 + y2^2) - y2*(X2*y2 - Y2*w - Y2*x2 + w*y2))/(w^2 + 2*w*x2 + x2^2 + y2^2))^2 + (Y2 - y2*(X2*w + X2*x2 + Y2*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq7
   -r^2 + y2^2*(Y3*w + Y3*x2 - w*y2)^2/(w^2 + 2*w*x2 + x2^2 + y2^2)^2 + (Y3 - y2*(Y3*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq8
   -r^2 + (X1 - (X1*(w^2 + 2*w*x2 + x2^2 + y2^2) - y2*(X1*y2 - r*w - r*x2 + w*y2))/(w^2 + 2*w*x2 + x2^2 + y2^2))^2 + (r - y2*(X1*w + X1*x2 + r*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq9
   -r^2 + (X1 - (X1*w^2 - 2*X1*w*x1 + X1*x1^2 - r*w*y1 + r*x1*y1 + w*y1^2)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2 + (r - y1*(-X1*w + X1*x1 + r*y1 + w^2 - w*x1)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2,  # eq10
   -r^2 + (X2 - (X2*(w^2 - 2*w*x1 + x1^2 + y1^2) - y1*(X2*y1 + Y2*w - Y2*x1 - w*y1))/(w^2 - 2*w*x1 + x1^2 + y1^2))^2 + (Y2 - y1*(-X2*w + X2*x1 + Y2*y1 + w^2 - w*x1)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2,  # eq11

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)
   (h, w, r, x1, y1, x2, y2, X1, X2, Y2, Y3) = u
   return [
       -2*r - 2*w + sqrt(y1^2 + (w + x1)^2) + 2,  # eq1
       -2*r - sqrt(y2^2 + (w - x2)^2) + sqrt((-x1 + x2)^2 + (-y1 + y2)^2) + 2,  # eq2
       -r*(2*w + sqrt(y1^2 + (w + x1)^2) + 2) + 2*sqrt(y1^2 + (w + x1)^2),  # eq3
       -r*(sqrt(y2^2 + (w - x2)^2) + sqrt((-x1 + x2)^2 + (-y1 + y2)^2) + 2) + 2*sqrt((-x1 + x2)^2 + (-y1 + y2)^2),  # eq4
       h^2*w^2*(-Y3 + h)^2/(h^2 + w^2)^2 - r^2 + (Y3 - h*(Y3*h + w^2)/(h^2 + w^2))^2,  # eq5
       -r^2 + (X2 - w*(X2*w - Y2*h + h^2)/(h^2 + w^2))^2 + (Y2 - h*(-X2*w + Y2*h + w^2)/(h^2 + w^2))^2,  # eq6
       -r^2 + (X2 - (X2*w^2 + 2*X2*w*x2 + X2*x2^2 + Y2*w*y2 + Y2*x2*y2 - w*y2^2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2 + (Y2 - y2*(X2*w + X2*x2 + Y2*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq7
       -r^2 + y2^2*(Y3*w + Y3*x2 - w*y2)^2/(w^2 + 2*w*x2 + x2^2 + y2^2)^2 + (Y3 - y2*(Y3*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq8
       -r^2 + (X1 - (X1*(w^2 + 2*w*x2 + x2^2 + y2^2) - y2*(X1*y2 - r*w - r*x2 + w*y2))/(w^2 + 2*w*x2 + x2^2 + y2^2))^2 + (r - y2*(X1*w + X1*x2 + r*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq9
       -r^2 + (X1 - (X1*(w^2 - 2*w*x1 + x1^2 + y1^2) - y1*(X1*y1 + r*w - r*x1 - w*y1))/(w^2 - 2*w*x1 + x1^2 + y1^2))^2 + (r - y1*(-X1*w + X1*x1 + r*y1 + w^2 - w*x1)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2,  # eq10
       -r^2 + (X2 - (X2*w^2 - 2*X2*w*x1 + X2*x1^2 - Y2*w*y1 + Y2*x1*y1 + w*y1^2)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2 + (Y2 - y1*(-X2*w + X2*x1 + Y2*y1 + w^2 - w*x1)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2,  # eq11
   ]
end;

iniv = [big"4.3", 1.3, 0.4, -0.3, 0.9, 0.5, 2.2, -0.2, 0.4, 1.2, 2.4]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[4.285714285714285714285714285714285714285714285714285714285574246330323963257686, 1.249999999999999999999999999999999999999999999999999999999958894684739139553026, 0.4999999999999999999999999999999999999999999999999999999999816989607612493273477, -0.3499999273461119877023840155306407964015055211277256589619837238743919315567991, 1.199999945509583990776788011647980597301129140845794244221458224297416065202069, 0.5499999999999999999999999999999999999999999999999999999999544520988613601022169, 2.399999999999999999999999999999999999999999999999999999999895684919240762427152, -0.2499999999999999999999999999999999999999999999999999999999983023984933316200398, 0.3499999999999999999999999999999999999999999999999999999999877259401613410750471, 1.299999999999999999999999999999999999999999999999999999999966810926003067216378, 2.499999999999999999999999999999999999999999999999999999999924978510268800929172], true)

names = ["h", "w", "r", "x1", "y1", "x2", "y2", "X1", "X2", "Y2", "Y3"]
for (i,  name) in enumerate(names)
   @printf("%2s = %.6f\n", name, res[1][i])
end

    h = 4.285714
    w = 1.250000
    r = 0.500000
   x1 = -0.350000
   y1 = 1.200000
   x2 = 0.550000
   y2 = 2.400000
   X1 = -0.250000
   X2 = 0.350000
   Y2 = 1.300000
   Y3 = 2.500000

円の直径は 2r = 1 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (h, w, r, x1, y1, x2, y2, X1, X2, Y2, Y3) = res[1]
   plot([w, 0, -w, w], [0, h, 0, 0], color=:black, lw=0.5)
   circle(X1, r, r)
   circle(X2, Y2, r)
   circle(0, Y3, r)
   segment(-w, 0, x2, y2)
   segment(x1, y1, w, 0)
   if more
       point(0, Y3, " Y3")
       point(X2, Y2, "(X2,Y2)\n", :green, :center, :bottom)
       point(X1, r, "(X1,r)\n", :green, :center, :bottom)
       point(-w, 0, "-w")
       point(w, 0, "w")
       point(0, h, " h")
       point(x1, y1, "(x1,y1)", :green, :right, :bottom)
       point(x2, y2, " (x2,y2)")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その335）

算額（その335）

石川県能登町 白山神社
【算額に挑戦】石川県の算額 －その１－
http://blog.livedoor.jp/enjoy_math/archives/51203926.html

等脚台形に甲円，乙円が入っている。それぞれは台形の辺に3点で内接し，また界斜にも内接している。
上頭は 12寸，甲円，乙円の直径はそれぞれ 24寸，15寸である。界斜の長さはいかほどか。

元の図を反時計回りに 90° 回転したほうが若干わかりやすいか。
甲円，乙円の半径を r2, r1 とする。等脚台形の脚を延長し，交点を原点 O とする。∠AOC を θとする。AC = 上頭/2 = 6 である。
AO = a, OC = b として，a, b を求める。

include("julia-source.txt");

using SymPy

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

(r1, r2, c) = (15//2, 24//2, 6)
sinθ = c/b
eq1 = (a + r1)sinθ - r1
eq2 = a^2 + c^2 - b^2

solve([eq1, eq2], (a, b))[1]

   (80/3, 82/3)

(a, b) = (80//3, 82//3)
(sinθ, tanθ, cosθ) = (c/b, c/a, a/b)

   (9//41, 9//40, 40//41)

界斜と脚の交点 (x1, y1), (x2, y2) を求める。

@syms x1::positive, y1::positive, x2::positive, y2::negative
y1 = x1*tanθ
y2 = -x2*tanθ
eq3 = distance(x1, y1, x2, y2, a + r1, 0) - r1^2
eq4 = distance(x1, y1, x2, y2, 2b, 0) - r2^2
solve([eq3, eq4])

   2-element Vector{Dict{Any, Any}}:
    Dict(x1 => 5200/123 - 40*sqrt(10)/41, x2 => 40*sqrt(10)/41 + 5200/123)
    Dict(x1 => 40*sqrt(10)/41 + 5200/123, x2 => 5200/123 - 40*sqrt(10)/41)

2 番目のものが適解である（1 番目のものと対称）。

(x1, x2) = (5200//123 + 40*sqrt(Sym(10))/41, 5200//123 - 40*sqrt(Sym(10))/41)
(y1, y2) = (x1, -x2) .* tanθ
(x1, y1, x2, y2) |> println

   (40*sqrt(10)/41 + 5200/123, 9*sqrt(10)/41 + 390/41, 5200/123 - 40*sqrt(10)/41, -390/41 + 9*sqrt(10)/41)

界斜の長さは，2 点間の距離を求めればよい。

界斜の長さ = sqrt((x1 - x2)^2 + (y1 - y2)^2)
@printf("界斜の長さ = %.6g\n", 界斜の長さ)

   界斜の長さ = 20

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   segment(0, 0, (2b + r2), (2b + r2)*tanθ, :green)
   segment(0, 0, (2b + r2), -(2b + r2)*tanθ, :green)
   segment((2b + r2), (2b + r2)*tanθ, (2b + r2), -(2b + r2)*tanθ)
   segment(a, c, a, -c)
   circle(a + r1, 0, r1)
   segment(a + r1, 0, (a + r1)*cosθ^2, (a + r1)*cosθ*sinθ)
   circle(2b, 0, r2)
   segment(2b, 0, 2b*cosθ^2, 2b*cosθ*sinθ)
   segment(x1, y1, x2, y2, :blue)
   if more
       point(a, 0, "a ", :green, :right)
       point(a, c, "(a,c) ", :green, :right, :bottom)
       point(a + r1, 0, " a+r1")
       point(2b, 0, "2b ", :green, :right)
       point(2b + r2, 0, "2b+r2 ", :green, :right)
       point(2b + r2, (2b + r2)*tanθ, "(2b+r2,(2b+r2)*tanθ", :green, :right, :bottom)
       point(x1, y1, "(x1,y1)", :blue, :right, :bottom)
       point(x2, y2, "(x2,y2)", :blue, :right, :top)
       point(0, 0, " O", :red)
       point(a, 0, " A", :red, :left)
       point(a, c, " C", :red, :left, :bottom)
       annotate!(0.5, 15, text("a = 80/3, b = (0,0)-(a,c) = 82/3, c = 6\n" *
               "∠AOC をθとおく\n2b+r2 = 200/3, (2b + r2)*tanθ = 15\n",
               :black, :left, :top, 9))
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;