藤田貞資：精要算法（下巻） 天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
キーワード：円5個，菱形
#Julia #SymPy #算額 #和算 #数学
菱形の中に大円，中円，小円が入っている。小円と中円の直径の積が 1 寸（平方寸），双股（小円と中円が菱形と接する 2 点間の距離）が 5 寸のとき，大円の直径はいかほどか。
菱形の対角線を \(2a,\ 2b\)
大円の半径と中心座標を \(r_1,\ (0,\ 0)\)
中円の半径と中心座標を \(r_2,\ (0,\ r_1 + r_2)\)
小円の半径と中心座標を \(r_3,\ (r_1 + r_3,\ 0)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
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θ;
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
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;
以下のアイコンをクリックして応援してください
