算額(その511)
長野県長野市安茂里正覚院 久保寺観音堂 文政13年(1830)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
半円内に2本の斜線を引き,甲円 2 個と乙円を入れる。乙円の径が 1 寸のとき,甲円の径はいかほどか。
各円の半径,中心座標を以下のようにする。
外円: R, (0, 0)
甲円: r1, (x1, y1) および (x2, r1)
乙円: r3, (x3, y3)
外円の円周上にある2本の斜線の交点座標を (x, y)
とおき,以下の連立方程式の数値解を求める。
include("julia-source.txt");
using SymPy
@syms x1, y1, r1, x2, r3, x3, y3, x, R
r3 = 1//2
y = sqrt(R^2 - x^2) # (x, y) は外円の円周上にある
eq1 = x1^2 + y1^2 - (R - r1)^2 # 右の甲円が外円に内接する
eq2 = x3^2 + y3^2 - (R - r3)^2 # 乙円が外円に内接する
eq3 = sqrt((R + x)^2 + y^2) + sqrt((R - x)^2 + y^2) - 2R - 2r1 # 鉤股弦に内接する円の直径
eq4 = (y3/x3)*(y/(x + R)) + 1
eq5 = numerator(distance(-R, 0, x, y, x2, r1) - r1^2) |> simplify
eq6 = numerator(distance(R, 0, x, y, x1, y1) - r1^2) |> simplify
eq7 = ((x - R)/2 - x3)^2 + (y/2 - y3)^2 - r3^2
eq8 = (y1/x1) * (y/(x - R)) + 1;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (x1, y1, r1, x2, x3, y3, x, R))
using Printf
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)
(x1, y1, r1, x2, x3, y3, x, R) = u
t = sqrt(R^2 - x^2)
return [
x1^2 + y1^2 - (R - r1)^2, # eq1
x3^2 + y3^2 - (R - 1/2)^2, # eq2
-2*R - 2*r1 + sqrt(R^2 - x^2 + (R - x)^2) + sqrt(R^2 - x^2 + (R + x)^2), # eq3
1 + y3*sqrt(R^2 - x^2)/(x3*(R + x)), # eq4
(-4*R^2*r1^2 + (R*r1 - R*t + r1*x - x2*t)^2 + (R^2 - R*x + R*x2 - r1*t - x*x2)^2)/(4*R^2), # eq5
(-4*R^2*r1^2 + (R*y1 - R*t - x*y1 + x1*t)^2 + (-R^2 - R*x + R*x1 + x*x1 + y1*t)^2)/(4*R^2), # eq6
(-y3 + t/2)^2 + (-R/2 + x/2 - x3)^2 - 1/4, # eq7
1 + y1*t/(x1*(-R + x)), # eq8
]
end;
iniv = BigFloat[1.7, 4.2, 2.1, -3.5, -5.8, 2.4, -4.6, 6.4]
res = nls(H, ini=iniv)
println(res)
(BigFloat[3.461538461538461538461538509169459561367959236632964526959234465030510845593508, 8.307692307692307692307692429999766197094115669634049393814154165178838427207838, 4.000000000000000000000000058695639811160390953981958412762403947086964336791015, -7.000000000000000000000000103675105612121327124114588925459978536567327493943071, -11.53846153846153846153846171778119379511106203957672690108302039508364569817559, 4.807692307692307692307692375277671136103845336879392455246078657573883943485054, -9.153846153846153846153846290689476553117442078831806765600936112963486841120747, 13.00000000000000000000000019241303418103091110455672021661459202029874280726291], true)
r3 = 0.5; x1 = 3.46154; y1 = 8.30769; r1 = 4; x2 = -7; x3 = -11.5385; y3 = 4.80769; x = -9.15385; y = 9.23077; R = 13
甲円の径 = 2r1 = 8
甲円の半径は 4 である。元の単位では甲円の直径は 8 寸である。
using Plots
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r3 = 1//2
(x1, y1, r1, x2, x3, y3, x, R) = res[1]
y = sqrt(R^2 - x^2)
@printf("r3 = %g; x1 = %g; y1 = %g; r1 = %g; x2 = %g; x3 = %g; y3 = %g; x = %g; y = %g; R = %g\n",
r3, x1, y1, r1, x2, x3, y3, x, y, R)
@printf("甲円の径 = 2r1 = %g\n", 2r1)
plot([R, x, -R, R], [0, y, 0, 0], color=:brown, lw=0.5)
circle(0, 0, R, beginangle=0, endangle=180)
circle(x1, y1, r1, :blue)
circle(x2, r1, r1, :green)
circle(x3, y3, r3, :magenta)
if more == true
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
point(x1, y1, " 甲円:r1\n (x1,y1)", :blue, :left, :vcenter)
point(x2, r1, " 甲円:r1\n (x2,r1)", :blue, :left, :top, delta=-delta)
point(x3, y3, " 乙円:r3,(x3,y3)", :magenta, :left, :vcenter)
point(x, y, "(x,y) ", :brown, :right, :vcenter)
else
plot!(showaxis=false)
end
end;
※コメント投稿者のブログIDはブログ作成者のみに通知されます