算額（その441）

算額（その441）

東都下谷稲荷社 寛政元年己酉3月(1624)
黒川喜兵衛盛榮

藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

円弧の中に甲円 1 個，乙円 2 個，丙円 2 個があり，それぞれの円は隣同士は外接し，円弧を構成する外円には内接し，弦に接している。
甲円，丙円の直径がそれぞれ 12 寸，3 寸のとき，乙円の直径はいかほどか。

以下のように記号を定め方程式を解く。

外円の中心を原点とする。
外円の半径，中心座標を r0, (0, 0)
甲円の半径，中心座標を r1, (0, a + r1)
乙円の半径，中心座標を r2, (x2, a + r2)
丙円の半径，中心座標を r3, (x3, a + r3)
弦が y 軸と交わる座標を (0, a)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, a::positive,
     r1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive;

eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq3 = x2^2 + (a + r2)^2 - (r0 -r2)^2
eq4 = x3^2 + (a + r3)^2 - (r0 -r3)^2
eq5 = a + 2r1 - r0

res = solve([eq1, eq2, eq3, eq4, eq5], (a, r0, r2, x2, x3))

   2-element Vector{NTuple{5, Sym}}:
    ((-2*r1^(3/2)*sqrt(r3) - r1^2 + 3*r1*r3)/(r1 - r3), (-2*r1^(3/2)*sqrt(r3) + r1^2 + r1*r3)/(r1 - r3), 2*(-r1^(7/2)*sqrt(r3) + 2*r1^(5/2)*r3^(3/2) - r1^(3/2)*r3^(5/2) - r1^3*r3 + 2*r1^2*r3^2 - r1*r3^3)/(r1^3 - 3*r1^2*r3 + 3*r1*r3^2 - r3^3), -2*sqrt(2)*r1^(3/2)*sqrt(r3)/sqrt(-r1^(3/2)*sqrt(r3) + r1*r3), sqrt(-8*r1^(3/2)*sqrt(r3) + 8*r1*r3))
    ((2*r1^(3/2)*sqrt(r3) - r1^2 + 3*r1*r3)/(r1 - r3), (2*r1^(3/2)*sqrt(r3) + r1^2 + r1*r3)/(r1 - r3), 2*(r1^(7/2)*sqrt(r3) - 2*r1^(5/2)*r3^(3/2) + r1^(3/2)*r3^(5/2) - r1^3*r3 + 2*r1^2*r3^2 - r1*r3^3)/(r1^3 - 3*r1^2*r3 + 3*r1*r3^2 - r3^3), 2*sqrt(2)*r1^(3/2)*sqrt(r3)/sqrt(r1^(3/2)*sqrt(r3) + r1*r3), sqrt(8*r1^(3/2)*sqrt(r3) + 8*r1*r3))

2 番めのものが適解である。

names = ("a", "r0", "r2", "x2", "x3")
for (i, name) in enumerate(names)
   @printf("%s = %g;  ", name, res[2][i](r1 => 6, r3 => 1.5).evalf())
end

   a = 6;  r0 = 18;  r2 = 4;  x2 = 9.79796;  x3 = 14.6969;  

甲円，丙円の直径が 12 寸，3 寸のとき，乙円の直径は 8 寸である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3) =(12, 3) ./ 2
   (a, r0, r2, x2, x3) = ((2*r1^(3/2)*sqrt(r3) - r1^2 + 3*r1*r3)/(r1 - r3), (2*r1^(3/2)*sqrt(r3) + r1^2 + r1*r3)/(r1 - r3), 2*(r1^(7/2)*sqrt(r3) - 2*r1^(5/2)*r3^(3/2) + r1^(3/2)*r3^(5/2) - r1^3*r3 + 2*r1^2*r3^2 - r1*r3^3)/(r1^3 - 3*r1^2*r3 + 3*r1*r3^2 - r3^3), 2*sqrt(2)*r1^(3/2)*sqrt(r3)/sqrt(r1^(3/2)*sqrt(r3) + r1*r3), sqrt(8*r1^(3/2)*sqrt(r3) + 8*r1*r3))
   @printf("r1 = %g; r3 = %g\n", r1, r3)
   @printf("a = %g;  r0 = %g;  r2 = %g;  x2 = %g;  x3 = %g\n", a, r0, r2, x2, x3)
   @printf("乙円の直径 = %g\n", 2r2)
   plot()
   circle(0, 0, r0, :green)
   circle(0, a + r1, r1, :magenta)
   circle(x2, a + r2, r2, :blue)
   circle(-x2, a + r2, r2, :blue)
   circle(x3, a + r3, r3, :red)
   circle(-x3, a + r3, r3, :red)
   b = sqrt(r0^2 - a^2)
   segment(-b, a, b, a, :orange)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r0, " r0", :green, :left, :bottom, delta=delta/2)
       point(0, a, " a", :orange, :left, :top, delta=-delta/2)
       point(0, a + r1, "甲円:r1,(0,a+r1)", :magenta, :center, delta=-delta/2)
       point(x2, a + r2, "乙円:r2,(x2,a+r2)", :black, :center, delta=-delta/2)
       point(x3, a + r3, "丙円:r3,(x3,a+r3)", :black, :center, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       plot!(ylims=(-3, r0 + 0.5))
   else
      plot!(showaxis=false)
   end
end;

算額（その440）

算額（その440）

東都下谷摩利支天堂 天明8年戊申3月(1788)
東都浅草新堀端住 三浦伴次郎成喜

藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

求めるものが小頭か大頭かの違いだけで，算額（その185）と同じ問題である。プログラミングスタイルを最近のものに変えた。

等脚台形の隔斜（対角線）で区切られた領域に甲，乙，丙の円を入れる。甲円の直径が 100 寸，乙円の直径が 28 寸，丙円の直径が 45 寸のとき，大頭（台形の下底）はいかほどか。

以下のように記号を定め方程式を解く。
台形の右下隅 A，右上隅 b， の座標をそれぞれ (a, 0)，(b, y) と置く。
甲円，乙円，丙円の半径をそれぞれ，r1，r2，r3 とする。
甲円，乙円，丙円の中心から隔斜（対角線）までの距離が円の半径に等しいという連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, x1::positive, y1::positive, y::positive,
     r1::positive, r2::positive, r3::positive;

eq1 = distance(b, y, a, 0, x1, y1) - r3^2
eq2 = distance(-b, y, a, 0, x1, y1) - r3^2
eq3 = distance(-a, 0, b, y, x1, y1) - r3^2
eq4 = distance(-a, 0, b, y, 0, y - r2) - r2^2
eq5 = distance(-b, y, a, 0, 0, r1) - r1^2;

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)
   (a, b, x1, y1, y) = u
   return [
       -r3^2 + (x1 - (a^2*x1 - 2*a*b*x1 + a*y^2 - a*y*y1 + b^2*x1 + b*y*y1)/(a^2 - 2*a*b + b^2 + y^2))^2 + (-y*(a^2 - a*b - a*x1 + b*x1 + y*y1)/(a^2 - 2*a*b + b^2 + y^2) + y1)^2,  # eq1
       -r3^2 + (x1 - (a^2*x1 + 2*a*b*x1 + a*y^2 - a*y*y1 + b^2*x1 - b*y*y1)/(a^2 + 2*a*b + b^2 + y^2))^2 + (-y*(a^2 + a*b - a*x1 - b*x1 + y*y1)/(a^2 + 2*a*b + b^2 + y^2) + y1)^2,  # eq2
       -r3^2 + (x1 - (x1*(a^2 + 2*a*b + b^2 + y^2) - y*(a*y - a*y1 - b*y1 + x1*y))/(a^2 + 2*a*b + b^2 + y^2))^2 + (-y*(a^2 + a*b + a*x1 + b*x1 + y*y1)/(a^2 + 2*a*b + b^2 + y^2) + y1)^2,  # eq3
       -r2^2 + y^2*(-a*r2 - b*r2 + b*y)^2/(a^2 + 2*a*b + b^2 + y^2)^2 + (-r2 - y*(a^2 + a*b - r2*y + y^2)/(a^2 + 2*a*b + b^2 + y^2) + y)^2,  # eq4
       -r1^2 + y^2*(-a*r1 + a*y - b*r1)^2/(a^2 + 2*a*b + b^2 + y^2)^2 + (r1 - y*(a^2 + a*b + r1*y)/(a^2 + 2*a*b + b^2 + y^2))^2,  # eq5
   ]
end;

(r1, r2, r3) =(100, 28, 45) ./ 2
iniv = [big"151.0", 44, 39, 112, 144]
res = nls(H, ini=iniv)
names = ("a", "b", "x1", "y1", "y")
if res[2]
   for (name, x) in zip(names, res[1])
       @printf("%s = %g;  ", name, round(Float64(x), digits=6))
   end
   println()
else
   println("収束していない")
end

   a = 150;  b = 42;  x1 = 37.5;  y1 = 112.5;  y = 144;  

大頭の長さは 2a = 300（寸）である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) =(100, 28, 45) ./ 2
   (a, b, x1, y1, y) = res[1]
   @printf("r1 = %g;  r2 = %g;  r3 = %g\n", r1, r2, r3)
   @printf("a = %g; b = %g;  x1 = %g; y1 = %g; y = %g\n", a, b, x1, y1, y)
   plot([a, b, -b, -a, a], [0, y, y, 0, 0], color=:orange, lw=0.5)
   circle(0, r1, r1, :green)
   circle(x1, y1, r3, :blue)
   circle(-x1, y1, r3, :blue)
   circle(0, y - r2, r2, :red)
   segment(-b, y, a, 0, :orange)
   segment(-a, 0, b, y, :orange)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r1, " 甲円:r1,(0,r1)", :black, :left, :vcenter)
       point(0, y - r2, " 乙円:r2,(0,y-r2)", :black, :left, :vcenter)
       point(x1, y1, " 丙円:r3,(x1,y1,r3)", :black, :left, :vcenter)
       point(b, y, " B(b,y)", :black, :left, :bottom, delta=delta/2)
       point(a, 0, " A(a,0)", :black, :left, :bottom, delta=delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その439）

算額（その439）

東都本郷真光寺天満宮 天明7年丁未11月(1787)
東都湯嶋御徒町 池田重次郎教政

藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

外円の一部分である円弧内に，大中小の 5 個の円が入っている。大円と小円は弦に接している。
弦の長さが 8 寸，小円の直径が 1 寸のとき，矢（弦と円の距離，図では r0 - a）の長さはいかほどか。

外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (x1, a + r1)
中円の半径と中心座標を r2, (0, r0 - r2)
小円の半径と中心座標を r3, (r3, a + r3)
弦と y 軸の交点座標を (0, a)
として以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive,
     r1::positive, x1::positive,
     r2::positive, r3::positive,
     弦::positive, a::positive;

eq1 = x1^2 + (a + r1)^2 - (r0 - r1)^2
eq2 = (x1 - r3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = x1^2 + (r0 - r2 - a - r1)^2 - (r1 + r2)^2
eq4 = r3^2 + (r0 - r2 - a - r3)^2 - (r2 + r3)^2 
eq5 = (r0^2 - (弦/2)^2) - a^2;

res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, x1, r2, a))

   2-element Vector{NTuple{5, Sym}}:
    ((1536*r3^6 + 288*r3^4*弦^2 - 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 84*r3^2*弦^4 + 52*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(32*r3*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), r3*(128*r3^4 - 8*r3^2*弦^2 - 16*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 3*弦^4 + 弦^3*sqrt(16*r3^2 + 5*弦^2))/(2*(256*r3^4 + 32*r3^2*弦^2 + 弦^4)), r3*弦*(2*弦 + sqrt(16*r3^2 + 5*弦^2))/(16*r3^2 + 弦^2), r3*(96*r3^4 + 44*r3^2*弦^2 + 20*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 7*弦^4 + 3*弦^3*sqrt(16*r3^2 + 5*弦^2))/(8*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), (-1536*r3^6 + 1440*r3^4*弦^2 + 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 180*r3^2*弦^4 + 20*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(4608*r3^5 - 1280*r3^3*弦^2 + 32*r3*弦^4))
    ((1536*r3^6 + 288*r3^4*弦^2 + 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 84*r3^2*弦^4 - 52*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 + 弦^5*sqrt(16*r3^2 + 5*弦^2))/(32*r3*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), r3*(128*r3^4 - 8*r3^2*弦^2 + 16*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 3*弦^4 - 弦^3*sqrt(16*r3^2 + 5*弦^2))/(2*(256*r3^4 + 32*r3^2*弦^2 + 弦^4)), r3*弦*(2*弦 - sqrt(16*r3^2 + 5*弦^2))/(16*r3^2 + 弦^2), r3*(96*r3^4 + 44*r3^2*弦^2 - 20*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 7*弦^4 - 3*弦^3*sqrt(16*r3^2 + 5*弦^2))/(8*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), (-1536*r3^6 + 1440*r3^4*弦^2 - 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 180*r3^2*弦^4 - 20*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 + 弦^5*sqrt(16*r3^2 + 5*弦^2))/(4608*r3^5 - 1280*r3^3*弦^2 + 32*r3*弦^4))

最初の組が適解である。

矢の長さは 2*r3*(3*r3^2 - 4*sqrt(r3^2 + 20) - 24)/(9*r3^2 - 16)

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

   r3*(24*r3^2 - 3*弦^2 - 弦*sqrt(16*r3^2 + 5*弦^2))/(36*r3^2 - 弦^2)

弦の長さが 8 寸，小円の直径が 1 寸 のとき，矢の長さは 3 寸である。

(r3, 弦) = (1/2, 8)
r3*(24*r3^2 - 3*弦^2 - 弦*sqrt(16*r3^2 + 5*弦^2))/(36*r3^2 - 弦^2)

   3.0

   r0 = 4.16667;  r1 = 1.125;  x1 = 2;  r2 = 1.04167;  a = 1.16667
   矢 = 3

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, 弦) = (1/2, 8)
   (r0, r1, x1, r2, a) = ((1536*r3^6 + 288*r3^4*弦^2 - 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 84*r3^2*弦^4 + 52*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(32*r3*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)),
       r3*(128*r3^4 - 8*r3^2*弦^2 - 16*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 3*弦^4 + 弦^3*sqrt(16*r3^2 + 5*弦^2))/(2*(256*r3^4 + 32*r3^2*弦^2 + 弦^4)),
       r3*弦*(2*弦 + sqrt(16*r3^2 + 5*弦^2))/(16*r3^2 + 弦^2),
       r3*(96*r3^4 + 44*r3^2*弦^2 + 20*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 7*弦^4 + 3*弦^3*sqrt(16*r3^2 + 5*弦^2))/(8*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)),
       (-1536*r3^6 + 1440*r3^4*弦^2 + 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 180*r3^2*弦^4 + 20*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(4608*r3^5 - 1280*r3^3*弦^2 + 32*r3*弦^4))
   @printf("r0 = %g;  r1 = %g;  x1 = %g;  r2 = %g;  a = %g\n",
           r0, r1, x1, r2, a)
   @printf("矢 = %g\n", r0 - a)
   plot()
   circle(0, 0, r0, :gray)
   circle(x1, a + r1, r1)
   circle(-x1, a + r1, r1)
   circle(0, r0 - r2, r2, :blue)
   circle(r3, a + r3, r3, :orange)
   circle(-r3, a + r3, r3, :orange)
   b = sqrt(r0^2 - a^2)
   segment(-b, a, b, a)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, a + r1, "大円:r1,(x1,a+r1)", :red, :center, :bottom, delta=delta/2)
       point(0, r0 - r2, "中円:r2,(0,r0-r2)", :blue, :center, :bottom, delta=delta/2)
       point(r3, a + r3, " 小円:r3,(r3,a+r3)", :black, :left, :vcenter)
       point(0, a, " a", :black, :left, :top, delta=-delta/2)
       point(弦/2, a, "(弦/2,a)", :black, :right, :top, delta=-delta/2)
       point(0, r0, " r0", :black, :left, :bottom, delta=delta/2)
       point(r0, 0, "r0 ", :black, :right, :bottom, delta=delta/2)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5, ylims=(-0.3, 4.3))
   else
       plot!(showaxis=false)
   end
end;