算額（その443）

算額（その443）

天明五年乙巳秋九月(1785) 藤田貞資門人 筑後州久留米藩 清水與市道香

藤田貞資(1789)：神壁算法巻下
http://www.wasan.jp/jinpeki/jinpekisanpo2.pdf

求めるものと条件として与えられる数値が異なるだけで，算額（その31）と基本的に同じである。

図のように，直線の上に大円，中円，小円，甲円，乙円がそれぞれ外接しあっている。甲円，乙円の径を98寸1厘，121寸 としたとき，中円の径を求めよ。

それぞれの円の中心座標と半径を (0, r1, r1), (x2, r2, r2), (x3, r3, r3), (x4, r4, r4), (x5, r5, r5) とする。

include("julia-source.txt")

using SymPy

@syms r1, r2, r3::positive, r4::positive, r5::positive
@syms x2::positive, x3::positive, x4::positive, x5::positive;

(r4, r5) = (981, 1210) .// 20;
eq1 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2;
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2;
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2;
eq4 = x5^2 + (r1 - r5)^2 - (r1 + r5)^2;
eq5 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2;
eq6 = (x2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2;
eq7 = (x5 - x3)^2 + (r3 - r5)^2 - (r3 + r5)^2;

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r2, r3, x2, x3, x4, x5))

   2-element Vector{NTuple{7, Sym}}:
    (3917133*sqrt(1090)/3682898 + 304705467/7365796, 121/2, 3917133*sqrt(1090)/52441 + 260073891/104882, 3993*sqrt(1090)/2714 + 118701/1357, 188259786/310753 + 11751399*sqrt(1090)/621506, 118701/1357 + 32373*sqrt(1090)/6785, 3993*sqrt(1090)/2714 + 118701/1357)
    (35254197*sqrt(1090)/3682898 + 2742349203/7365796, 70508394*sqrt(1090)/14891881 + 6218626689/29783762, 2340665019/104882 - 35254197*sqrt(1090)/52441, 2340665019/5236663 + 176270985*sqrt(1090)/10473326, -401684184/310753 + 35254197*sqrt(1090)/621506, 356103/1357 + 97119*sqrt(1090)/6785, 11979*sqrt(1090)/2714 + 356103/1357)

2番めのものが適解である。

name = ["r1", "r2", "r3", "x2", "x3", "x4", "x5"]
j = 2
for i in 1:7
   println("$(name[i]) = $(res[j][i].evalf())")
end

   r1 = 688.343020749222
   r2 = 365.108908025960
   r3 = 122.232157448001
   x2 = 1002.63686078867
   x3 = 580.129821644955
   x4 = 734.990886123080
   x5 = 408.140920542540

中円の直径は 730.21781605192 である。算額の答えでは 729 となっており，かなり違いがある。

2res[2][2].evalf()

   730.21781605192

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r4, r5) = (981, 1210) ./ 20
   r1, r2, r3, x2, x3, x4, x5 = (35254197*sqrt(1090)/3682898 + 2742349203/7365796, 70508394*sqrt(1090)/14891881 + 6218626689/29783762, 2340665019/104882 - 35254197*sqrt(1090)/52441, 2340665019/5236663 + 176270985*sqrt(1090)/10473326, -401684184/310753 + 35254197*sqrt(1090)/621506, 356103/1357 + 97119*sqrt(1090)/6785, 11979*sqrt(1090)/2714 + 356103/1357)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  r4 = %g;  r5 = %g\n", r1, r2, r3, r4, r5)
   @printf("中円の直径 = %g\n", 2r2)
   plot()
   circle(0, r1, r1, :red)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :brown)
   circle(x4, r4, r4, :green)
   circle(x5, r5, r5, :red)
   hline!([0], color=:black, linewidth=0.25)
   if more
       point(0, r1, "大円 r1")
       point(x2, r2,"中円 r2")
       point(x3, r3,"小")
       point(x4, r4,"乙")
       point(x5, r5,"甲")
   end
end;