算額（その1626） - 算額あれこれ

長野県下高井郡木島平村往郷水穂神社寛政12年(1800)

中村信弥「改訂増補長野県の算額」県内の算額(P.62)
http://www.wasan.jp/zoho/zoho.html
キーワード：長方形，菱形
#Julia #SymPy #算額 #和算 #数学

長方形の中に，大小の矢によって菱形ができている。直長が 480 寸，直平が 360 寸，大矢が 400，小矢が 351 寸のとき，菱面を求めよ。

注：直長，直平は長方形の長辺と短辺，菱面は菱形の一辺の長さ。大矢，小矢は図の DE, BC

三角形 ACB において，∠ACB を θ として，第二余弦定理を使う。
三角形 BDE において，∠BDE を 180° - θ として，第二余弦定理を使う。
この 2 つの方程式から，菱面と θ を求める。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy
@syms 大矢::positive, 小矢::positive, 菱面::positive,
直長::integer, 直平::positive, θ::positive
eq1 = (大矢 - 菱面)^2 + 小矢^2 - 2(大矢 - 菱面)*小矢*cos(θ) - 直平^2
eq2 = 大矢^2 + (小矢 - 菱面)^2 - 2大矢*(小矢 - 菱面)*cos(PI - θ) - 直長^2;
eq2 = 大矢^2 + (小矢 - 菱面)^2 + 2大矢*(小矢 - 菱面)*cos(θ) - 直長^2;

方程式自体は簡潔なので，直長，直平，大矢，小矢に実際の数値を与えれば容易に解が得られる。

菱面 = 175, θ = acos(7/25) である。

eq11 = eq1(直長 => 480, 直平 => 360, 大矢 => 400, 小矢 => 351)
eq12 = eq2(直長 => 480, 直平 => 360, 大矢 => 400, 小矢 => 351)
res = solve([eq11, eq12], (菱面, θ))[1]

(175, acos(7/25))

しかし，それらを変数のまま数式解を得ようとすると，得られる解は SymPy で簡約化しきれない複雑な式になる。

@syms x
eq11 = eq1(cos(θ) => x)
eq11 |> println

-x*小矢*(2*大矢 - 2*菱面) + 小矢^2 - 直平^2 + (大矢 - 菱面)^2

ans_x = solve(eq11, x)[1]
ans_x |> println

(小矢^2 - 直平^2 + (大矢 - 菱面)^2)/(2*小矢*(大矢 - 菱面))

eq12 = eq2(cos(θ) => ans_x) |> simplify |> numerator;
eq12 |> println

大矢*(小矢 - 菱面)*(小矢^2 - 直平^2 + (大矢 - 菱面)^2) + 小矢*(大矢 - 菱面)*(大矢^2 - 直長^2 + (小矢 - 菱面)^2)

ans_菱面 = solve(eq12, 菱面)[1];
ans_菱面 |> println

-(-3*(大矢^3 + 3*大矢^2*小矢 + 3*大矢*小矢^2 - 大矢*直平^2 + 小矢^3 - 小矢*直長^2)/(大矢 + 小矢) + (-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)^2/(大矢 + 小矢)^2)/(3*(sqrt(-4*(-3*(大矢^3 + 3*大矢^2*小矢 + 3*大矢*小矢^2 - 大矢*直平^2 + 小矢^3 - 小矢*直長^2)/(大矢 + 小矢) + (-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)^2/(大矢 + 小矢)^2)^3 + (27*(-2*大矢^3*小矢 - 2*大矢*小矢^3 + 大矢*小矢*直平^2 + 大矢*小矢*直長^2)/(大矢 + 小矢) - 9*(-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)*(大矢^3 + 3*大矢^2*小矢 + 3*大矢*小矢^2 - 大矢*直平^2 + 小矢^3 - 小矢*直長^2)/(大矢 + 小矢)^2 + 2*(-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)^3/(大矢 + 小矢)^3)^2)/2 + 27*(-2*大矢^3*小矢 - 2*大矢*小矢^3 + 大矢*小矢*直平^2 + 大矢*小矢*直長^2)/(2*(大矢 + 小矢)) - 9*(-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)*(大矢^3 + 3*大矢^2*小矢 + 3*大矢*小矢^2 - 大矢*直平^2 + 小矢^3 - 小矢*直長^2)/(2*(大矢 + 小矢)^2) + (-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)^3/(大矢 + 小矢)^3)^(1/3)) - (sqrt(-4*(-3*(大矢^3 + 3*大矢^2*小矢 + 3*大矢*小矢^2 - 大矢*直平^2 + 小矢^3 - 小矢*直長^2)/(大矢 + 小矢) + (-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)^2/(大矢 + 小矢)^2)^3 + (27*(-2*大矢^3*小矢 - 2*大矢*小矢^3 + 大矢*小矢*直平^2 + 大矢*小矢*直長^2)/(大矢 + 小矢) - 9*(-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)*(大矢^3 + 3*大矢^2*小矢 + 3*大矢*小矢^2 - 大矢*直平^2 + 小矢^3 - 小矢*直長^2)/(大矢 + 小矢)^2 + 2*(-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)^3/(大矢 + 小矢)^3)^2)/2 + 27*(-2*大矢^3*小矢 - 2*大矢*小矢^3 + 大矢*小矢*直平^2 + 大矢*小矢*直長^2)/(2*(大矢 + 小矢)) - 9*(-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)*(大矢^3 + 3*大矢^2*小矢 + 3*大矢*小矢^2 - 大矢*直平^2 + 小矢^3 - 小矢*直長^2)/(2*(大矢 + 小矢)^2) + (-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)^3/(大矢 + 小矢)^3)^(1/3)/3 - (-2*大矢^2 - 2*大矢*小矢 - 2*小矢^2)/(3*(大矢 + 小矢))

菱面の式に実直を代入すると虚数解が得られるが虚数部は誤差程度なので，実数部を取れば数値解は正確に求まる。

ans_菱面(大矢 => 400, 小矢 => 351, 直長 => 480, 直平 => 360).evalf() |> println

175.0 + 0.e-19*I

数式は長いが以下のようにすれば短い記述で数値解を求めることができる。

function ans2_菱面(直長, 直平, 大矢, 小矢)
s = -2*大矢^3*小矢 - 2*大矢*小矢^3 + 大矢*小矢*直平^2 + 大矢*小矢*直長^2
t = 大矢 + 小矢
u = -2大矢^2 - 2大矢*小矢 - 2小矢^2
v = 大矢^3 + 3大矢^2*小矢 + 3大矢*小矢^2 - 大矢*直平^2 + 小矢^3 - 小矢*直長^2
w = (27s/(2t) + sqrt(Complex(4(3v/t - u^2/t^2)^3 +
(27s/t - 9u*v/t^2 + 2u^3/t^3)^2))/2 - 9u*v/(2*t^2) + u^3/t^3)^(1/3)
abs( (3v/t - u^2/t^2)/(3w) - w/3 - u/(3t))
end
ans2_菱面(480, 360, 400, 351)

175.0000000000001

θ は以下の方程式を解いて求めることができる。

eq = cos(θ) - ans_x
ans_θ = solve(eq, θ)[2];
ans_θ |> println

acos( (大矢^2 - 2*大矢*菱面 + 小矢^2 - 直平^2 + 菱面^2)/(2*小矢*(大矢 - 菱面)))

ans_θ(直長 => 480, 直平 => 360, 大矢 => 400, 小矢 => 351, 菱面 => 175).evalf()*(180/π) |> println

73.7397952916880

∠DBE は以下のようにして求めることができる。

@syms φ
eq3 = (小矢 - 菱面)^2 + 直長^2 - 2(小矢 - 菱面)*直長*cos(φ) - 大矢^2
ans_φ = solve(eq3, φ)[2]
ans_φ |> println

acos( (-大矢^2 + 小矢^2 - 2*小矢*菱面 + 直長^2 + 菱面^2)/(2*直長*(小矢 - 菱面)))

ans_φ(直長 => 480, 直平 => 360, 大矢 => 400, 小矢 => 351, 菱面 => 175).evalf()*(180/π) |> println

53.1301023541560

描画関数プログラムのソースを見る

function draw(直長, 直平, 大矢, 小矢, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
菱面 = ans2_菱面(直長, 直平, 大矢, 小矢)
θ = acos( (大矢^2 - 2*大矢*菱面 + 小矢^2 - 直平^2 + 菱面^2)/(2*小矢*(大矢 - 菱面)))*(180/π)
φ = acos( (-大矢^2 + 小矢^2 - 2*小矢*菱面 + 直長^2 + 菱面^2)/(2*直長*(小矢 - 菱面)))*(180/π)
@printf("直長 = %g; 直平 = %g; 大矢 = %g; 小矢 = %g; 菱面 = %g; θ = %g; φ = %g\n", 直長, 直平, 大矢, 小矢, 菱面, θ, φ)
plot([0, 直長, 直長, 0, 0], [0, 0, 直平, 直平, 0], color=:green, lw=0.5)
(Cx, Cy) = 小矢.*[cosd(φ), sind(φ)]
(Dx, Dy) = (小矢 - 菱面).*[cosd(φ), sind(φ)]
segment(0, 0, Cx, Cy, :blue)
segment(直長, 直平, 直長 - Cx, 直平 - Cy, :blue)
segment(Dx, Dy, 直長, 0, :blue)
segment(直長 - Dx, 直平 - Dy, 0, 直平, :blue)
circle(Cx, Cy, 0.04直長, beginangle=180 - θ + φ, endangle= 180 + φ)
circle(Cx, Cy, 0.045直長, beginangle=180 - θ + φ, endangle= 180 + φ)
circle(0, 0, 0.04直長, beginangle=0, endangle= φ)
circle(0, 0, 0.045直長, beginangle=0, endangle= φ)
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(0, 直平, " A", :green, :left, :bottom, delta=delta/2)
point(0, 0, " B")
point(0, 0, "φ", :red, :left, :bottom, delta=2delta, deltax=5delta, mark=false)
point(Cx, Cy, "C", :green, :left, delta=-delta)
point(Cx, Cy, "θ", :red, :right, :vcenter, deltax=-5delta, mark=false)
point(Dx, Dy, "D", :green, :left, delta=-delta)
point(直長, 0, " E")
point(直長 - Cx, 直平 - Cy, "F", :green, :left, delta=-delta)
point(直長 - Dx, 直平 - Dy, "G", :green, :left, delta=-delta)
ylims!(-5delta, 直平 + 2delta)
end
end;