裏 RjpWiki

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

算額(その72)

2022年12月24日 | Julia

算額(その72)

四二 浦和市西堀(現さいたま市桜区西堀) 氷川神社 嘉永5年(1852)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

埼玉県浦和市西堀 氷川神社 嘉永5年
http://www.wasan.jp/saitama/uhikawa.html

正方形の中に大円,中円,小円が互いに接して入っている。中円の径が 10 寸のとき,小円の径はいかほどか。

下図のように,大円,中円,小円の半径と中心座標をそれぞれ {1, (0, 1)}, {r1, (1-r1, r1)}, {r2,(x2, r2)} とし,方程式を立てて解く。

using SymPy
@syms r1::positive, x2::positive, r2::positive
eq1 = (1 - r1)^2 + (1 - r1)^2 - (1 + r1)^2;
eq2 = x2^2 + (1 - r2)^2 - (1 + r2)^2;
eq3 = (1 - r1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2;
results = solve([eq1, eq2, eq3], (r1, r2, x2));
results[2]

   (3 - 2*sqrt(2), 3/2 - sqrt(2), 2 - sqrt(2))

小円の径は中円の径の r2/r1

(3/2 - sqrt(2)) / (3 - 2*sqrt(2))  # r2/r1

   0.5

よって、中円の径が 10 寸ならば、小円の径は 5 寸である。

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; fontsize=10, mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, fontsize, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="")#, fontfamily="IPAMincho")
    (r1, r2, x2) = (3 - 2*sqrt(2), 3/2 - sqrt(2), 2 - sqrt(2))
    plot([-1, 1, 1, -1, -1], [0, 0, 2, 2, 0], color=:black, lw=0.5)
    circle(0, 1, 1)
    circle(1 - r1, r1, r1, :green)
    circle(x2, r2, r2, :blue)
    if more
       println("r1 = $r1")
       println("r2 = $r2")
       println("x2 = $x2)")
       println("r2/r1 = $(r2/r1)")
        point(0, 1, "1 ", :red, :right)
        point(1 - r1, r1, "\n(1-r1,r1)", :green, :center) 
        point(x2, r2, "(x2,r2)", :blue, :center) 
       vline!([0], color=:black, lw=0.5, xlims=(-0.1, 1.1), ylims=(-0.1, 1.1))
    end
end;

   r1 = 0.1715728752538097
   r2 = 0.08578643762690485
   x2 = 0.5857864376269049)
   r2/r1 = 0.5

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

算額(その71)

2022年12月22日 | Julia

算額(その71)

福島県田村郡三春町山中 田村大元神社 明治34年(1901)5月
http://www.wasan.jp/fukusima/tamuradaigen.html

矩形内に対角線に接して甲円 2 個,乙円 1 個がある。甲円の径が 3 寸のとき,乙円の径はいかほどか。

図のように記号を定め,方程式を解く。

using SymPy
@syms r1::positive, r2::positive;

算額では矩形を「直」と書いているが,これは,鉤 = 3k,股 = 4k,弦 = 5k の直角三角形を 2 つ組み合わせたものである。

大円の半径 r1 は「鉤股弦 ABC と内接円の定理」により,r1*(3 + 4 - 5)/2 = 3 を解けばよい。

r1 = solve(r1*(3 + 4 - 5)//2 - 3)[1]

    3

小円の半径 r2 は,大円と小円が外接することを用いて (x1//2-r1)^2 + (y1-r2-r1)^2 = (r1 + r2)^2 を解けばよい。(x1, y1) は (12, 9) である。

(x1, y1) = (12, 9)
r2 = solve((x1//2-r1)^2 + (y1-r2-r1)^2 - (r1 + r2)^2)[1]

    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 segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="")#, fontfamily="IPAMincho")
   # 矩形の左下角を原点,右上角の座標 (x1, y1)として図を描く
    (x1, y1) = (12, 9)
    (r1, r2) = (3, 2)
    plot([0, x1, x1, 0, 0], [0, 0, y1, y1, 0], color=:black, showaxis=false)
    segment(0, 0, x1, y1)
    segment(0, y1, x1, 0)
    circle(r1, r1, r1)
    circle(x1-r1, r1, r1)
    circle(x1/2, y1-r2, r2, :blue)
   if more
       println("r1 = $r1;  r2 = $r2")
       plot!(xlims=(-0.5, x1+0.5), ylims=(-0.5, y1+0.5), showaxis=true)
       point(0, y1, "A ", :black, :right)
       point(0, 0, "B ", :black, :right)
       point(x1, 0, "C ", :black, :right)
       point(r1, r1, "(r1, r1)", :red, :center)
       point(x1/2, y1-r2, "(x1/2, y1-r2)", :blue, :center)
   end
end;

   r1 = 3;  r2 = 2

 

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

算額(その70)

2022年12月22日 | Julia

算額(その70)

東京都府中市 大國魂神社
http://www2.ttcn.ne.jp/~nagai/sangaku/sangakumondai1.htm

直角三角形 ABC 内に正方形,長方形,等円 2 個がある。
AB が 27 寸,BC が 36 寸のとき,等円の直径を求めよ。

図のように記号を定め,方程式を解く。

eq6 が https://blog.goo.ne.jp/r-de-r/e/b63974576ba4fdb13955aa322b1ad901 で示した,「鉤股弦と内接円の径」に基づくものである。

using SymPy
@syms R::positive, AB::positive, AE::positive, AD::positive, ae::positive, ad::positive, de::positive, ae::positive;

DE = ae
BD = ae
AB = 27
eq1 = AE - 3//5*AD;
eq2 = ae - 3//5*ad;
eq3 = de - 4//5*ad;
eq4 = AE - DE/de*ae;
eq5 = AD + BD - AB;
eq6 = ae + de - ad - R;

solve([eq1, eq2, eq3, eq4, eq5, eq6], (R, ae, de, ad, AE, AD))

   1-element Vector{NTuple{6, Sym}}:
    (8, 12, 16, 20, 9, 15)

等円の径は 8 である。また,正方形の一辺は 12,長方形の長辺は 20,短辺は 4.8 である。

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 distance(x1, y1, x2, y2, x0, y0)
 p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
 l = sympy.Line(p1, p2)
 l.distance(sympy.Point(x0, y0))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, ae, de, ad, AE, AD) = (8, 12, 16, 20, 9, 15)
   r = R/2
   cosine, sine = (4, 3) ./ 5
   Ax, Ay = (0, 27)
   Bx, By = (0, 0)
   Cx, Cy = (36, 0)
   plot([Ax, Bx, Cx, Ax], [Ay, By, Cy, Ay], lw=0.5, showaxis=false)
   Dx, Dy = (0, 27 - AD)
   Ex, Ey = AE*cosine, AB - AE*sine
   DE = AD*cosine
   ex, ey = (Ex + DE*cosine, Ey - DE*sine)
   ax, ay = (ex - DE*sine, ey - DE*cosine)
   dx, dy = (ex + de*cosine, ey - de*sine)
   plot!([Ex, Dx, ax, ex, Ex], [Ey, Dy, ay, ey, Ey], color=:black, lw=0.5)
   plot!([ax, ax, dx, dx, ax], [ay, 0, 0, ay, ay], color=:black, lw=0.5)
   circle(r, r, r)
   # 右の等円の中心座標 -- 2辺までの距離が等しいとの方程式を解く
   @syms Ox, Oy
   eq1 = distance(ax, ay, ex, ey, Ox, Oy) - r
   eq2 = distance(dx, dy, ex, ey, Ox, Oy) - r
   res = solve([eq1, eq2])[3]
   circle(res[Ox], res[Oy], r)
   if more
       println("正方形の一辺 = $DE")
       println("長方形の辺 = $ad, $dy")
       println("等円の径と中心座標 = $R, $((res[Ox], res[Oy]))")
       point(Ax, Ay, "A ", :black, :right)
       point(Bx, By, "B ", :black, :right)
       point(Cx, Cy, " C", :black)
       point(Dx, Dy, "  D", :black)
       point(Ex, Ey, "  E", :black)
       point(ex, ey, "  e", :black)
       point(ax, ay, " a", :black)
       point(dx, dy, "  d", :black)
       point(r, r, "(r,r)", :red, :center)
       point(res[Ox], res[Oy], "(Ox,Oy)", :red, :center)
       point((Dx + ex)/2, (Dy + ey)/2, "正方形", :center, mark=false)
       point((ax + dx)/2, dy/2, "長方形", :center, mark=false)
       point(r, 1.5r, "等円", :center, mark=false)
       point(res[Ox], 0.5r + res[Oy], "等円", :center, mark=false)
   end
end;

   正方形の一辺 = 12.0
   長方形の辺 = 20, 4.8000000000000025
   等円の径と中心座標 = 8, (17.6000000000000, 8.80000000000000)

 

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

和算の心(その001)

2022年12月22日 | Julia

和算の心(その001)

算額の問題を解くための定石をまとめておく。

直角三角形において,直角を挟む短い辺(鉤:こう)と長い辺(股:こ)と直角に対する辺(弦:げん)をいう(鉤は勾とも書かれる)。また「鉤股弦」は直角三角形のこと。ちなみに,「鉤股弦の定理」は三平方の定理(ピタゴラスの定理)のこと。

この直角三角形に内接する円の径(直径)は,鉤,股,弦を用いて以下のように表される。

径 = 鉤 + 股 - 弦

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 segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end

function kou_ko_gen(鉤, 股, 弦)
   println("鉤 = $鉤;  股 = $股;  弦 = $弦")
   isapprox(鉤^2 + 股^2, 弦^2) || (println("鉤,股,弦の直角三角形ではない"); return)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   径 = 鉤 + 股 - 弦
   println("径 = $径")
   r = 径/2
   plot(showaxis=false)
   ox, oy = 股-r, r
   x, y = (股 - r) .* ([股, 鉤] ./ 弦)
   circle(股-r, r, r, :cyan)
   # 内接円の中心から,鉤,股,弦への垂線
   segment(ox, oy, ox, 0, :blue)
   segment(ox, oy, 股, oy, :blue)
   segment(ox, oy, x, y, :blue)
   # 鉤,股,弦
   segment(0, 0, x, y, :green)
   segment(0, 0, 股-r, 0, :green)
   segment(股, 0, 股-r, 0, :blue)
   segment(股, 0, 股, r, :blue)
   segment(股, 鉤, 股, r, :red)
   segment(股, 鉤, x, y, :red)
   point(0.2, 鉤, "鉤 = 赤 + 青\n股 = 緑 + 青\n弦 = 緑 + 赤\n" *
       "辺々足して\n鉤 + 股 + 弦\n = 2赤 + 2緑 + 2青\n" *
       " = 2(鉤 - 半径) + 2(股 - 半径) + 2半径\n" *
       "移項,整理して,\n鉤 + 股 - 弦 = 直径\n", :black, mark=false)
end;

kou_ko_gen(3, 8, sqrt(9+64))

   鉤 = 3;  股 = 8;  弦 = 8.54400374531753
   径 = 2.4559962546824696

kou_ko_gen(3, 4, 6)

   鉤 = 3;  股 = 4;  弦 = 6
   鉤,股,弦の直角三角形ではない

 

 

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

算額(その69)

2022年12月21日 | Julia

算額(その69)

長野県長野市篠井 長谷観音堂 亨和3年(1803)
http://www.wasan.jp/nagano/hasekannon.html

直角三角形の中に大円と小円が斜線で区切られて入っている。大円,小円の径がそれぞれ 32 尺,24 尺である。 

鉤,股,弦を求めよ。

注:直角三角形の直角を挟む短い辺(鉤)と長い辺(股)と直角に対する辺(弦)

図のように記号を定め,方程式を解く。

using SymPy
@syms AB::positive, BC::positive, BD::positive, AD::positive, CD::positive;

各線分が整数になるように BD=12 と仮置して割合を求める。

eq1 = BD^2 + AD^2 - AB^2;
eq2 = BD^2 + CD^2 - BC^2;
eq3 = 16BC - 12AB;
eq4 = 16BD - 12AD;
eq5 = BD - 12;

solve([eq1, eq2, eq3, eq4, eq5], (AB, BC, BD, AD, CD))

   1-element Vector{NTuple{5, Sym}}:
    (20, 15, 12, 16, 9)

三角形 ABD において,AB, BD, AD の長さは実際は 20n, 12n, 16n (n は整数) であり,その中に半径 16 の円が内接している。

Ag = Ae, Bg = Bf, De = Df = 16

Ag + Bg = AB = 20n

Ae + De = Ae + 16 = AD = 16n より,Ae = 16n - 16 = Ag

Bf + Df = Bf + 16 = BD = 12n より,Bf = 12n - 16 = Bg

Ag + Bg = 28n - 32 = AB = 20n より,8n = 32 で,n = 4

すなわち,AB, BD, AD の長さは実際は 20n=80, 12n=48, 16n=64

三角形 ABD とBCD は相似比 4:3 なので BC = 80*3/4 = 60, CD = 48*3/4 = 36

以上から,直角三角形の直角を挟む短い辺(鉤)の長さは BC = 60, 長い辺(股)の長さは AB = 80,直角に対する辺(弦)の長さは AD + CD = 64 + 36 = 100 であることがわかる。

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 segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot([0, 80, 80, 0], [0, 0, 60, 0], xlims=(-5, 85), ylims=(-5, 65))
   circle(48, 16, 16)
   circle(80 - 12, 36, 12)
   x3, y3 = 48 .* [0.6, 0.8]
   segment(80, 0, 80 - x3, y3, :green)
   if more
       plot!(showaxis=false)
       x1, y1 = 48 .* [0.8, 0.6]
       segment(48, 16, x1, y1)
       segment(48, 16, 48, 0)
       x2, y2 = 32 .* [0.6, 0.8]
       segment(48, 16, 80 - x2, y2)
       x4, y4 = 12 .* [0.8, 0.6]
       segment(80 - 12, 36, 80, 36)
       segment(80 - 12, 36, 80 - 12 - x4, 36 - y4)
       x5, y5 = 12 .* [0.6, 0.8]
       segment(80 - 12, 36, 80 - 12 - x5, 36 + y5)
       point(x1, y1, " e", :blue)
       point(80-x2, y2, "f  ", :blue, :right)
       point(48, 0, "g  ", :blue, :right)
       point(80-x3, y3, "D  ", :blue, :right)
       point(80 - 12 - x4, 36 - y4, "h  ", :blue, :right)
       point(80 - 12 - x5, 36 + y5, " i", :blue)
       point(80, 36, " k", :blue)
       point(0, 0, "A ", :blue, :right)
       point(80, 0, " B", :blue)
       point(80, 60, " C", :blue)
       point(48, 16, " O1\n R1=16", :blue)
       point(80-12, 36, " O2\n R2=12", :blue)
   end
end;

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

算額(その68)

2022年12月21日 | Julia

算額(その68)

新潟県柏崎市 椎谷観音堂 寛政7年(1795)5月
http://www.wasan.jp/niigata/shiiyahukugen.png

外円の中に、弦と外円に接する甲円 1 個,乙円 2 個がある。外円の径が 5 寸であるとき,弦と甲円の径の差が最も大きくなるような乙円の径を求めよ。

答えは 1 寸

外円の半径を 5 とし,図のように記号を定め,方程式を解く。

弦(の 1/2)と r2 の差は以下である。

using SymPy
@syms r2::positive
gen = sqrt((5^2 - (5 - 2r2)^2)) - r2
gen |> println

   -r2 + sqrt(25 - (5 - 2*r2)^2)

これが最も大きくなるのは gen の導関数が 0 になるときの r2 である。

r2 = solve(diff(gen))[1]
println("r2 = $r2 = $(r2.evalf())")

   r2 = 5/2 - sqrt(5)/2 = 1.38196601125011

using Plots
plot(gen, color=:black, aspectratio=:none, ylims=(0, 3.5), size=(500, 300), label="")
xlabel!("y2")
ylabel!("gen")
vline!([r2], label="")



r2 = 5/2 - sqrt(5)/2 として,二本の方程式を立て,解く。

@syms x1, r1
r2 = 5//2 - sqrt(Sym(5))/2
eq1 = x1^2 +(5 - 2r2 + r1)^2 - (5 - r1)^2;
eq2 = x1^2 + (r2 - r1)^2 - (r2 + r1)^2;

solve([eq1, eq2], (x1, r1))

   2-element Vector{Tuple{Sym, Sym}}:
    (-sqrt(10 - 2*sqrt(5)), 1)
    (sqrt(10 - 2*sqrt(5)), 1)

2 組目の解が適切な解である。すなわち,乙円の径は 1 である
また,乙円の中心座標は,(sqrt(10 - 2*sqrt(5)), 1+sqrt(5)) である。

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")
   plot(ylims=(-0.2, 5.1))
   r2 = 5/2 - sqrt(5)/2
   x1, r1 = (sqrt(10 - 2*sqrt(5)), 1)
   gen = sqrt(5^2 - (5 - 2r2)^2)
   println("x1 = $x1;  r1 = $r1;  r2 = $r2")
   println("弦 = $gen;  差 = $(gen - r2)")
   circle(0, 0, 5)
   circle(0, 5-r2, r2, :blue)
   circle(x1, 5-2r2+r1, r1, :magenta)
   hline!([5-2r2])
   if more
       point(0, 5-r2, "5-r2 ", :blue, :right)
       point(0, 5-2r2, "5-2r2 ", :blue, :right)
       point(x1, r1+√5, "(x1,r1+√5))", :blue, :center)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   x1 = 2.3511410091698925;  r1 = 1;  r2 = 1.381966011250105
   弦 = 4.47213595499958;  差 = 3.0901699437494745

 

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

算額(その67)

2022年12月20日 | Julia

算額(その67)

新潟県柏崎市 椎谷観音堂 寛政7年(1795)5月
http://www.wasan.jp/niigata/shiiyahukugen.png

長野県下高井郡 水穂神社 寛政12年正月
http://www.wasan.jp/nagano/mizuho1.html

外円の中に、5 本の斜線と甲円 1 個、乙円 5 個がある。甲円の径が 2889 寸であるとき、乙円の径を求めよ。

答えは 2584 寸

図のように記号を定め,方程式を解く。

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 point0(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 distance(x1, y1, x2, y2, x0, y0)
  p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
  l = sympy.Line(p1, p2)
  l.distance(sympy.Point(x0, y0))
end;

function pentagon(x, y)
   for i in 1:5
       plot!([x[i], x[(i+2)%5+1]], [y[i], y[(i+2)%5+1]], color=:magenta, lw=0.5)
   end
end

using SymPy
甲円,乙円の中心座標と線分 (x[2], y2]), (x[4], y[4]) の距離がそれぞれの円の半径である
function get_r1_r2(x, y) # 2 円の半径を求める
   r1 = distance(x[2], y[2], x[4], y[4], 0, 0)
   @syms r2::poitive
   eq1 = distance(x[2], y[2], x[4], y[4], 0, 1-r2) - r2
   res = solve(eq1)
   return r1.evalf(), res[1].evalf()
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   point = -90:72:270
   x = cosd.(point)
   y = sind.(point)
   pentagon(x, y)
   r1, r2 = get_r1_r2(x, y)
   println("r1 = $r1;  r2 = $r2;  1-r2 = $(1-r2)")
   circle(0, 0, 1)
   circle(0, 0, r1)
   for i = 90:72:378
       circle((1-r2)*cosd(i), (1-r2)*sind(i), r2, :blue)
   end
   if more
       point0(x[2], y[2], "(x[2],y[2])", :blue, :center)
       point0(x[4], y[4], "(x[4],y[4])   ", :blue, :right)
       point0(0, r1, "r1 ", :blue, :right)
       point0(0, 0, "0 ", :blue, :right)
       point0(0, 0, " 甲円", :blue)
       point0(0, 1-r2, "1-r2 ", :blue, :right)
       point0(0, 1-r2, " 乙圓", :blue,)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5, xlims=(-1.15, 1.15))
   end
end;

   r1 = 0.309016994374948;  r2 = 0.276393202250021;  1-r2 = 0.723606797749979

甲円の径が 2889 なら,乙円の径は 2584.0001547987517 となる。

r1 = 0.309016994374948;  r2 = 0.276393202250021
2889/r1 * r2


   2584.0001547987517

 

 

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

Windows の JupyterLab desktop で使用する Julia のバージョンを上げる

2022年12月17日 | ブログラミング

まずは、Julia, JupyterLab desktop のバージョンを上げる

次に、バージョンアップ後の Julia を立ち上げ、

pkg > add IJulia

pkg > build IJulia

そのあと、Juila を抜ける(その状態で作業したければどうぞ)

JupyterLab desktop を立ち上げ、Select Kernel で新しいバージョンを選択できるようになっているはず。

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

算額(その66)

2022年12月17日 | Julia

算額(その66)

東京都中央区茅場町 薬師堂 寛政元年(1789)6月
http://www.wasan.jp/tokyo/yakusi.html

福島県白河市鹿島 鹿島神社
http://www.wasan.jp/fukusima/kashima1.html
に,寛政3年(1791)に奉納されたものの復元算額の中に同じ問題がある。

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2(147)
長野県長野市鬼無里 松厳寺 天保10年(1839)
にもある。

甲円,乙円の径がそれぞれ 69 寸,23 寸である。外円の径を求めよ。

図のように記号を定め,方程式を解く。

注:外円の中心は乙円の円周上にあるのではない。

using SymPy
@syms y::negative, y1::negative, y2::negative, R::positive;
r = 23
eq0 = y2 - y1 - r*sqrt(Sym(3));
eq1 = 4r^2 + y2^2 - 16r^2  # 甲円と乙円2が外接する
eq2 = r^2 + (y - y1)^2 - (R - r)^2  # 外円と乙円1が内接する
eq3 = 4r^2 + y^2 - (R - 3r)^2;  # 外円と甲円が内接する
res = solve([eq0, eq1, eq2, eq3], (y, y1, y2, R))

   1-element Vector{NTuple{4, Sym}}:
    (-14*sqrt(3), -69*sqrt(3), -46*sqrt(3), 121)

R は 121 となった。

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")
   plot()
   r = 23
   (y, y1, y2, R) = (-14*sqrt(3), -69*sqrt(3), -46*sqrt(3), 121)
   println("R = $R;  y = $y, y1 = $y1;  y2 = $y2")
   circle(0, y, R)
   circle(2r, 0, 3r, :green)
   circle(-2r, 0, 3r, :green)
   circle(0, 0, r, :blue)
   circle(0, y2, r, :blue)
   circle(r, y1, r, :blue)
   circle(-r, y1, r, :blue)
   if more
       point(0, 0, "0 ", :blue, :right)
       point(r, 0, "r ", :blue, :right)
       point(0, y, "外円の中心 y", :red, :center)
       point(2r, 0, "2r ", :green, :right)
       point(5r, 0, "5r ", :green, :right)
       point(r, y1, "(r,y1)", :blue, :center)
       point(0, y2, "y2 ", :blue, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(25, -80, "小円2", :blue, mark=false)
       point(43, -100, "小円1", :blue, mark=false)
   end
end;

   R = 121;  y = -24.24871130596428, y1 = -119.51150572225252;  y2 = -79.67433714816835

 

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

算額(その65)

2022年12月16日 | Julia

算額(その65)

長野県木曽郡上松町寝覚 臨川寺弁財天社 文政13年
http://www.wasan.jp/nagano/rinsenji.html

3 個の甲円が交差している部分に内接する乙円の径を求めよ。

甲円の半径を 1 として,図のように記号を定め,方程式を解く。
x, y は簡単に決められる。方程式は,甲円と乙円が内接することに基づく。

using SymPy
@syms x::positive, y::positive, r::positive;
x = sqrt(Sym(3))/4
y = 1//4
eq = x^2 + (1 - y)^2 - (1 - r)^2;

res = solve(eq) |> println

   Sym[1 - sqrt(3)/2, sqrt(3)/2 + 1]

解は 2 つ得られるが,1 - sqrt(3)/2 = 0.1339745962155614 が適切な解である。

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")
   plot()
   r = 1 - sqrt(3)/2
   println("r = $r")
   circle(0, 1, 1, :green)
   circle(cos(pi/6), -sin(pi/6), 1, :green)
   circle(-cos(pi/6), -sin(pi/6), 1, :green)
   circle(x, y, r)
   circle(-x, y, r)
   circle(0, -1/2, r)
   if more
       point(x, y, "(x,y)", :red, :center)
       point(0, 1, "1 ", :green, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   r = 0.1339745962155614

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

算額(その64)

2022年12月16日 | Julia

算額(その64)

長野県飯山市滝澤家所蔵
http://www.wasan.jp/nagano/takizawa.html

大円の中に正方形,中円,小円が収まっている。大円の径を 5 寸としたとき,中円の中心を結ぶ正方形の一辺の長さと小円の径を求めよ。

大円の半径を5として,図のように記号を定め,方程式を解く。

using SymPy
@syms r1::positive, r2::positive;
eq1 = 2r2 - 5/sqrt(Sym(2));
eq2 = 2r2^2 - (r1 + r2)^2;
res = solve([eq1, eq2], (r1, r2));
res[1][1] |> println  # r1
res[1][2] |> println  # r2

   5/2 - 5*sqrt(2)/4
   5*sqrt(2)/4

正方形の一辺は 2r2。算額では 3寸5分3厘としている。

2 * 5*sqrt(2)/4

   3.5355339059327378

小円の径は 5/2 - 5*sqrt(2)/4。算額では 7分8厘としている。

5/2 - 5*sqrt(2)/4 |> N

   0.7322330470336311

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")
   plot()
   (r1, r2) = (5/2 - 5*sqrt(2)/4, 5*sqrt(2)/4)
   println("r1 = $r1, r2 = $r2")
   circle(0, 0, 5, :black)
   circle(0, 0, r1)
   circle(r2, r2, r2, :brown)
   circle(r2, -r2, r2, :brown)
   circle(-r2, r2, r2, :brown)
   circle(-r2, -r2, r2, :brown)
   x = 5/sqrt(2)
   plot!([-x, x, x, -x, -x], [-x, -x, x, x, -x], color=:black, lw=0.5)
   plot!([-r2, r2, r2, -r2, -r2], [-r2, -r2, r2, r2, -r2], color=:blue, lw=0.5, linestyle=:dot)
   if more
       point(0, 0, "0 ", :black, :right)
       point(r1, 0, "r1 ", :black, :right)
       point(r2, 0, "r2 ", :brown, :right)
       point(5/√2, 0, "5/√2 ", :brown, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.7322330470336311, r2 = 1.7677669529663689

 

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

算額(その63)

2022年12月16日 | Julia

算額(その63)

飯山市常磐柳新田 瀧澤治夫宅所蔵 天保11年(1840)7月
http://www.wasan.jp/nagano/takizawa2.html

半円の中に,3 種類の円が収まっている。各円の径を求めよ。

半円の半径を1として,図のように記号を定め,方程式を解く。

using SymPy
@syms x1::positive, r1::positive, x2::positive, r2::positive;
eq1 = x1^2 + r1^2 - (1 - r1)^2;
eq2 = x2^2 + r2^2 - (1 - r2)^2;
eq3 = x1^2 + (1//2 - r1)^2 - (1//2 + r1)^2;
eq4 = (x2 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2;

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

   1-element Vector{NTuple{4, Sym}}:
    (sqrt(2)/2, 1/4, 2*sqrt(2)/3, 1/18)

算額では,甲円は1尺2寸(12寸),乙円は8寸,丙円は結果として 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 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")
   plot()
   (x1, r1, x2, r2) = (sqrt(2)/2, 1/4, 2*sqrt(2)/3, 1/18)
   println("x1 = $x1, r1 = $r1\nx2 = $x2, r2 = $r2")
   circle(0, 0, 1, :black)
   circle(0, 1/2, 1/2)
   circle(x1, r1, r1, :brown)
   circle(-x1, r1, r1, :brown)
   circle(x2, r2, r2, :green)
   circle(-x2, r2, r2, :green)
   hline!([0], color=:black, lw=0.5)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, 1/2, "1/2 ", :black, :right)
       point(x1, r1, "(x1,r1)", :brown, :center)
       point(x2, r2, "\n    (x2,r2)", :green, :center)
       vline!([0], color=:black, lw=0.5)
   end
end;

   x1 = 0.7071067811865476, r1 = 0.25
   x2 = 0.9428090415820635, r2 = 0.05555555555555555

 

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

算額(その62)

2022年12月16日 | Julia

算額(その62)

長野市篠ノ井 布制神社 文化3年(1806)
http://www.wasan.jp/nagano/fuse.html

玄が12寸,甲円の径が2寸5分のとき丙円の径を求めよ。

玄を 240,甲円の半径を 25 とし,図のように記号を定め,方程式を解く。

using SymPy
@syms x1::positive, r1::positive, x2::positive, r2::positive, y0::negative, r0::positive;
eq1 = x1^2 + (25 - r1)^2 - (25 + r1)^2;
eq2 = x2^2 + (25 - r2)^2 - (25 + r2)^2;
eq3 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2;
eq4 = y0^2 + 120^2 - r0^2;
eq5 = x1^2 + (r1 - y0)^2 - (r0 - r1)^2;
eq6 = 50 - y0 - r0;

res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (x1, r1, x2, r2, y0, r0))

   2-element Vector{NTuple{6, Sym}}:
    (600/13, 3600/169, 24, 144/25, -119, 169)
    (600/13, 3600/169, 600, 3600, -119, 169)

丙円の半径は算額では 5分7厘5毛となっているが,正しくは 5分7厘6毛(5.76)となる。

また,算額の各円の見かけの大きさは,計算結果とはかなり異なっている。

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")
   plot()
   (x1, r1, x2, r2, y0, r0) = (600/13, 3600/169, 24, 144/25, -119, 169)
   println("x1 = $x1, r1 = $r1\nx2 = $x2, r2 = $r2\ny0 = $y0, r0 = $r0")
   circle(0, y0, r0, :black)
   circle(0, 25, 25)
   circle(x1, r1, r1, :brown)
   circle(-x1, r1, r1, :brown)
   circle(x2, r2, r2, :green)
   plot!([-120, 120], [0, 0], color=:black, lw=0.5)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, 25, "25 ", :black, :right)
       point(120, 0, " 120", :black)
       point(x1, r1, "(x1,r1)", :brown, :center)
       point(x2, r2, "\n    (x2,r2)", :green, :center)
       point(0, y0, "y0 ", :blue, :right)
       vline!([0], color=:black, lw=0.5)
   end
end;

   x1 = 46.15384615384615, r1 = 21.301775147928993
   x2 = 24, r2 = 5.76
   y0 = -119, r0 = 169

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

算額(その61)

2022年12月16日 | Julia

算額(その61)

長野市鬼無里 智光山文珠堂 明治11年(1878)8月
http://www.wasan.jp/nagano/monjudo1.html

半円の中に 5 個の円がある。それぞれの径を求めよ。

半円の半径を 1 とし,図のように記号を定め,方程式を解く。


using SymPy
@syms x1::positive, r1::positive, x2::negative, r2::positive,
   x3::positive, r3::positive, x4::negative, r4::positive,
   x5::positive, r5::positive;
eq1 = x1^2 + r1^2 - (1 - r1)^2
eq2 = x3^2 + r3^2 - (1 - r3)^2
eq3 = x2^2 + r2^2 - (1 - r2)^2
eq4 = (x1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq5 = (x1 - x5)^2 + (r1 - r5)^2 - (r1 + r5)^2
eq6 = (x5 - x3)^2 + (r5 - r3)^2 - (r5 + r3)^2
eq7 = (x2 - x1)^2 + (r2 - r1)^2 - (r2 + r1)^2
eq8 = (x4 - x1)^2 + (r4 - r1)^2 - (r4 + r1)^2
eq9 = (x2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq10 = 1/r1 + 1/r3 + 2sqrt(1/r1/r3) - 1/r5;  # デカルトの円定理

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10])

この問題の場合も solve() では解けないので nlsolve() を用いる。

eq1 |> expand |> println
eq2 |> expand |> println
eq3 |> expand |> println
eq4 |> expand |> println
eq5 |> expand |> println
eq6 |> expand |> println
eq7 |> expand |> println
eq8 |> expand |> println
eq9 |> expand |> println
eq10 |> expand |> println


   2*r1 + x1^2 - 1
   2*r3 + x3^2 - 1
   2*r2 + x2^2 - 1
   -4*r1*r3 + x1^2 - 2*x1*x3 + x3^2
   -4*r1*r5 + x1^2 - 2*x1*x5 + x5^2
   -4*r3*r5 + x3^2 - 2*x3*x5 + x5^2
   -4*r1*r2 + x1^2 - 2*x1*x2 + x2^2
   -4*r1*r4 + x1^2 - 2*x1*x4 + x4^2
   -4*r2*r4 + x2^2 - 2*x2*x4 + x4^2
   -1/r5 + 1/r3 + 1/r1 + 2/(sqrt(r1)*sqrt(r3))

using NLsolve

function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params...), [ini], ftol=1e-14)
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
v = r.zero
end
return v, r.f_converged
end;

function H(u)
x1, r1, x2, r2, x3, r3, x4, r4, x5, r5 = u
return [
       2*r1 + x1^2 - 1
       2*r3 + x3^2 - 1
       2*r2 + x2^2 - 1
       -4*r1*r3 + x1^2 - 2*x1*x3 + x3^2
       -4*r1*r5 + x1^2 - 2*x1*x5 + x5^2
       -4*r3*r5 + x3^2 - 2*x3*x5 + x5^2
       -4*r1*r2 + x1^2 - 2*x1*x2 + x2^2
       -4*r1*r4 + x1^2 - 2*x1*x4 + x4^2
       -4*r2*r4 + x2^2 - 2*x2*x4 + x4^2
       -1/r5 + 1/r3 + 1/r1 + 2/(sqrt(r1)*sqrt(r3))
       ];
end;

iniv = [0.01, 0.5, -0.7, 0.22, 0.8, 0.19, -0.5, 0.05, 0.4, 0.05];
res = nls(H, ini=iniv)

   ([0.21236918980148334, 0.47744966361153107, -0.5821590777643445, 0.3305454040882841, 0.7994277491932246, 0.18045763690992736, -0.22131238148404836, 0.0984814314508113, 0.5759211572513911, 0.06920626565999469], true)

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")
   plot(ylims=(-0.07, 1.07))
   x1, r1, x2, r2, x3, r3, x4, r4, x5, r5 = [0.21236918980148334, 0.47744966361153107, -0.5821590777643445, 0.3305454040882841, 0.7994277491932246, 0.18045763690992736, -0.22131238148404836, 0.0984814314508113, 0.5759211572513911, 0.06920626565999469]
   println("x1 = $x1, r1 = $r1\nx2 = $x2, r2 = $r2\nx3 = $x3, r3 = $r3\nx4 = $x4, r4 = $r4\nx5 = $x5, r5 = $r5")
   circle(0, 0, 1, :black)
   circle(x1, r1, r1, :brown)
   circle(x2, r2, r2, :green)
   circle(x3, r3, r3, :blue)
   circle(x4, r4, r4, :red)
   circle(x5, r5, r5, :magenta)
   hline!([0], color=:black, lw=0.5)
   if more
       point(0, 0, " 0", :black)
       point(x1, r1, "(x1,r1)", :brown, :center)
       point(x2, r2, "(x2,r2)", :green, :center)
       point(x3, r3, "(r3,r3)", :blue, :center)
       point(x4, r4, "(x4,r4)", :red, :center)
       point(x5, r5, "(x5,r5)", :magenta, :center)
       vline!([0], color=:black, lw=0.5)
   end
end;

   x1 = 0.21236918980148334, r1 = 0.47744966361153107
   x2 = -0.5821590777643445, r2 = 0.3305454040882841
   x3 = 0.7994277491932246, r3 = 0.18045763690992736
   x4 = -0.22131238148404836, r4 = 0.0984814314508113
   x5 = 0.5759211572513911, r5 = 0.06920626565999469

 

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

算額(その60)

2022年12月15日 | Julia

算額(その60)

石川県松任市若宮 若宮八幡宮 文政13年吉祥日
http://www.wasan.jp/isikawa/wakamiyahatiman1.html

石川県鳳珠郡穴水町 美麻奈比古神社 文政13年8月
http://www.wasan.jp/isikawa/mimanahiko.html

群馬県高崎市 祖師堂 明治16年
http://www.wasan.jp/gunma/sosido.html

外円の中に 8 個の小円と 1 個の大円がある。それぞれの径を求めよ。

外円の半径を 1 とし,図のように記号を定め,方程式を解く。

using SymPy
@syms r1::positive, x2::positive, y2::positive, r2::positive;

x2 = (1 - r2) / sqrt(Sym(2))
y2 = x2
eq1 = r1 + 2r2 - 1 |> expand
eq2 = x2^2 + y2^2 - (1 - r2)^2 |> expand
eq3 = x2^2 + (1 - r2 - y2)^2 - 4r2^2 |> expand;

res = solve([eq1, eq3], (r1, r2))

   2-element Vector{Tuple{Sym, Sym}}:
    (-4*sqrt(2) - 2*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 7, -3 + sqrt(2)*sqrt(10 - 7*sqrt(2)) + 2*sqrt(2))
    (-4*sqrt(2) + 2*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 7, -3 - sqrt(2)*sqrt(10 - 7*sqrt(2)) + 2*sqrt(2))

解が 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")
   plot()
   (r1, r2) = (-4*sqrt(2) - 2*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 7, -3 + sqrt(2)*sqrt(10 - 7*sqrt(2)) + 2*sqrt(2))
   x2 = y2 = (1 - r2)/sqrt(2)
   println("r1 = $r1, x2 = $x2, y2 = $y2, r2 = $r2")
   circle(0, 0, 1, :black)
   circle(0, 0, r1, :brown)
   circle(x2, y2, r2)
   circle(x2, -y2, r2)
   circle(-x2, y2, r2)
   circle(-x2, -y2, r2)
   circle(0, 1-r2, r2, :green)
   circle(0, r2-1, r2, :green)
   circle(1-r2, 0, r2, :magenta)
   circle(r2-1, 0, r2, :magenta)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, 1-r2, "1-r2 ", :green, :right)
       point(0, r1, "r1 ", :brown, :right)
       point(1-r2, 0, "1-r2 ", :green, :right)
       point(r1, 0, "r1 ", :brown, :right)
       point(x2, y2, "(x2,y2)", :red, :center)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.4464626921716892, x2 = 0.5114017891839755, y2 = 0.5114017891839755, r2 = 0.2767686539141554

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

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

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