長野県松本市筑摩 筑摩神社 明治6年(1873)
中村信弥「改訂増補 長野県の算額」県内の算額(P.240)
http://www.wasan.jp/zoho/zoho.html
キーワード：円2個，直角三角形，斜線
#Julia #SymPy #算額 #和算 #数学
鈎股弦（直角三角形）に 2 本の斜線を容れ，区切られた領域に 2 個の等円を容れる。鈎が 5 寸のとき，等円の直径を求めよ。
斜線と底辺，斜辺との交点座標を \( (a,\ 0),\ (a + b,\ 0),\ (x,\ y)\)
等円の半径と中心座標を \(r,\ (r,\ y_1),\ (x_1,\ r)\)
とおき，連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, b::positive, r::positive, x1::positive, y1::positive,
x::positive, y::positive;
eq1 = distance(0, 5, a, 0, r, y1) - r^2
eq2 = distance(0, 5, a, 0, x1, r) - r^2
eq3 = distance(0, 0, x, y, r, y1) - r^2
eq4 = distance(0, 0, x, y, x1, r) - r^2
eq5 = distance(0, 5, a + b, 0, x1, r) - r^2
eq6 = (5 - y)*(a + b) - 5x
eq7 = x*y + x*(5-y)/2 + (a + b - x)*y/2 - 5(a + b)/2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])
nlsove() で解く。
println(eq1, ", # eq1")
println(eq2, ", # eq2")
println(eq3, ", # eq3")
println(eq4, ", # eq4")
println(eq5, ", # eq5")
println(eq6, ", # eq6")
println(eq7, ", # eq7")
-r^2 + (y1 - 5*(a^2 - a*r + 5*y1)/(a^2 + 25))^2 + (-a*(a*r - 5*y1 + 25)/(a^2 + 25) + r)^2, # eq1
-r^2 + (r - 5*(a^2 - a*x1 + 5*r)/(a^2 + 25))^2 + (-a*(a*x1 - 5*r + 25)/(a^2 + 25) + x1)^2, # eq2
-r^2 + (r - x*(r*x + y*y1)/(x^2 + y^2))^2 + (-y*(r*x + y*y1)/(x^2 + y^2) + y1)^2, # eq3
-r^2 + (r - y*(r*y + x*x1)/(x^2 + y^2))^2 + (-x*(r*y + x*x1)/(x^2 + y^2) + x1)^2, # eq4
-r^2 + (r - 5*(a^2 + 2*a*b - a*x1 + b^2 - b*x1 + 5*r)/(a^2 + 2*a*b + b^2 + 25))^2 + (x1 - (a + b)*(a*x1 + b*x1 - 5*r + 25)/(a^2 + 2*a*b + b^2 + 25))^2, # eq5
-5*x + (5 - y)*(a + b), # eq6
-5*a/2 - 5*b/2 + x*y + x*(5 - y)/2 + y*(a + b - x)/2, # eq7
function H(u)
(a, b, r, x1, y1, x, y) = u
return [
-r^2 + (y1 - 5*(a^2 - a*r + 5*y1)/(a^2 + 25))^2 + (-a*(a*r - 5*y1 + 25)/(a^2 + 25) + r)^2, # eq1
-r^2 + (r - 5*(a^2 - a*x1 + 5*r)/(a^2 + 25))^2 + (-a*(a*x1 - 5*r + 25)/(a^2 + 25) + x1)^2, # eq2
-r^2 + (r - x*(r*x + y*y1)/(x^2 + y^2))^2 + (-y*(r*x + y*y1)/(x^2 + y^2) + y1)^2, # eq3
-r^2 + (r - y*(r*y + x*x1)/(x^2 + y^2))^2 + (-x*(r*y + x*x1)/(x^2 + y^2) + x1)^2, # eq4
-r^2 + (r - 5*(a^2 + 2*a*b - a*x1 + b^2 - b*x1 + 5*r)/(a^2 + 2*a*b + b^2 + 25))^2 + (x1 - (a + b)*(a*x1 + b*x1 - 5*r + 25)/(a^2 + 2*a*b + b^2 + 25))^2, # eq5
-5*x + (5 - y)*(a + b), # eq6
-5*a/2 - 5*b/2 + x*y + x*(5 - y)/2 + y*(a + b - x)/2, # eq7
]
end;
iniv = [big"2.8", 2.9, 0.9, 3.3, 1.6, 3.4, 2.1]
res = nls(H, ini=iniv);
println(res);
([2.800221924622119, 2.914870271332841, 0.8974814920915777, 3.326236908516906, 1.560761971632354, 3.4319778903414297, 1.9974431096925063], true)
\(a = 2.800222; b = 2.914870; r = 0.897481; x_1 = 3.326237\)
\(y_1 = 1.560762; x = 3.431978; y = 1.997443\)
\(2r = 1.794963\)
等円の直径は 1.794963 寸である。
算額の術では 1.830127 寸としているが，図に描く限り上の方が正しそう。
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b, r, x1, y1, x, y) = res[1]
@printf("a = %.6f; b = %.6f; r = %.6f; x1 = %.6f; y1 = %.6f; x = %.6f; y = %.6f\n", a, b, r, x1, y1, x, y)
@printf("2r = %.6f\n", 2r)
plot([0, a + b, 0, 0, x], [0, 0, 5, 0, y], color=:black, lw=0.5)
segment(0, 5, a, 0)
circle(r, y1, r, :blue)
circle(x1, r, r, :blue)
if more
point(x, y, " (x,y)", :green, :left, :bottom)
point(x1, r, " (x1,r)")
point(r, y1, " (r,y1)")
point(a, 0, "a ", :green, :right, :bottom)
point(a + b, 0, " a+b", :green, :left, :bottom)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
end
end;
以下のアイコンをクリックして応援してください
