算額（その803）

算額（その803）

藤田貞資：精要算法（下巻），天明元年(1781) http://www.wasan.jp/seiyou/seiyou.html

三辺（大斜，中斜，小斜）が 28 寸，25 寸，17 寸の不等辺三角形の内部の 1 点から 3 つの頂点を結ぶ線分を引き，できる三角形の面積が等しくなるようにするとき，大斜と小斜の交点から件の点までの距離（甲斜）を求めよ。

大斜と小斜の交点を原点とし，大斜，中斜，小斜の長さを L, M, S とおき，三角形内部の 1 点の座標を (x1, y1)，中斜と小斜の交点を (x2, y2) として，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms x1, y1, x2, y2, L, M, S, d
eq1 = x2^2 + y2^2 - S^2
eq2 = (L - x2)^2 + y2^2 - M^2
eq3 = S*sqrt(dist(0, 0, x2, y2, x1, y1)) - L*y1
eq3 = dist(0, 0, x2, y2, x1, y1) - (L*y1/S)^2
eq3 = numerator(apart(eq3, d))
eq4 = M*sqrt(dist(L, 0, x2, y2, x1, y1)) - L*y1
eq4 = dist(L, 0, x2, y2, x1, y1) - (L*y1/M)^2
eq4 = numerator(apart(eq4, d))
res = solve([eq1, eq2, eq3, eq4], (x1, y1, x2, y2));
res[2]

   ((3*L^2 - M^2 + S^2)/(6*L), sqrt(-L^4 + 2*L^2*M^2 + 2*L^2*S^2 - M^4 + 2*M^2*S^2 - S^4)/(6*L), (L^2 - M^2 + S^2)/(2*L), sqrt(-(L - M - S)*(L - M + S)*(L + M - S)*(L + M + S))/(2*L))

2 組の解が得られるが， 2 番目のものが適解である。

甲斜はピタゴラスの定理で計算する。
甲斜は，大斜の 2 乗の 2 倍と，小斜の 2 乗の 2 倍の和から中斜の 2 乗を引き，平方根を取ったものを 3 で割ることで得られる。

甲斜 = sqrt(res[2][1]^2 + res[2][2]^2) |> simplify
甲斜 |> println

   sqrt(2*L^2 - M^2 + 2*S^2)/3

大斜，中斜，小斜が 28 寸，25 寸，17 寸のとき，甲斜は 13 寸である。

甲斜(L => 28, M => 25, S => 17).evalf() |> println

   13.0000000000000

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (L, M, S) = (28, 25, 17)
   (x1, y1, x2, y2) = ((3*L^2 - M^2 + S^2)/(6*L), sqrt(-L^4 + 2*L^2*M^2 + 2*L^2*S^2 - M^4 + 2*M^2*S^2 - S^4)/(6*L), (L^2 - M^2 + S^2)/(2*L), sqrt(-(L - M - S)*(L - M + S)*(L + M - S)*(L + M + S))/(2*L))
   plot([0, 28, x2, 0], [0, 0, y2, 0], color=:magenta, lw=0.5)
   segment(0, 0, x1, y1, :red)
   segment(28, 0, x1, y1, :blue)
   segment(x2, y2, x1, y1, :blue)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(x1, y1, "(x1,y1)", :blue, :left, :bottom, delta=delta)
       point(x2, y2, "(x2,y2)", :magenta, :left, :bottom, delta=delta)
       point(L/2, 0, "大斜", :black, :center, :bottom, delta=delta, mark=false)
       point((L + x2)/2, y2/2, "中斜", :black, :left, :bottom, delta=delta, mark=false)
       point(x2/2, y2/2, "小斜", :black, :right, :bottom, delta=delta, mark=false)
       point(x1/2, y1/2, "甲斜", :red, :right, :bottom, delta=delta, mark=false)
   end
end;

算額（その802）

算額（その802）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

菱形の中に大円，中円，小円が入っている。小円と中円の直径の積が 1 寸（平方寸），双股（小円と中円が菱形と接する 2 点間の距離）が 5 寸のとき，大円の直径はいかほどか。

菱形の対角線を 2a, 2b
大円の半径と中心座標を r1, (0, 0)
中円の半径と中心座標を r2, (0, r1 + r2)
小円の半径と中心座標を r3, (r1 + r3, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive, r3::positive
c = sqrt(a^2 + b^2)
sinθ = b/c
cosθ = a/c
eq1 = 2r2*2r3 - 1
eq2 = c - (c - r1 - r2)*sinθ - (c -r1 - r3)*cosθ - 5//2
eq3 = r1 - a*sinθ
eq4 = r2 - (b - r1 - r2)*cosθ
eq5 = r3 - (a - r1 - r3)*sinθ;

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 Float64.(v), r.f_converged
end;

function H(u)
   (a, b, r1, r2, r3) = u
   return [
       4*r2*r3 - 1,  # eq1
       -a*(-r1 - r3 + sqrt(a^2 + b^2))/sqrt(a^2 + b^2) - b*(-r1 - r2 + sqrt(a^2 + b^2))/sqrt(a^2 + b^2) + sqrt(a^2 + b^2) - 2.5,  # eq2
       -a*b/sqrt(a^2 + b^2) + r1,  # eq3
       -a*(b - r1 - r2)/sqrt(a^2 + b^2) + r2,  # eq4
       -b*(a - r1 - r3)/sqrt(a^2 + b^2) + r3,  # eq5
   ]
end;

iniv = BigFloat[4, 5, 0.3, 3, 0.74]
res = nls(H, ini=iniv)

   ([5.0, 3.75, 3.0, 0.3333333333333333, 0.75], true)

大円の直径は 6 寸である。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, r1, r2, r3) = res[1]
   @printf("大円の直径 = %g\n", 2r1)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:magenta, lw=0.5)
   circle(0, 0, r1)
   circle(0, r1 + r2, r2, :blue)
   circle(r1 + r3, 0, r3, :green)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
   end
end;

算額（その801）

算額（その801）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

大円の周りに中円と小円が，大円と同時に隣の円とも外接するように配置されている。
大円と小円の直径がそれぞれ 21707 寸，5759 寸のとき，中円の直径はいかほどか。
なお，大円と小円の直径が上述の数値の場合には精要算法に示された図とはイメージがかなり異なる。小円の直径が 6800 寸程度の場合によく似た図（下図）になる。

大円の半径と中心座標を r1, (0, 0)
中円の半径と中心座標を r2, (0, r1 + r2)
小円の半径と中心座標を r3, (x3, y3), (x4, y4), (x5, y5)
中円と小円の中心から大円の中心へ引いた2本の直線のなす角を θ
とおき，以下の連立方程式を解く。
eq3，eq4, eq5 は第二余弦定理を援用した。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1, r2, r3, x3, y3, x4, y4, x5, y5, θ
eq1 = x3^2 + y3^2 - (r1 + r3)^2
eq2 = x3^2 + (r1 + r2 - y3)^2 - (r2 + r3)^2
eq3 = (x4 - x3)^2 + (y3 - y4)^2 - 4r3^2
eq4 = (r1 + r2)^2 + (r1 + r3)^2 - 2(r1 + r2)*(r1 + r3)cos(θ) - (x3^2 + (r1 + r2 - y3)^2)
eq5 = 2(r1 + r3)^2 - 2(r1 + r3)^2*cos(PI/3 - θ) - ((x3 - x4)^2 + (y3 - y4)^2)
eq6 = (r1 + r2)^2 + (r1 + r3)^2 - 2(r1 + r2)*(r1 + r3)cos(PI/3) - (x4^2 + (r1 + r2 - y4)^2);

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (r2, x3, y3, x4, y4, θ))

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 Float64.(v), r.f_converged
end;

function H(u)
   (r2, x3, y3, x4, y4, θ) = u
   return [
       x3^2 + y3^2 - (r1 + r3)^2,  # eq1
       x3^2 - (r2 + r3)^2 + (r1 + r2 - y3)^2,  # eq2
       -4*r3^2 + (-x3 + x4)^2 + (y3 - y4)^2,  # eq3
       -x3^2 + (r1 + r2)^2 + (r1 + r3)^2 - (r1 + r3)*(2*r1 + 2*r2)*cos(θ) - (r1 + r2 - y3)^2,  # eq4
       -2*(r1 + r3)^2*sin(θ + pi/6) + 2*(r1 + r3)^2 - (x3 - x4)^2 - (y3 - y4)^2,  # eq5
       -x4^2 - (r1/2 + r3/2)*(2*r1 + 2*r2) + (r1 + r2)^2 + (r1 + r3)^2 - (r1 + r2 - y4)^2,  # eq6
   ]
end;

(r1, r3) = (21707, 5759) .// 2
# (r1, r3) = (21707, 6800) .// 2
r1r3 = r1 + r3
iniv = BigFloat[r1r3, r1r3*cos(pi/3), r1r3*sin(pi/3), r1r3*cos(pi/6), r1r3*sin(pi/6), pi/6]
res = nls(H, ini=iniv)

   ([8893.500000000004, 8031.882766640407, 11139.306451612903, 11893.126870171644, 6866.499999999951, 0.624707483913989], true)

中円の直径は 17787 寸である。

その他のパラメータは以下のとおりである。

r2 = 8893.5;  x3 = 8031.88;  y3 = 11139.3;  x4 = 11893.1;  y4 = 6866.5;  θ = 0.624707

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, x3, y3, x4, y4, θ) = res[1]
   @printf("中円の直径 = %g\n", 2r2)
   @printf("r2 = %g;  x3 = %g;  y3 = %g;  x4 = %g;  y4 = %g;  θ = %g\n", r2, x3, y3, x4, y4, θ)
   plot()
   circle(0, 0, r1)
   rotate(0, r1 + r2, r2, :blue)
   #circle2(x3, y3, r3, :green)
   #circle2(x4, y4, r3, :green)
   rotate(x3, y3, r3, :green)
   rotate(x4, y4, r3, :green)
   x5 = (r1 + r3)*cos(θ - pi/6)
   y5 = (r1 + r3)*sin(θ - pi/6)
   rotate(x5, y5, r3, :green)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
   end
end;