算額（その261）

算額（その261）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（71）
長野県長野市 久保寺観音堂 享和3年(1803)

十七 群馬県高崎市八幡町 文化7年(1810)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

キーワード：円2個，直角三角形，中鈎
#Julia, #SymPy, #算額, #和算

鈎股弦（直角三角形）の中鈎を隔てて，大円，小円が入っている。大円，小円の直径が 68.8 寸，51.6 寸であるとき，鈎，股，弦を求めよ。

鈎，股，弦，中鈎，短弦，長弦 をそのまま「鈎」，「股」，「弦」，「中鈎」，「短弦」，「長弦」
大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (股 - r2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive,
     中鈎::positive, 長弦::positive, 短弦::positive,
     x1::positive, y2::positive,
     r1::positive, r2::positive;

eq1 = 鈎^2 + 股^2 - 弦^2

# eq2 = 短弦^2 + 中鈎^2 - 鈎^2
# eq3 = 長弦^2 + 中鈎^2 - 股^2
eq2 = 短弦/鈎 - 鈎/弦
eq3 = 中鈎/股 - 鈎/弦

eq4 = 弦 - 長弦 - 短弦

# eq5 = dist2(長弦*股/弦, 長弦*鈎/弦, 股, 0, x1, r1, r1)
# eq6 = dist2(長弦*股/弦, 長弦*鈎/弦, 股, 0, 股 - r2, y2, r2)
eq5 = 長弦 + 中鈎 - 股 - 2r1
eq6 = 短弦 + 中鈎 - 鈎 - 2r2

eq7 = dist2(0, 0, 股, 鈎, x1, r1, r1)
eq8 = dist2(0, 0, 股, 鈎, 股 - r2, y2, r2);

SymPy の性能上，一度に解けないので，まず eq1 〜 eq6 を解いて 鈎, 股, 弦, 中鈎, 長弦, 短弦 を求める。

res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (鈎, 股, 弦, 中鈎, 長弦, 短弦))[2]  # 2 of 2

   ((r1^4 + r1^3*r2 + r1^3*sqrt(r1^2 + r2^2) + 3*r1^2*r2^2 + 2*r1^2*r2*sqrt(r1^2 + r2^2) + r1*r2^3 + r1*r2^2*sqrt(r1^2 + r2^2) + 2*r2^4 + 2*r2^3*sqrt(r1^2 + r2^2))/(r1*(r1^2 + r2^2 + r2*sqrt(r1^2 + r2^2))), (r1^4 + r1^3*r2 + r1^3*sqrt(r1^2 + r2^2) + 3*r1^2*r2^2 + 2*r1^2*r2*sqrt(r1^2 + r2^2) + r1*r2^3 + r1*r2^2*sqrt(r1^2 + r2^2) + 2*r2^4 + 2*r2^3*sqrt(r1^2 + r2^2))/(r2*(r1^2 + r2^2 + r2*sqrt(r1^2 + r2^2))), (r1^2*(r1 + sqrt(r1^2 + r2^2)) + r1*r2*(r1 + r2) + r2^2*(r2 + sqrt(r1^2 + r2^2)))/(r1*r2), (r1^3 + 2*r1^2*r2 + r1^2*sqrt(r1^2 + r2^2) + r1*r2^2 + r1*r2*sqrt(r1^2 + r2^2) + 2*r2^3 + 2*r2^2*sqrt(r1^2 + r2^2))/(r1^2 + r2^2 + r2*sqrt(r1^2 + r2^2)), r1*(r1 + r2 + sqrt(r1^2 + r2^2))/r2, r2*(r1 + r2 + sqrt(r1^2 + r2^2))/r1)

2 組の解が得られるが，2 番目のものが適解である。弦を表す式は以下のように簡約化できる。

ans_鈎 = res[1] |> factor
ans_鈎 |> println

   (r1^2 + r2^2)*(r1^2 + r1*r2 + r1*sqrt(r1^2 + r2^2) + 2*r2^2 + 2*r2*sqrt(r1^2 + r2^2))/(r1*(r1^2 + r2^2 + r2*sqrt(r1^2 + r2^2)))

ans_股 = res[2] |> factor
ans_股 |> println

   (r1^2 + r2^2)*(r1^2 + r1*r2 + r1*sqrt(r1^2 + r2^2) + 2*r2^2 + 2*r2*sqrt(r1^2 + r2^2))/(r2*(r1^2 + r2^2 + r2*sqrt(r1^2 + r2^2)))

ans_弦 = res[3] |> factor
ans_弦 |> println

   (r1^2 + r2^2)*(r1 + r2 + sqrt(r1^2 + r2^2))/(r1*r2)

大円，小円の直径が 68.8 寸，51.6 寸であるとき，鈎 = 129 寸，股 = 172 寸，弦 = 215 寸である。

ans_鈎(r1 => 68.8/2, r2 => 51.6/2).evalf() |> println
ans_股(r1 => 68.8/2, r2 => 51.6/2).evalf() |> println
ans_弦(r1 => 68.8/2, r2 => 51.6/2).evalf() |> println

   129.000000000000
   172.000000000000
   215.000000000000

鈎，股の解が得られたので，eq7, eq8 の鈎，股に代入し，連立方程式を解いて x1, y2 を求める。

eq17 = eq7(鈎 => res[1], 股 => res[2]) |> simplify |> numerator;
eq18 = eq8(鈎 => res[1], 股 => res[2]) |> simplify |> numerator;
res2 = solve([eq17, eq18], (x1, y2))[3]  # 3 of 4

   (r1^2/r2 + r1*sqrt(r1^2 + r2^2)/r2, (r1^4 + 2*r1^3*r2 + r1^3*sqrt(r1^2 + r2^2) + 4*r1^2*r2^2 + 3*r1^2*r2*sqrt(r1^2 + r2^2) + 2*r1*r2^3 + 2*r1*r2^2*sqrt(r1^2 + r2^2) + 2*r2^4 + 2*r2^3*sqrt(r1^2 + r2^2) - r2*sqrt(r1^6 + 9*r1^4*r2^2 + 4*r1^4*r2*sqrt(r1^2 + r2^2) + 16*r1^2*r2^4 + 12*r1^2*r2^3*sqrt(r1^2 + r2^2) + 8*r2^6 + 8*r2^5*sqrt(r1^2 + r2^2)))/(r1^3 + 2*r1*r2^2 + 2*r1*r2*sqrt(r1^2 + r2^2)))

4 組の解が得られるが 3 番目のものが適解である。

x
res2[1](r1 => 68.8/2, r2 => 51.6/2) |>  println

   103.200000000000

y
res2[2](r1 => 68.8/2, r2 => 51.6/2) |>  println

   77.4000000000000

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = sqrt(r1^2 + r2^2)
   (鈎, 股, 弦, 中鈎, 長弦, 短弦) = ((r1^4 + r1^3*r2 + r1^3*t + 3*r1^2*r2^2 + 2*r1^2*r2*t + r1*r2^3 + r1*r2^2*t + 2*r2^4 + 2*r2^3*t)/(r1*(r1^2 + r2^2 + r2*t)), (r1^4 + r1^3*r2 + r1^3*t + 3*r1^2*r2^2 + 2*r1^2*r2*t + r1*r2^3 + r1*r2^2*t + 2*r2^4 + 2*r2^3*t)/(r2*(r1^2 + r2^2 + r2*t)), (r1^2*(r1 + t) + r1*r2*(r1 + r2) + r2^2*(r2 + t))/(r1*r2), (r1^3 + 2*r1^2*r2 + r1^2*t + r1*r2^2 + r1*r2*t + 2*r2^3 + 2*r2^2*t)/(r1^2 + r2^2 + r2*t), r1*(r1 + r2 + t)/r2, r2*(r1 + r2 + t)/r1)
   (x1, y2) = (r1^2/r2 + r1*t/r2, (r1^4 + 2*r1^3*r2 + r1^3*t + 4*r1^2*r2^2 + 3*r1^2*r2*t + 2*r1*r2^3 + 2*r1*r2^2*t + 2*r2^4 + 2*r2^3*t - r2*sqrt(r1^6 + 9*r1^4*r2^2 + 4*r1^4*r2*t + 16*r1^2*r2^4 + 12*r1^2*r2^3*t + 8*r2^6 + 8*r2^5*t))/(r1^3 + 2*r1*r2^2 + 2*r1*r2*t))
   @printf("鈎 = %g;  股 = %g;  弦 = %g;  中鈎 = %g;  長弦 = %g;  短弦 = %g; x1 = %g;  y2 = %g\n", 鈎, 股, 弦, 中鈎, 長弦, 短弦, x1, y2)
   # println("2r1 + 4r2 = $(2r1 + 4r2)")
   plot([股, 股, 0, 股, 長弦*股/弦], [ 0, 鈎, 0, 0, 長弦*鈎/弦], linecolor=:black, linewidth=0.5)
   circle(x1, r1, r1)
   circle(股 - r2, y2, r2, :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, r1, "大円:r1\n(x1,r1)", :red, :center, delta=-delta)
       point(股 - r2, y2, "小円:r2\n(股-r2,y2)", :blue, :center, delta=-delta)
       point(股, 鈎/2, " 鈎", mark=false)
       point(股/2, 0, "股  ", :green, :right, :bottom, mark=false)
       point((股 + 長弦*股/弦)/2, 長弦*鈎/弦/2, "  中鈎", mark=false)
       point(長弦*股/弦/2, 長弦*鈎/弦/2, "長弦   \n", :green, :right, mark=false)
       point((股 + 長弦*股/弦)/2, (鈎 + 長弦*鈎/弦)/2, "短弦\n", :green, :right, :bottom, mark=false)
       point(股/2, 鈎/2, "弦 = 長弦 + 短弦   ", :green, :right, mark=false)
       xlims!(-3delta, 股 + 5delta)
   end
end;

draw(68.8/2, 51.6/2, true)

 

算額（その260）

算額（その260）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（70）
長野県長野市 久保寺観音堂 享和3年(1803)

圭（二等辺三角形）内に，大円 1 個，小円 2 個が入っている。
中鉤（圭の高さ）が 36 寸，大円と小円 2 個の径の和が 44寸のとき，大円の径を求めよ。

大円と小円の半径をそれぞれ r1, r2，圭の底辺の長さを 2x とする。
以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

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

l = sqrt(36^2 + x^2)
eq1 = 2r1 + 4r2 - 44  # 大円の径と小円の径の 2 倍（2 個ぶん）の和が 44 寸
eq2 = (l + x)r1 - 36x  # 面積との関係
eq3 = (36 + l + x)r2 - 36x  # 面積との関係

res = solve([eq1, eq2, eq3], (r1, r2, x))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (10, 6, 15)

大円の径 = 20 寸，小円の径 = 12 寸，底辺の長さ = 30 寸である。

using Printf
using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x) = res[1]
   @printf("r1 = %.5f;  r2 = %.5f;  x = %.5f\n", r1, r2, x)
   println("2r1 + 4r2 = $(2r1 + 4r2)")
   plot([x, 0, -x, x], [0, 36, 0, 0], linecolor=:black, linewidth=0.5)
   circle(0, r1, r1)
   circle(r2, r2, r2, :blue)
   circle(-r2, r2, r2, :blue)
   if more == true
       point(0, 36, " 36")
       point(x, 0, " x", :green, :left, :bottom)
       point(0, r1, " 大円:r1")
       point(r2, r2, " 小円:r1\n (r2,r2)")
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

A = 中鉤, B = 径の和 として A, B を変数のまま eq1, eq2, eq3 を解く。4 組の解が得られるが，4 番目のものが適切である。

using SymPy

@syms l::positive, x::positive, r1::positive, r2::positive;
@syms A, B

l = sqrt(A^2 + x^2)
eq1 = 2r1 + 4r2 - B  # 大円の径と小円の径の 2 倍（2 個ぶん）の和が 44 寸
eq2 = (l + x)r1 - A*x
eq3 = r1*(x - r2) - r2*x

res = solve([eq1, eq2, eq3], (r1, r2, x))[4]

   ((-3*A*(4*A - B) + 3*sqrt(A*(4*A - B))*(2*A - B) + (3*A - B)*(B + sqrt((B^2*(3*A - B)^2 - 2*B*(3*A - B)*(A*(4*A - B) - sqrt(A*(4*A - B))*(2*A - B)) + 9*(A*(4*A - B) - sqrt(A*(4*A - B))*(2*A - B))^2)/(3*A - B)^2)))/(4*(3*A - B)), (3*A*(4*A - B) + 3*sqrt(A*(4*A - B))*(-2*A + B) + (3*A - B)*(B - sqrt((B^2*(3*A - B)^2 - 2*B*(3*A - B)*(A*(4*A - B) - sqrt(A*(4*A - B))*(2*A - B)) + 9*(A*(4*A - B) - sqrt(A*(4*A - B))*(2*A - B))^2)/(3*A - B)^2)))/(8*(3*A - B)), (A*(4*A - B) + sqrt(A*(4*A - B))*(-2*A + B))/(2*(3*A - B)))

同じ項が何箇所にも現れるのでそれをまとめると文字数としては少なくなるが，複雑な式であることに変わりはない。

A = 36; B = 44
t1 = A*(4*A - B)
t2 = (2*A - B)
t3 = (3*A - B)
t4 = sqrt(t1)*t2
t5 = (t1 - t4);

((-3t1 + 3t4 + t3*(B + sqrt((B^2*t3^2 - 2B*t3*t5 + 9t5^2)/t3^2)))/(4t3),
(3t1 + 3sqrt(t1)*(-2A + B) + t3*(B - sqrt((B^2*t3^2 - 2B*t3*t5 + 9t5^2)/t3^2)))/(8t3),
(t1 + sqrt(t1)*(-2A + B))/(2t3))

   (10.0, 6.0, 15.0)

 

算額（その259）

算額（その259）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（70）
長野県長野市 久保寺観音堂 享和3年(1803)

鉤股弦（直角三角形）において，
   (1) 鉤^3 + 股^3 = 223.894125
   (2) 中鉤^3 + 短弦^3 = 48.361131
のとき，鉤，股，弦を求めよ。

「中鉤」，「短弦」は図の如し。
以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms 鉤::positive, 股::positive, 弦::positive, 中鉤::positive, 短弦::positive, 長弦::positive;

eq1 = 鉤^3 + 股^3 - 223.894125
eq2 = 中鉤^3 + 短弦^3 - 48.361131
eq3 = 鉤^2 + 股^2 - 弦^2
eq4 = 短弦^2 + 中鉤^2 - 鉤^2
eq5 = 長弦^2 + 中鉤^2 - 股^2
eq6 = 長弦 + 短弦 - 弦

res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (鉤, 股, 弦, 中鉤, 短弦, 長弦))

   1-element Vector{NTuple{6, Sym}}:
    (4.05000000000000, 5.40000000000000, 6.75000000000000, 3.24000000000000, 2.43000000000000, 4.32000000000000)

   鉤 = 4.05000;  股 = 5.40000;  弦 = 6.75000;  中鉤 = 3.24000;  短弦 = 2.43000;  長弦 = 4.32000

鉤 = 4.05, 股 = 5.40, 弦 = 6.75寸 である。

using Printf
using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, 弦, 中鉤, 短弦, 長弦) = res[1]
   @printf("鉤 = %.5f;  股 = %.5f;  弦 = %.5f;  中鉤 = %.5f;  短弦 = %.5f;  長弦 = %.5f\n", 鉤, 股, 弦, 中鉤, 短弦, 長弦)
   plot([股, 股, 0, 股, 長弦*股/弦], [ 0, 鉤, 0, 0, 長弦*鉤/弦], linecolor=:black, linewidth=0.5)
   if more == true
       point(股, 鉤/2, " 鉤", mark=false)
       point(股/2, 0, " 股", :green, :center, :bottom, mark=false)
       point((股 + 長弦*股/弦)/2, 長弦*鉤/弦/2, "  中鉤", mark=false)
       point(長弦*股/弦/2, 長弦*鉤/弦/2, " 長弦", mark=false)
       point((股 + 長弦*股/弦)/2, (鉤 + 長弦*鉤/弦)/2, " 短弦", mark=false)
       point(股/2, 鉤/2, "弦 = 長弦 + 短弦   ", :green, :right, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その258）

算額（その258）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（70）
長野県長野市 久保寺観音堂 享和3年(1803)

菱形の中に大円 1 個と小円 2 個が入っている。菱形の辺の長さの和と菱形の横の対角線の長さの積が 160 平方寸，菱形の辺の長さの和と大円の直径の積が 96 平方寸のとき，小円の直径を求めよ。

菱形の中心を原点とし，長い方と短い方の対角線の長さを 2x, 2y とおく。
大円の半径，中心座標を r1, (0, 0)，
小円の半径，中心座標を r2, (r2, 0) とおく。

以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x::positive, y::positive, sinθ::positive;

xy = sqrt(x^2 + y^2)
sinθ = y / xy
eq1 = 2x * 4xy - 160
eq2 = 2r1 * 4xy - 96
eq3 = x*sinθ - r1
eq4 = (x - r2)sinθ - r2;

# res = solve([eq1, eq2, eq3, eq4], (r1, r2, x, y))

なぜか solve() では解けないので nlsolve() を用いる。

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)
   (r1, r2, x, y) = u
   return [
       8*x*sqrt(x^2 + y^2) - 160,  # eq1
       8*r1*sqrt(x^2 + y^2) - 96,  # eq2
       -r1 + x*y/sqrt(x^2 + y^2),  # eq3
       -r2 + y*(-r2 + x)/sqrt(x^2 + y^2),  # eq4
   ]
end;

iniv = [1.0, 1, 1, 1]
res = nls(H, ini=iniv);
println(res);

   ([2.4, 1.5, 4.0, 2.9999999999999996], true)

小円の直径は 3 寸である。

using Plots
using Printf

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x, y) = res[1]
   @printf("r1 = %.5f;  r2 = %.5f;  x = %.5f;  y = %.5f\n", r1, r2, x, y)
   println("小円の直径 = $(2r2)")
   plot([x, 0, -x, 0, x], [0, y, 0, -y, 0], linecolor=:black, linewidth=0.5)
   circlef(0, 0, r1, :red)
   circlef(r2, 0, r2, :yellow)
   circlef(-r2, 0, r2, :yellow)
   circle(0, 0, r1, :black)
   circle(r2, 0, r2, :black)
   circle(-r2, 0, r2, :black)
   if more == true
       point(0, y, " y", :black, :left, :bottom)
       point(x, 0, " y", :black, :left, :bottom)
       point(0.2, 1.8, "大円:r1", :black, mark=false)
       point(r2, 0, "小円:r2\n", :black, :center, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;