算額（その423）

算額（その423）

長野県諏訪市 諏訪大明神 天保12年(1841)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

直角三角形（鈎股弦）で，その面積が 210 歩，周と弦の和が 99 寸のとき，各辺の長さを求めよ。

鈎，股，弦の長さを 鈎, 股, 弦とし，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive,
     a::positive, b::positive;

(a, b) = (210, 99)
eq1 = 股*鈎/2 - a
eq2 = 鈎 + 股 + 2弦 - b
eq3 = 鈎^2 + 股^2 - 弦^2

res = solve([eq1, eq2, eq3], (鈎, 股, 弦))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    (-b/6 - sqrt(12*a + b^2)/3 - sqrt(-24*a + 5*b^2 + 4*b*sqrt(12*a + b^2))/6, -b/6 + sqrt(-72*a + (b + 2*sqrt(12*a + b^2))^2)/6 - sqrt(12*a + b^2)/3, 2*b/3 + sqrt(12*a + b^2)/3)
    (-b/6 - sqrt(12*a + b^2)/3 + sqrt(-24*a + 5*b^2 + 4*b*sqrt(12*a + b^2))/6, -b/6 - sqrt(-72*a + (b + 2*sqrt(12*a + b^2))^2)/6 - sqrt(12*a + b^2)/3, 2*b/3 + sqrt(12*a + b^2)/3)
    (-b/6 + sqrt(12*a + b^2)/3 - sqrt(-24*a + 5*b^2 - 4*b*sqrt(12*a + b^2))/6, -b/6 + sqrt(-72*a + (b - 2*sqrt(12*a + b^2))^2)/6 + sqrt(12*a + b^2)/3, 2*b/3 - sqrt(12*a + b^2)/3)
    (-b/6 + sqrt(12*a + b^2)/3 + sqrt(-24*a + 5*b^2 - 4*b*sqrt(12*a + b^2))/6, -b/6 - sqrt(-72*a + (b - 2*sqrt(12*a + b^2))^2)/6 + sqrt(12*a + b^2)/3, 2*b/3 - sqrt(12*a + b^2)/3)

鈎 < 股 なので，3 番目の組が適解である。

(a, b) = (210, 99)
(-b/6 - sqrt(12*a + b^2)/3 - sqrt(-24*a + 5*b^2 + 4*b*sqrt(12*a + b^2))/6, -b/6 + sqrt(-72*a + (b + 2*sqrt(12*a + b^2))^2)/6 - sqrt(12*a + b^2)/3, 2*b/3 + sqrt(12*a + b^2)/3) |> println
(-b/6 - sqrt(12*a + b^2)/3 + sqrt(-24*a + 5*b^2 + 4*b*sqrt(12*a + b^2))/6, -b/6 - sqrt(-72*a + (b + 2*sqrt(12*a + b^2))^2)/6 - sqrt(12*a + b^2)/3, 2*b/3 + sqrt(12*a + b^2)/3) |> println
(-b/6 + sqrt(12*a + b^2)/3 - sqrt(-24*a + 5*b^2 - 4*b*sqrt(12*a + b^2))/6, -b/6 + sqrt(-72*a + (b - 2*sqrt(12*a + b^2))^2)/6 + sqrt(12*a + b^2)/3, 2*b/3 - sqrt(12*a + b^2)/3) |> println
(-b/6 + sqrt(12*a + b^2)/3 + sqrt(-24*a + 5*b^2 - 4*b*sqrt(12*a + b^2))/6, -b/6 - sqrt(-72*a + (b - 2*sqrt(12*a + b^2))^2)/6 + sqrt(12*a + b^2)/3, 2*b/3 - sqrt(12*a + b^2)/3) |> println

   (-102.91912585224469, -4.080874147755303, 103.0)
   (-4.080874147755303, -102.91912585224469, 103.0)
   (20.0, 21.0, 29.0)
   (21.0, 20.0, 29.0)

(a, b) = (210, 99) のとき，(鈎, 股, 弦) = (20, 21, 29) である。

 

算額（その422）

算額（その422）

長野県 諏訪社 文化2年(1805)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

鈎股弦に大円と菱形が内接し，菱形に小円が内接している。小円の直径が　83.7 寸のとき，大円の直径を求めよ。

鈎，股，弦の長さを 鈎, 股, 弦
大円の半径と中心座標を r1, (股 - r1, r1)
小円の半径と中心座標を r2, (x2, r2)
菱形の一辺と股の交点座標を (a, 0) とする（図参照）
菱形の高さは 2r2 である（小円は菱形の上辺と下辺に接している）。
以下の連立方程式を立て，数値解を求める。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive,
     r1::positive, r2::positive, x2::positive,
     a::positive, b::positive;

eq1 = 鈎^2 + 股^2 - 弦^2
eq2 = 鈎 + 股 - 弦 - 2r1
eq3 = distance(0, 0, 股, 鈎, x2, r2) - r2^2
eq4 = distance(a, 0, 股, 2r2, x2, r2) - r2^2
eq5 = (股 - r1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq6 = 2r2/(股 - a) - 鈎/股;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (鈎, 股, 弦, r1, x2, a))

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)
   (鈎, 股, 弦, r1, x2, a) = u
   return [
       -弦^2 + 股^2 + 鈎^2,  # eq1
       -2*r1 - 弦 + 股 + 鈎,  # eq2
       -r2^2 + (r2 - 鈎*(r2*鈎 + x2*股)/(股^2 + 鈎^2))^2 + (x2 - 股*(r2*鈎 + x2*股)/(股^2 + 鈎^2))^2,  # eq3
       -r2^2 + (r2 - 2*r2*(a^2 - a*x2 - a*股 + 2*r2^2 + x2*股)/(a^2 - 2*a*股 + 4*r2^2 + 股^2))^2 + (x2 - (2*r2^2*(a - 2*x2 + 股) + x2*(a^2 - 2*a*股 + 4*r2^2 + 股^2))/(a^2 - 2*a*股 + 4*r2^2 + 股^2))^2,  # eq4
       (r1 - r2)^2 - (r1 + r2)^2 + (-r1 - x2 + 股)^2,  # eq5
       2*r2/(-a + 股) - 鈎/股,  # eq6
   ]
end;

r2 = 83.7/2
iniv = [big"180.0", 330.0, 378.0, 72.0, 150.0, 192.0]

res = nls(H, ini=iniv);
println(res);

   (BigFloat[177.7149360644993599352411051684090983213846737911418418853773103721415723113209, 347.4194693522861942938268302046185973814888150527085239226608605207901268261221, 390.2343990288826671178217820524254737950070308416933506657933666693103715275904, 67.4500031939514435556230766603011109539332290010785075711224021118106638049274, 173.709734676143097146913415102309298777291728644346174703749593832642708454533, 183.7922007121731238055940109044496826135959195331646665299783932361428838118555], true)

r2 = 41.85
r1 = 67.45;  x2 = 173.71;  鈎 = 177.715;  股 = 347.419;  弦 = 390.234;  a = 183.792

小円の直径が 83.7 寸のとき，大円の直径は 134.9 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 83.7/2
   (鈎, 股, 弦, r1, x2, a) = res[1]
   @printf("r2 = %g\n", r2)
   @printf("r1 = %g;  x2 = %g;  鈎 = %g;  股 = %g;  弦 = %g;  a = %g\n", r1, x2, 鈎, 股, 弦, a)
   @printf("大円の直径 = %g\n", 2r1)
   plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(股 - r1, r1, r1)
   circle(x2, r2, r2, :blue)
   segment(a, 0, 股, 2r2)
   segment(股 - a, 2r2, 股, 2r2)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(a, 0, " a", :black, :right, :bottom, delta=delta)
       point(股 - a, 2r2, "(股-a,2r2)", :black, :right, :bottom, delta=delta)
       point(股 - r1, r1, "大円:r1,(股-r1,r1)", :red, :center, :bottom, delta=2delta)
       point(x2, r2, "小円:r2,(x2,r2)", :blue, :center, :bottom, delta=2delta)
       point(股, 2r2, " (股,2r2)", :black, :left, :vcenter)
       point(股, 鈎, " (股,鈎)", :black, :left, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;