80 岩手県藤沢町保呂羽字二本柳 保呂羽神社 明治26年(1893)
安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円3個，直線上
#Julia #SymPy #算額 #和算 #数学
直線の上に甲円，乙円が互いに接して載っている。その 2 つの円の上に丙円が載っている。3 個の円の直径が与えられたとき，丙円のてっぺんまでの高さを求めよ。
甲円の半径と中心座標を \(r_1, (x_1, r_1)\)
乙円の半径と中心座標を \(r_2, (0, r_2)\)
丙円の半径と中心座標を \(r_3, (x_3, y_3)\)
とおき，以下の連立方程式を解く。
高さは \(y_3 + r_3\) である。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r1::positive, x1::positive, r2::positive,
r3::positive, x3::positive, y3
eq1 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand
eq2 = (x1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2 |> expand
eq3 = x3^2 + (y3 - r2)^2 - (r2 + r3)^2 |> expand
#res = solve([eq1, eq2, eq3], (x1, x3, y3))[1]; # 1 of 2
SymPy で解くと，\(x_1\) の式が信じられないくらい長くなってしまうので，別途解く。
\(x_1 = 2\sqrt{r_1 r_2}\) である。ちなみに，「算法助術」の公式40ではある。
ans_x1 = solve(eq1, x1)[1]
ans_x1 |> println
2*sqrt(r1)*sqrt(r2)
\(x_3, y_3\) も別々に解いたほうが若干短くなるが，連立方程式として解く。eq2 は \(x_1\) を含むので，代入しておく。
eq2 = eq2(x1 => ans_x1);
(ans_x3, ans_y3) = solve([eq2, eq3], (x3, y3))[2]; # 2 of 2
ans_x3 = ans_x3 |> factor
ans_x3 |> println
2*sqrt(r1)*sqrt(r2)*(r1*r2 - r1*sqrt(r3)*sqrt(r1 + r2 + r3) - r1*r3 + r2^2 + r2*sqrt(r3)*sqrt(r1 + r2 + r3) + r2*r3)/(r1 + r2)^2
ans_y3 = ans_y3 |> factor
ans_y3 |> println
(2*r1^2*r2 - r1^2*r3 + 2*r1*r2^2 + 4*r1*r2*sqrt(r3)*sqrt(r1 + r2 + r3) + 2*r1*r2*r3 - r2^2*r3)/(r1 + r2)^2
ans_x1(r1 => 5/2, r2 => 3/2, r3 => 2/2).evalf() |> println
ans_x3(r1 => 5/2, r2 => 3/2, r3 => 2/2).evalf() |> println
ans_y3(r1 => 5/2, r2 => 3/2, r3 => 2/2).evalf() |> println
3.87298334620742
0.669041418324543
3.90881372890605
高さは \(ans_{y3} + r_3 = 3.90881372890605 + 1 = 4.9088137289060505\) である。
高さだけを求めるならば，「算法助術」の公式47を適用する。
@syms h
公式47 = r1*h + r2*h - 2r1*r2 - 2sqrt(2r1*r2*r3*h)
res = solve(公式47, h)[2]
res |> println
2*r1*r2*(r1 + r2 + 2*sqrt(r3)*sqrt(r1 + r2 + r3) + 2*r3)/(r1^2 + 2*r1*r2 + r2^2)
res(r1 => 5/2, r2 => 3/2, r3 => 2/2).evalf() |> println
4.90881372890605
「術」も，簡約化した結果の見た目は異なるが，公式47と等価である。
@syms r1, r2, r3
天 = 2(r1 + r2 + r3)
A = sqrt(天 * 2r3)*2
B = 天 - A + 2r3
高 = 2r1*2r2/B |> simplify
高 |> println
2*r1*r2/(r1 + r2 + 2*r3 - 2*sqrt(r3*(r1 + r2 + r3)))
高(r1 => 5/2, r2 => 3/2, r3 => 2/2).evalf() |> println
4.90881372890605
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3) = (5, 3, 2) ./ 2
x1 = 2*sqrt(r1)*sqrt(r2)
x3 = 2*sqrt(r1)*sqrt(r2)*(r1*r2 - r1*sqrt(r3)*sqrt(r1 + r2 + r3) - r1*r3 + r2^2 + r2*sqrt(r3)*sqrt(r1 + r2 + r3) + r2*r3)/(r1 + r2)^2
y3 = (2*r1^2*r2 - r1^2*r3 + 2*r1*r2^2 + 4*r1*r2*sqrt(r3)*sqrt(r1 + r2 + r3) + 2*r1*r2*r3 - r2^2*r3)/(r1 + r2)^2
@printf("甲円，乙円，丙円の直径がそれぞれ %g, %g, %g のとき，盤面から丙円のてっぺんまでの高さは %g である。\n", 2r1, 2r2, 2r3, y3 + r3)
@printf("x1 = %g; x3 = %g; y3 = %g\n",x1, x3, y3)
plot()
circle(x1, r1, r1)
circle(0, r2, r2, :green)
circle(x3, y3, r3, :blue)
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)
point(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=delta)
point(0, r2, "乙円:r2,(0,r2)", :green, :center, delta=delta)
point(x3, y3, "丙円:r3,(x3,y3)", :blue, :center, delta=delta)
point(x3, y3 + r3, "高さ=y3+r3", :magenta, :center, :bottom, delta=delta)
end
end;
以下のアイコンをクリックして応援してください
