裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

算額(その31)

2022年11月19日 | Julia

算額(その31)

大阪府茨木市 井於神社 弘化3年(1846)
http://www.wasan.jp/osaka/iyo.html

(3) 奈良県大和郡山市小泉町 耳成山口神社 嘉永7年(1854)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

キーワード:円5個,直線上

図のように,大円,中円とそれらに接する 3 個の円(甲,乙,丙)がある。大円,中円の径を三寸五分,二寸 としたとき,甲,乙,丙の円の径を求めよ。

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

using SymPy

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

r1 = 35//10;
r2 = 2;
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 = x4^2 + (r1 - r4)^2 - (r1 + r4)^2;
eq5 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2;
eq6 = (x2 - x5)^2 + (r2 - r5)^2 - (r2 + r5)^2;
eq7 = (x5 - x3)^2 + (r3 - r5)^2 - (r3 + r5)^2;

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

   3-element Vector{NTuple{7, Sym}}:
    (154/9 - 56*sqrt(7)/9, -505064*sqrt(7)*sqrt(296 - 110*sqrt(7))/177147 - 1326104*sqrt(296 - 110*sqrt(7))/177147 - 56*sqrt(7)/81 + 242/81 + 513812*sqrt(2072 - 770*sqrt(7))/177147 + 195692*sqrt(7)*sqrt(2072 - 770*sqrt(7))/177147, 7/9 - 7*sqrt(7)/36, 2*sqrt(7), -28/3 + 14*sqrt(7)/3, 2*sqrt(7)/9 + 4*sqrt(14)*sqrt(148 - 55*sqrt(7))/(9*(-1 + sqrt(7))) + 28/9, -7/3 + 7*sqrt(7)/3)
    (154/9 - 56*sqrt(7)/9, -195692*sqrt(7)*sqrt(2072 - 770*sqrt(7))/177147 - 513812*sqrt(2072 - 770*sqrt(7))/177147 - 56*sqrt(7)/81 + 242/81 + 1326104*sqrt(296 - 110*sqrt(7))/177147 + 505064*sqrt(7)*sqrt(296 - 110*sqrt(7))/177147, 7/9 - 7*sqrt(7)/36, 2*sqrt(7), -28/3 + 14*sqrt(7)/3, -4*sqrt(14)*sqrt(148 - 55*sqrt(7))/(9*(-1 + sqrt(7))) + 2*sqrt(7)/9 + 28/9, -7/3 + 7*sqrt(7)/3)
    (56*sqrt(7)/9 + 154/9, -195692*sqrt(7)*sqrt(770*sqrt(7) + 2072)/177147 - 505064*sqrt(7)*sqrt(110*sqrt(7) + 296)/177147 + 56*sqrt(7)/81 + 242/81 + 1326104*sqrt(110*sqrt(7) + 296)/177147 + 513812*sqrt(770*sqrt(7) + 2072)/177147, 7*sqrt(7)/36 + 7/9, 2*sqrt(7), 28/3 + 14*sqrt(7)/3, -28/9 + 2*sqrt(7)/9 + 4*sqrt(14)*sqrt(55*sqrt(7) + 148)/(9*(1 + sqrt(7))), 7/3 + 7*sqrt(7)/3)

題意を満たすのは res[2] である。それにしても,r4 の式はものすごい。

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

   r3 = 154/9 - 56*sqrt(7)/9 = 0.648658508931436
   r4 = -195692*sqrt(7)*sqrt(2072 - 770*sqrt(7))/177147 - 513812*sqrt(2072 - 770*sqrt(7))/177147 - 56*sqrt(7)/81 + 242/81 + 1326104*sqrt(296 - 110*sqrt(7))/177147 + 505064*sqrt(7)*sqrt(296 - 110*sqrt(7))/177147 = 0.316985841490936
   r5 = 7/9 - 7*sqrt(7)/36 = 0.263326133959663
   x2 = 2*sqrt(7) = 5.29150262212918
   x3 = -28/3 + 14*sqrt(7)/3 = 3.01350611830142
   x4 = -4*sqrt(14)*sqrt(148 - 55*sqrt(7))/(9*(-1 + sqrt(7))) + 2*sqrt(7)/9 + 28/9 = 2.10660907167730
   x5 = -7/3 + 7*sqrt(7)/3 = 3.84008639248404

算額の解答では,
甲: 0.625,
乙: 0.36525,
丙: 0.19315 
となっているが,かなりの誤差がある。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1, r2 = 35/10, 2
   r3, r4, r5, x2, x3, x4, x5 = (154/9 - 56*sqrt(7)/9, -195692*sqrt(7)*sqrt(2072 - 770*sqrt(7))/177147 - 513812*sqrt(2072 - 770*sqrt(7))/177147 - 56*sqrt(7)/81 + 242/81 + 1326104*sqrt(296 - 110*sqrt(7))/177147 + 505064*sqrt(7)*sqrt(296 - 110*sqrt(7))/177147, 7/9 - 7*sqrt(7)/36, 2*sqrt(7), -28/3 + 14*sqrt(7)/3, -4*sqrt(14)*sqrt(148 - 55*sqrt(7))/(9*(-1 + sqrt(7))) + 2*sqrt(7)/9 + 28/9, -7/3 + 7*sqrt(7)/3)
   println("r1 = $r1;  r2 = $r2;  r3 = $r3;  r4 = $r4;  r5 = $r5")
   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;

なお,フォードの円デカルトの円定理を簡約化したものであるが,それを順次用いることで円の径を求めることができる。しかし,この方法では円の中心座標についての情報は得られないので,図を描くためには別途中心座標を求める必要がある。

using SymPy

@syms r1, r2, r3::positive, r4::positive, r5::positive
r1 = 35//10
r2 = 2
eq1 = 1/sqrt(r1) + 1/sqrt(r2) - 1/sqrt(r3);
r3 = solve(eq1)[1]
    0.648658508931432

eq2 = 1/sqrt(r1) + 1/sqrt(r3) - 1/sqrt(r4);
r4 = solve(eq2)[1]
    0.316985841490935

eq3 = 1/sqrt(r2) + 1/sqrt(r3) - 1/sqrt(r5);
r5 = solve(eq3)[1]
    0.263326133959661

デカルトの円定理では「曲率=円の径の逆数」を用いる。

using SymPy
@syms r1, r2, r3::positive, r4::positive, r5::positive
@syms k1, k2, k3, k4, k5
r1 = 35//10
r2 = 2
k1, k2, k3, k4, k5 = 1 ./ [r1, r2, r3, r4, r5]

まず r1, r2 から r3 を求める。

k3 =  k1 + k2 + 2√(k1*k2)
1/k3.evalf()
    0.648658508931436

ついで,r1, r3 から r4 を求める。

k4 = k1 + k3 + 2√(k1*k3)
1/k4.evalf()
    0.316985841490937

更に r2,r3 から r5 を求める。

k5 = k2 + k3 + 2√(k2*k3)
1/k5.evalf()
    0.263326133959663

 

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その30)

2022年11月19日 | Julia

算額(その30)

奈良県奈良市 円満寺 天保 15 年
http://www.wasan.jp/nara/enman.html

奈良県奈良市下山町 円満寺 天保15年(1844)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

図のように,1 個の大円と,6 個の小円がある。それぞれの径を求めよ。

見ただけでわかる。小円の径は大円の径の 1/3 である。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1, r2 = 1, 1/3
   println("r1 = $r1;  r2 = $r2")
   plot()
   circle(0, 0, r1)
   for degree = 30:60:330
       circle(-4r1/3*cosd(degree), 4r1/3*sind(degree), r2, :blue)
   end
   plot!([0, 2r1*cosd(210), 2r1*cosd(330), 0], [2r1, 2r1*sind(210), 2r1*sind(330), 2r1], color=:black, linewidth=0.25)
   plot!([0, 2r1*cosd(210), 2r1*cosd(330), 0], [-2r1, -2r1*sind(210), -2r1*sind(330), -2r1], color=:black, linewidth=0.25)
end;

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その29)

2022年11月19日 | Julia

算額(その29)

京都府東山区 安井金比羅宮 平成 7 年
http://www.wasan.jp/kyoto/yasuikonpira3.html

(2) 大阪府豊中市服部元町 服部天神社 天保14年(1843)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

図のように,円の中に大円が 5 個,小円が 1 個ある。それぞれの円の径を求めよ。

外円のないものは「算額(その682)

服部天神社のものは「答」の誤差が大きすぎる。

外円の半径を 1 とし,以下のように記号を定め,方程式を立てる。

using SymPy

@syms r1::positive, r2::positive, x2::positive, y2::positive, k, m;

以下の 4 方程式を立て,解く。

eq1 = x2 - (r1 + r2)cos(PI/10);
eq2 = y2 - (r1 + r2)sin(PI/10);
eq3 = r1 + 2r2 - 1;
eq4 = x2^2 + (r1 + r2 - y2)^2 - 4r2^2;

res = solve([eq1, eq2, eq3, eq4], (r1, r2, x2, y2));

name = ["r1", "r2", "x2", "y2"]
for i in 1:4
   println("$(name[i]) = $(simplify(res[1][i])) = $(res[1][i].evalf())")
end

   r1 = -4*sqrt(5) - 3*sqrt(10 - 2*sqrt(5)) + sqrt(50 - 10*sqrt(5)) + 11 = 0.259616183682500
   r2 = -5 - sqrt(50 - 10*sqrt(5))/2 + 3*sqrt(10 - 2*sqrt(5))/2 + 2*sqrt(5) = 0.370191908158750
   x2 = -sqrt(10*sqrt(5) + 50)/2 - 3*sqrt(5)/2 + 5/2 + 3*sqrt(2*sqrt(5) + 10)/2 = 0.598983089761037
   y2 = -4 - sqrt(5)*sqrt(10 - 2*sqrt(5))/2 + sqrt(10 - 2*sqrt(5)) + 2*sqrt(5) = 0.194621403573804

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = -4*sqrt(5) - 3*sqrt(10 - 2*sqrt(5)) + sqrt(50 - 10*sqrt(5)) + 11
   r2 = -5 - sqrt(50 - 10*sqrt(5))/2 + 3*sqrt(10 - 2*sqrt(5))/2 + 2*sqrt(5)
   x2 = -sqrt(10*sqrt(5) + 50)/2 - 3*sqrt(5)/2 + 5/2 + 3*sqrt(2*sqrt(5) + 10)/2
   y2 = -4 - sqrt(5)*sqrt(10 - 2*sqrt(5))/2 + sqrt(10 - 2*sqrt(5)) + 2*sqrt(5)
   println("r1 = $r1;  r2 = $r2")
   println("x2 = $x2;  y2 = $y2")
   plot()
   r0 = r1 + r2
   circle(0, 0, r1 + 2r2, :black)
   circle(0, 0, r1, :green)
   circle(0, r0, r2)
   circle(r0*cos(pi/10), r0*sin(pi/10), r2)
   circle(-r0*cos(pi/10), r0*sin(pi/10), r2)
   circle(r0*cos(17pi/10), r0*sin(17pi/10), r2)
   circle(-r0*cos(17pi/10), r0*sin(17pi/10), r2)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, r1 + r2, "r1+r2 ", :red, :right)
       point(0, r1, "r1 ", :green, :right, :bottom)
       point(x2, y2, "(x2,y2)", :magenta, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その28)

2022年11月19日 | Julia

算額(その28)

(5) 京都府宮津市天橋立文殊 知恩寺文殊堂 文政元年(1818)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

京都府宮津市天橋立 知恩寺文殊堂 文政元年
http://www.wasan.jp/kyoto/tionji1.html

キーワード:円7個

図のように,4 個の大円と,3 個の小円がある。それぞれの径を求めよ。

見ただけでわかる。大円の径は小円の径の 2 倍である。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1, r2 = 2, 1
   println("r1 = $r1;  r2 = $r2")
   plot()
   circle(0, 0, r1)
   circle(0, r1, r1)
   circle( r1*cos(pi/6), -r1*sin(pi/6), r1)
   circle(-r1*cos(pi/6), -r1*sin(pi/6), r1)
   circle(0, r1 + r2, r2, :blue)
   circle((r1 + r2)*cos(pi/6), -(r1 + r2)*sin(pi/6), r2, :blue)
   circle(-(r1 + r2)*cos(pi/6), -(r1 + r2)*sin(pi/6), r2, :blue)
end;

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村