裏 RjpWiki

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

算額(その724)

2024年02月24日 | Julia

算額(その724)

群馬県桐生市梅田町 日枝神社 明治12年
https://gunmawasan.web.fc2.com/files/sangak-corner/hiejin.htm

長方形と輪(円)が交わってできる領域に 3 個の小円がある。
また,下部に 2 ほんの斜線を入れ,正方形を 3 個入れる。左右の正方形の 2 つの頂点は輪の円周上と斜線の上にある。中央の正方形の上辺は輪と接している。
小円の著系が 1 寸のとき,正方形の一辺の長さはいかほどか。

長方形の短辺と長辺を 2a, b とする。
輪の半径と中心座標を r1, (0, 2c + r1); a = r1
小円の半径と中心座標を r2, (a - r2, b - r2)
正方形の一辺の長さを 2c
として,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive, r1::positive, r2::positive, y0::positive
a = r1
r2 = 1//2
eq1 = 2c + 2r1 - 2r2 - b
eq2 = y0/2c - 2c/(a - c)
eq3 = (a - 2c)^2 + (y0 + 2c - 2c - r1)^2 - r1^2
eq4 = ( a - r2)^2 + (b - r2 - 2c - r1)^2 - (r1 + r2)^2;
res = solve([eq1, eq2, eq3, eq4], (b, c, r1, y0))

   4-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (2*r2/5, r2/5, r2, r2/5)
    (r2, r2/2, r2, 2*r2)
    (98*r2/5, 9*r2/5, 9*r2, 9*r2/5)
    (25*r2, 9*r2/2, 9*r2, 18*r2)

4 組の解が得られるが,3 番目のものが適解である。

r2 = 1/2 のとき,
(b, c, r1, y0) =(98*r2/5, 9*r2/5, 9*r2, 9*r2/5) = (9.8, 0.9, 4.5, 0.9) である。
すなわち,正方形の一辺の長さは小円の直径の 1.8 倍,すなわち 1 寸 8 分である。

r2 = 1/2
(98*r2/5, 9*r2/5, 9*r2, 9*r2/5)

   (9.8, 0.9, 4.5, 0.9)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   (b, c, r1, y0) = (98*r2/5, 9*r2/5, 9*r2, 9*r2/5)
   @printf("正方形の一辺の長さ = %g\n", 2c)
   @printf("b = %g;  c = %g;  r1 = %g;  y0 = %g\n", b, c, r1, y0)
   a = r1
   plot([a, a, -a, -a, a], [0, b, b, 0, 0], color=:black, lw=0.5)
   circle(0, 2c + r1, r1)
   circle(0, b + r2, r2, :blue)
   circle(a - r2, b - r2, r2, :blue)
   circle(r2 - a, b - r2, r2, :blue)
   rect(-c, 0, c, 2c, :green)
   rect(a - 2c, y0, a, y0 + 2c, :green) 
   rect(2c - a, y0, -a, y0 + 2c, :green)
   segment(a, 0, c, 2c)
   segment(-a, 0, -c, 2c)
   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(a - 2c, y0 + 2c, " (a-2c,y0+2c)", delta=-delta/2)
       point(a - 2c, y0, " (a-2c,y0)", :green, :left, :bottom, delta=delta/2)
       point(c, 0, " c", :green, :left, :bottom, delta=delta/2)
       point(0, 2c, " 2c", :green, :left, :bottom, delta=delta/2)
       point(0, 2c + r1, " 輪:r1,(0,2c+r1)", :red, :left, :vcenter)
       point(a - r2, b - r2, "小円:r2,(a-r2,b-r2)    ", :blue, :right, :vcenter)
       point(0, b + r2, " b+r2", :black, :center, :bottom, delta=delta/2)
       point(a, b, "(a,b)", :black, :center, :bottom, delta=delta/2)
   end
end;

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

算額(その723)

2024年02月24日 | Julia

算額(その723)

群馬県桐生市天神町 桐生天満宮 明治11年
https://gunmawasan.web.fc2.com/files/sangak-corner/tenman-g1.htm

外円内に甲円,乙円,丙円,丁円,戊円が互いに接して入っている。乙円,丙円,戊円の直径がそれぞれ 15寸,6 寸,4 寸のとき,外円,甲円,丁円の直径はいかほどか。

算額(その127)を解き直した。
https://blog.goo.ne.jp/r-de-r/e/401c8ab1043169fd4660aec699096fc4

パラメータを減らすために,甲円の中心が y 軸上にあるように回転して考える。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (x4, y4)
戊円の半径と中心座標を r5, (x5, y5)
とおき,以下の連立方程式を解く。
なお,丙円と丁円の外接条件として eq6 を用いようとしたが,感度不足で思うような収束解が得られなかった。
算法助術に活路を見出そうとして,公式77から eq12(または eq13)を得て,やっと解が得られた。

今回の算額もまた,図と実際の解は全く異なるものである。最も,問と答を見れば,図における円の比率がおかしいのはすぐに分かるのではあるが。

include("julia-source.txt");

using SymPy
@syms R, r1, r2, x2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5, t
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq3 = x5^2 + (R - r1 - y5)^2 - (r1 + r5)^2
eq4 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq5 = (x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2
# eq6 = (x3 - x4)^2 - (y3 - y4)^2 - (r3 + r4)^2  # 感度不足
eq7 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq8 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq9 = x2^2 + y2^2 - (R - r2)^2
eq10 = x3^2 + y3^2 - (R - r3)^2
eq11 = x4^2 + y4^2 - (R - r4)^2;
# 算法助術 公式77 より  eq12 か eq13 いずれか
eq12 = r4*r5*t^2 + r1*r5*t^2 + 2r1*r4*t^2 - 2r1*r4*r5*(r2 + r3) -4r2*r1*r3*r4;
# eq13 = -R*r4*t^2 - R*r1*t^2 + 2r1*r4*t^2 +2R*r1*r4*(r2 + r3) - 4r2*r1*r3*r4;

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (R, r1, x2, y2, x3, y3, r4, x4, y4, x5, y5) = u
   return [
       x2^2 - (r1 + r2)^2 + (R - r1 - y2)^2,  # eq1
       x3^2 - (r1 + r3)^2 + (R - r1 - y3)^2,  # eq2
       x5^2 - (r1 + r5)^2 + (R - r1 - y5)^2,  # eq3
       -(r2 + r4)^2 + (x2 - x4)^2 + (y2 - y4)^2,  # eq4
       -(r2 + r5)^2 + (x2 - x5)^2 + (y2 - y5)^2,  # eq5
       -(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq7
       -(r4 + r5)^2 + (x4 - x5)^2 + (y4 - y5)^2,  # eq8
       x2^2 + y2^2 - (R - r2)^2,  # eq9
       x3^2 + y3^2 - (R - r3)^2,  # eq10
       x4^2 + y4^2 - (R - r4)^2,  # eq11
       -4*r1*r2*r3*r4 - 2*r1*r4*r5*(r2 + r3) + 2*r1*r4*t^2 + r1*r5*t^2 + r4*r5*t^2,  # eq12
   ]
end;

(r2, r3, r5) = (15, 6, 4) .// 2
t = sqrt(Sym(2)r2*r3)
iniv = BigFloat[31, 15.5, 21, 7.4, 17, 21, 2.7, 21.5, 17.5, 17, 16]
res = nls(H, ini=iniv)

   ([30.0, 15.0, 21.213203435596427, 7.5, 16.97056274847714, 21.0, 2.5, 21.213203435596427, 17.5, 16.97056274847714, 16.0], true)

外円の直径 = 60;  甲円の直径 = 30;  丁円の直径 = 5 である。
なお,乙円はx 軸に接する。すなわち,y2 = r2 であった。

その他のパラメータは以下のとおりである。

r2 = 7.5;  r3 = 3;  r5 = 2
外円の直径 = 60;  甲円の直径 = 30;  丁円の直径 = 5
R = 30;  r1 = 15;  x2 = 21.2132;  y2 = 7.5;  x3 = 16.9706;  y3 = 21;  r4 = 2.5;  x4 = 21.2132;  y4 = 17.5;  x5 = 16.9706;  y5 = 16

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, r5) = (15, 6, 4) .// 2
   (R, r1, x2, y2, x3, y3, r4, x4, y4, x5, y5) = res[1]
   @printf("r2 = %g;  r3 = %g;  r5 = %g\n", r2, r3, r5)
   @printf("外円の直径 = %g;  甲円の直径 = %g;  丁円の直径 = %g\n", 2R, 2r1, 2r4)
   @printf("R = %g;  r1 = %g;  x2 = %g;  y2 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x4 = %g;  y4 = %g;  x5 = %g;  y5 = %g\n", R, r1, x2, y2, x3, y3, r4, x4, y4, x5, y5)
   plot()
   circle(0, 0, R, :orange)
   circle(0, R - r1, r1)
   circle(x2, y2, r2, :blue)
   circle(x3, y3, r3, :magenta)
   circle(x4, y4, r4, :green)
   circle(x5, y5, r5, :brown)
   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, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, delta=-delta/2)
       point(x3, y3, "丙円:r3,(x3,y3) ", :black, :right, :vcenter)
       point(x4, y4, "丁円:r4,(x4,y4) ", :black, :right, :bottom, delta=delta/2)
       point(x5, y5, "戊円:r5,(x5,y5) ", :black, :right, :vcenter)
   end
end;

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

算額(その722)

2024年02月23日 | Julia

算額(その722)

群馬県桐生市天神町 桐生天満宮 明治11年
https://gunmawasan.web.fc2.com/files/sangak-corner/tenman-g1.htm

正三角形内に 3 個の楕円が入っている。三角形の一辺の長さと楕円の短径が与えられたとき,楕円の長径を求めよ。

正三角形の重心を原点として,
正三角形の一辺の長さを 2c
楕円の長半径と短半径を 2a, 2b
下側の楕円と左上の楕円の接点座標を (x0, y0)
として,下側の楕円についての以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a, b, c, x0, y0

s3 = sqrt(Sym(3))
eq1 = x0^2/a^2 + (y0 - b + c/s3)^2/b^2 - 1
eq2 = -b^2*x0/(a^2*(y0 - b + c/s3)) - 1/s3
eq3 = y0/x0 - 1/s3
res = solve([eq1, eq2, eq3], (a, x0, y0))

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (sqrt(c*(-6*sqrt(3)*b^2 + 9*b*c - sqrt(3)*c^2)/(3*b - sqrt(3)*c)), c*(-6*b + sqrt(3)*c)/(3*b - sqrt(3)*c), c*(-2*sqrt(3)*b + c)/(3*b - sqrt(3)*c))
    (-sqrt(-c*(6*sqrt(3)*b^2 - 9*b*c + sqrt(3)*c^2)/(3*b - sqrt(3)*c)), c*(-6*b + sqrt(3)*c)/(3*b - sqrt(3)*c), c*(-2*sqrt(3)*b + c)/(3*b - sqrt(3)*c))

2 組の解が得られるが,最初のものが適解である。長半径 a の式は以下のようになる。

a = res[1][1]
a |> println

   sqrt(c*(-6*sqrt(3)*b^2 + 9*b*c - sqrt(3)*c^2)/(3*b - sqrt(3)*c))

この式を簡約化するためにルートの中身について考える。

a2 = a^2
a2 |> println

   c*(-6*sqrt(3)*b^2 + 9*b*c - sqrt(3)*c^2)/(3*b - sqrt(3)*c)

分数式であるが,割り切れる。結果の第1要素が商,第2要素が剰余であり,確かに剰余は0になっており,割り切れている。

temp = sympy.div(numerator(a2), denominator(a2))

   (-2*sqrt(3)*b*c + c^2, 0)

楕円の長半径は商の平方根を取れば求まる。

楕円の長半径 = sqrt(temp[1])
楕円の長半径 |> println

   sqrt(-2*sqrt(3)*b*c + c^2)

楕円の長半径は sqrt(-2*sqrt(3)*b*c + c^2) であるが,もう少しわかりやすく書くと,
sqrt(c*(c - √12*b))
となる。この式は「術」に述べられている「12の平方根に短径を掛けて,正三角形の一辺の長さから引き,短径を掛けて,開平する」というのと同じである。
なお,この解法の最初で「正三角形の一辺の長さを 2c,楕円の長半径と短半径を 2a, 2b」としているので,注意が必要である。

正三角形の一辺の長さが 1,楕円の短径が 0.2 のとき,c = 1/2, b = 0.2/1 を与えて長半径を得て,2倍して長径を得る。

楕円の長半径(c => 1/2, b => 0.2/2).evalf()*2 |> println

   0.554237924516582

c = 1, b = 0.2 を代入すると,(2倍しなくても)長径が得られる。

楕円の長半径(c => 1, b => 0.2).evalf() |> println

   0.554237924516582

パラメータをまとめると以下のようになる。

楕円の長径 = 0.554238;  長半径 = 0.277119
c = 0.5
b = 0.1
a = 0.2771189622582913
x0 = -0.2349942264126275
y0 = -0.13567397987733842

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   s3 = sqrt(3)
   (c, b) = (1/2, 2/20)
   (a, x0, y0) = (sqrt(-2*sqrt(3)*b*c + c^2), c*(-6*b + sqrt(3)*c)/(3*b - sqrt(3)*c), c*(-2*sqrt(3)*b + c)/(3*b - sqrt(3)*c))
   @printf("楕円の長径 = %g;  長半径 = %g\n", 2a, a)
   @show(c, b, a, x0, y0)
   plot([c, 0, -c, c], [-c/s3, 2c/s3, -c/s3, -c/s3], color=:blue, lw=0.5)
   ellipse(0, b - c/s3, a, b, color=:red)
   d = c/s3 - b
   ellipse(d*cosd(30), d*sind(30), a, b, color=:red, φ=120)
   ellipse(-d*cosd(30), d*sind(30), a, b, color=:red, φ=240)
   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, 2c/√3, " 2c/√3", :blue, :left, :vcenter)
       point(c, -c/√3, " c", :blue, :left, :bottom, delta=delta/2)
       point(x0, y0, "(x0,y0)", :red, :left, delta=-delta/2)
       point(0, b - c/√3, " 楕円:a,b,(0,b-c/√3)", :red, :center, delta=-delta/2)
   end
end;

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

算額(その721)

2024年02月23日 | Julia

算額(その721)

長野県諏訪市中洲 諏訪大明神社 天保12年(1841)
中村信弥「改訂増補 長野県の算額

http://www.wasan.jp/zoho/zoho.html

正方形,菱形,正三角形の一辺の長さを a, b, c とおく。
a^3 + b^3 = 281827, b^3 + c^3 = 23779, a^(1/2) + b^(1/3) + c^(1/4) = 13 である。
a, b, c 特に,整数解を求めよ。

a, b は最大 281827^(1/3) < 66
b, c は最大 23779^(1/3) < 29
c は最大 28561

以下でブルートフォースで解を求める。

using Printf
for a = 1:65, b = 1:28, c = 1:28
   if isapprox(a^3 + b^3, 281827) && isapprox(b^3 + c^3, 23779) && isapprox(a^(1/2) + b^(1/3) + c^(1/4), 13)
       @printf("a = %g;  b = %g;  c = %g\n", a, b, c)
       break
   end
end   

   a = 64;  b = 27;  c = 16

a = 64,b = 27,c = 16 は条件を満たす整数解である。

実数解を求めてみる。

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (a, b, c) = u
   return [
       a^3 + b^3 - big"281827",  # eq1
       b^3 + c^3 - big"23779",  # eq2
       sqrt(a) + b^(1/3) + c^(1/4) -13.0,  # eq3
   ]
end;

iniv = BigFloat[64.4, 3.03, 230]
res = nls(H, ini=iniv)

   ([64.97494444657413, 19.591250379477188, 25.33394250747033], true)

(a, b, c) = [64.97494444657413, 19.591250379477188, 25.33394250747033]
は条件を満たすので,解である。

function check(a, b, c)
   @show(isapprox(a^3 + b^3, 281827), isapprox(b^3 + c^3, 23779), isapprox(a^(1/2) + b^(1/3) + c^(1/4), 13))
end;

check(64.97494444657413, 19.591250379477188, 25.33394250747033);

   isapprox(a ^ 3 + b ^ 3, 281827) = true
   isapprox(b ^ 3 + c ^ 3, 23779) = true
   isapprox(a ^ (1 / 2) + b ^ (1 / 3) + c ^ (1 / 4), 13) = true

その他の実数解があるか,初期値をランダムに変えて実行してみるが,整数解以外には 1 組しかないようだ。

results = Set()
for i in 1:100
   iniv = BigFloat[(rand(1)*66)[1], (rand(1)*29)[1], (rand(1)*29)[1]]
   res = nls(H, ini=iniv)
   if res[2]
       push!(results, res[1])
   end
end
println(results)

   Set(Any[[64.0, 26.999999999999996, 16.000000000000007], [64.97494444657413, 19.591250379477188, 25.33394250747033]])

 

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

算額(その720)

2024年02月23日 | Julia

算額(その720)

長野県諏訪市中洲 諏訪大明神社 天保12年(1841)
中村信弥「改訂増補 長野県の算額」

http://www.wasan.jp/zoho/zoho.html

正五角形に 2 本の対角線を引き,区画された領域に大円,中円,小円,正方形を入れる。対角線の長さの 8 乗が 256 となるときの,それぞれの大きさはいかほどか。

図形は相似なので,R = 1 として問題を解き,必要な場合に「対角線の長さの 8 乗が 256 となるとき」という条件下の値を提示する。

まず,基本的な座標を定義する。

正五角形が内接する円の直径を R として 5 個の頂点の座標を以下のように定義する。

include("julia-source.txt");

using SymPy
@syms R::positive, r1, ox1, r2, ox2, r3, ox3, xa, a
R = 1
(x1, y1) = R.*(cos(2PI/20), sin(2PI/20))
(x2, y2) = (0, R)
(x3, y3) = (-x1, y1)
(x4, y4) = R.*(-sin(2PI/10), -cos(2PI/10))
(x5, y5) = (-x4, y4);

それぞれの図形のパラメータは独立して求めることができるので,まず大円の直径 r1 と中心座標 (ox1, y1 - r1) を求める。

eq1 = (y1 - y4)/2 - r1
solve(eq1, r1)[1] |> println
eq2 = (x3 + x5)/2 - ox1;
res1 = solve(eq2, ox1)[1] |> println

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

次に,中円の直径 r2 と中心座標 (ox2, y1 + r2) を求める。

@syms d
eq3 = dist(x2, y2, x3, y3, ox2, y1 + r2) - r2^2
eq3 = eq3 = apart(eq3, d) |> simplify
eq4 = dist(x2, y2, x5, y5, ox2, y1 + r2) - r2^2
eq4 = apart(eq4, d) |> simplify
res2 = solve([eq3, eq4], (r2, ox2));
res2[3][1].evalf(), res2[3][2].evalf()

   (0.263932022500210, -0.138757275712888)

次に,小円の直径 r3 と中心座標 (ox3, y1 + r3) を求める。

eq5 = dist(x2, y2, x5, y5, ox3, y1 + r3) - r3^2
eq5 = apart(eq5, d) |> simplify
eq6 = dist(x2, y2, x1, y1, ox3, y1 + r3) - r3^2
eq6 = apart(eq6, d) |> simplify;
res3 = solve([eq5, eq6], (r3, ox3));
res3[1][1].evalf(), res3[1][2].evalf()

   (0.190983005625053, 0.363271264002680)

正方形の一辺の長さを a,右上の頂点の座標を (xa, y1) とおき,以下の連立方程式から a, xa を求める。

eq7 = (x1 - xa)/a - (x1 - x5)/(y1 - y5)
eq8 = (x5 - xa + a)/(y1 - a - y5) - x5/(y2 - y5)
res4 = solve([eq7, eq8], (xa, a))
res4[a].evalf(), res4[xa].evalf()

   (0.440371669705066, 0.807971087145006)

以上により計算された数値は R = 1 のときのものである。
「問」の条件は「対角線の長さの 8 乗が 256 となるとき」,つまり「対角線の長さが 2 のとき」というものである。
R = 1 のとき,対角線の長さは,2x1 = 1.902113032590307 なので 問の条件に従うためには R = 1 で計算された数値を 2/1.902113032590307 = 1.0514622242382672 倍すればよいということになる。

大円の直径 = 1.1755705045849463
中円の直径 = 0.5550291028515503
小円の直径 = 0.40162283177245545
正方形の一辺の長さ = 0.4630341753196083

なお,「答」で提示された数値には誤差があり,また,この数値を 10 倍したものが提示されている。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 1
   (x1, y1) = R.*(cos(2pi/20), sin(2pi/20))
   (x2, y2) = (0, R)
   (x3, y3) = (-x1, y1)
   (x4, y4) = R.*(-sin(2pi/10), -cos(2pi/10))
   (x5, y5) = (-x4, y4)
   plot([x1, x2, x3, x4, x5, x1], [y1, y2, y3, y4, y5, y1], color=:blue, lw=0.5)
   segment(x2, y2, x5, y5, :green)
   segment(x1, y1, x3, y3, :green)
   circle(0, 0, R, :orange)
   r1 = sqrt(5)*R/4
   ox1 = sqrt(2)*R*(-sqrt(sqrt(5) + 5) + sqrt(5 - sqrt(5)))/8
   circle(ox1, y1 - r1, r1, :magenta)
   (r2, ox2) = (0.263932022500210, -0.138757275712888)
   circle(ox2, y1 + r2, r2, :brown)
   (r3, ox3) = (0.190983005625053, 0.363271264002680)
   circle(ox3, y1 + r3, r3, :red)
   (a, xa) = (0.440371669705066, 0.807971087145006)
   rect(xa - a, y1-a, xa, y1, :purple)
   factor = 2/1.902113032590307
   println("大円の直径 = $(r1*2factor)")
   println("中円の直径 = $(r2*2factor)")
   println("小円の直径 = $(r3*2factor)")
   println("正方形の一辺の長さ = $(a*factor)")
   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, R, " R", :red, :left, :bottom, delta=delta/2)
       point(xa, y1, "(xa,y1)", :purple, :right, :bottom, delta=delta/2)
       point(xa - a, y1-a, " (xa-a,y1-a)", :purple, :left, :bottom, delta=delta/2)
       point(ox1, y1 - r1, "大円:r1,(ox1,y1-r1)", :magenta, :center, delta=-delta/2)
       point(ox2, y1 + r2, "中円:r2\n(ox2,y1+r2)", :brown, :center, delta=-delta/2)
       point(ox3, y1 + r3, "小円:r3\n(ox3,y1+r3)", :red, :center, delta=-delta/2)
   end
end;

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

今流行りの ChatGPT に算額問題を解かせてみた

2024年02月23日 | Julia

算額(その682)を解かせてみました。

ChatGPT は JPEG や PNG ファイルを読む能力はまだないので,文章でくどく説明しなければならないようです。

ChatGPT 3.5

User
原点を中心とする円を描きます。次に,5個の円を描き足しますが,次のような条件があります。
5個の円は互いに外接します。
それぞれの円はすべて最初に述べた原点を中心とする円に外接します。
最初の円(A とします)の半径を r1,中心座標を (0, 0) とします。
5個の円の半径は全部同じで r2 とします。
5個の円のうちの一つの円(B とします)の中心座標を (x2, y2)とします。
5個の円のうちのもう一つの円(C とします)の中心座標を (0, r1+r2)とします。
すなわち,円C はy軸上にあり,円Aに外接しています。
この状態で,円B が円A と円C に同時に外接するときの中心座標 (x2, y2) を求めてください。

ChatGPT
まず、円Aと円Cが外接する条件から、円Aの半径r1と円Cの中心座標(0, r1+r2)を考慮して、円Cの半径r2は以下のように求められます。

円Aの半径 + 円Cの半径 = 円Aと円Cの中心間の距離

r1 + r2 = √((0 - 0)^2 + (r1 + r2 - r1)^2)
r1 + r2 = √(0 + r2^2)
r1 + r2 = r2

したがって、r1 = 0 となります。

次に、円Bが円Aと円Cに同時に外接するという条件を考慮して、円Bの中心座標(x2, y2)を求めます。

円Bが円Aに外接するための条件から、円Aの中心から円Bの中心までの距離はr1 + r2となります。また、円Bが円Cに外接するための条件から、円Cの中心から円Bの中心までの距離もr1 + r2となります。

これらの条件から、円Bの中心座標(x2, y2)を求めることができます。

円Aと円Bの中心間の距離の式は以下のようになります。

√((0 - x2)^2 + (0 - y2)^2) = r1 + r2

円Cと円Bの中心間の距離の式は以下のようになります。

√((0 - x2)^2 + ((r1 + r2) - y2)^2) = r1 + r2

これらの2つの式を解くことで、円Bの中心座標(x2, y2)を求めることができます。

User
そもそも円なので,円A の半径 r1 が 0というのは間違いです。なので,その後の答えも間違いです。

ChatGPT
申し訳ありません、誤解がありました。では、円Aの半径をr1として、円Cの中心座標を(0, r1+r2)として考えましょう。

円Aと円Cが外接する条件より、円Aの半径r1と円Cの半径r2を用いて次の式が成り立ちます。

r1 + r2 = 円Aと円Cの中心間の距離

円Aと円Bが外接する条件から、円Aの中心から円Bの中心までの距離もr1 + r2になります。同様に、円Bと円Cが外接する条件から、円Cの中心から円Bの中心までの距離もr1 + r2になります。

これらの条件を用いて、円Bの中心座標(x2, y2)を求めることができます。

円Aと円Bの中心間の距離の式は以下のようになります。

√((0 - x2)^2 + (0 - y2)^2) = r1 + r2

円Bと円Cの中心間の距離の式は以下のようになります。

√((0 - x2)^2 + ((r1 + r2) - (r1 + r2))^2) = r1 + r2

これらの2つの式を解くことで、円Bの中心座標(x2, y2)を求めることができます。

User
中心間の距離の計算が違います。

ChatGPT
申し訳ありません、誤りがありました。中心間の距離を正しく計算しましょう。

円Aの中心を原点とすると、円Bの中心座標を(x2, y2)とすると、円Bの中心から円Aの中心までの距離は、円Aの半径r1に円Bの半径r2を加えたものになります。

同様に、円Bの中心から円Cの中心までの距離も、円Cの半径r2になります。

したがって、円Aと円Bの中心間の距離と円Bと円Cの中心間の距離はそれぞれ次のように表されます。

円Aと円Bの中心間の距離:

sqrt((x2 - 0)^2 + (y2 - 0)^2) = r1 + r2

円Bと円Cの中心間の距離:

sqrt((x2 - 0)^2 + (y2 - (r1 + r2))^2) = r2 + r2

これらの式を解いて、円Bの中心座標(x2, y2)を求めることができます。

User
正解です。r1 = 8.775 のとき,r2 はいくつになりますか?

ChatGPT
円Aの半径がr1 = 8.775の場合、円Aと円Cが外接しているので、円Cの中心座標は(0, r1 + r2)となります。

円Cの中心座標から円Aの中心座標までの距離は円Aの半径r1と円Cの半径r2の和に等しいので、以下の式が成り立ちます。

r1 + r2 = sqrt((0 - 0)^2 + (r1 + r2 - r1)^2)

これを解くと、

8.775 + r2 = sqrt(r*2^2)
8.775 + r2 = r2

よって、

r2=8.775 となります。

User
r2を求めるには円Aと円Bの中心間距離を考えないといけません。

ChatGPT
申し訳ありません、その通りです。r2を求めるには、円Aと円Bの中心間の距離を考慮する必要があります。

円Aと円Bが外接するための条件から、円Aの中心から円Bの中心までの距離はr1 + r2に等しいと述べました。

sqrt((x2 - 0)^2 + (y2 - 0)^2) = r1 + r2

ここで、r1は与えられているので、r2を求めるためには、x2とy2の値が必要です。しかし、x2とy2の値が与えられていません。したがって、現時点ではr2を具体的に求めることはできません。

ただし、r1が8.775であるとき、r2を求める必要がある場合には、x2とy2の値が提供される必要があります。それらの値が提供されると、上記の式を使用してr2を求めることができます。

User
x2 と y2は,前に求めましたよ。

ChatGPT
申し訳ありません、先に述べられた条件に従って、円Bの中心座標(x2, y2)を使用することで、r2を求めることができます。

円Bの中心座標が(x2, y2)であるという条件を考慮に入れると、円Aと円Bの中心間の距離と円Bと円Cの中心間の距離は以下のように表されます。

円Aと円Bの中心間の距離:

sqrt((x2 - 0)^2 + (y2 - 0)^2) = r1 + r2

円Bと円Cの中心間の距離:

sqrt((x2 - 0)^2 + (y2 - (r1 + r2))^2) = r2 + r2

これらの式を解いて、r2を求めることができます。既に(x2, y2)が与えられている場合、これらの式を解いてr2を求めることができます。

User
正解です。ちょっと回り道をしましたけど。

ChatGPT
申し訳ありません、回り道をしてしまいましたが、最終的に正しい答えにたどり着けて良かったです。これで問題が解決し、r2を求めることができました。もし他に何か質問があれば、お気軽にどうぞ。

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

算額(その719)

2024年02月22日 | Julia

算額(その719)

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

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

直角三角形内に斜線を引き一方の区画に内接する円を描く。直角を挟む二辺の短い方(鈎)と長い方(股)の辺の長さがそれぞれ 8 寸,15 寸,円の直径が 4 寸のとき,斜線の長さはいかほどか。

鈎,股をそれぞれ「鈎」,「股」
斜線と股の交点座標を (a, 0)
円の半径と中心座標を r, (x, r)
とおき,以下の連立方程式を解く。
斜線の長さ length は後で計算してもよいが,一緒に求めておく。

include("julia-source.txt");

using SymPy
@syms r::positive, x::positive, a::positive,
     length::positive,
     鈎::positive, 股::positive
eq1 = dist(0, 鈎, a, 0, x, r) - r^2
eq2 = dist(0, 鈎, 股, 0, x, r) - r^2
eq3 = sqrt(鈎^2 + a^2) - length
res = solve([eq1, eq2, eq3], (a, x, length));
("a:", res[3][1]) |> println  # a
("x:", res[3][2]) |> println  # x
("length:", res[3][3]) |> println  # length

   ("a:", (-2*r^2*股 - 2*r^2*sqrt(股^2 + 鈎^2) + 2*r*股*鈎 + 2*r*鈎*sqrt(股^2 + 鈎^2) - 股*鈎^2)/(鈎*(2*r - 鈎)))
   ("x:", (-r*sqrt(股^2 + 鈎^2) - 股*(r - 鈎))/鈎)
   ("length:", sqrt(鈎^4*(2*r - 鈎)^2 + (2*r^2*股 + 2*r^2*sqrt(股^2 + 鈎^2) - 2*r*股*鈎 - 2*r*鈎*sqrt(股^2 + 鈎^2) + 股*鈎^2)^2)/(鈎*Abs(2*r - 鈎)))

4 組の解が得られるが,3 番目のものが適解である。
sqrt(股^2 + 鈎^2) を 「弦」とおくと少しきれいになるか。
a: (-2*r^2*股 - 2*r^2*弦 + 2*r*股*鈎 + 2*r*鈎*弦 - 股*鈎^2)/(鈎*(2*r - 鈎))
x: (-r*弦 - 股*(r - 鈎))/鈎
length: sqrt(鈎^4*(2*r - 鈎)^2 + (2*r^2*股 + 2*r^2*弦 - 2*r*股*鈎 - 2*r*鈎*弦 + 股*鈎^2)^2)/(鈎*Abs(2*r - 鈎)))

鈎,股,円の半径が 8,15,4/2 のとき斜線の長さは 10 寸である。
他のパラメータは a = 6,x = 7 である。

(鈎, 股) = (8, 15)
r = 4//2
Abs = abs
((-2*r^2*股 - 2*r^2*sqrt(股^2 + 鈎^2) + 2*r*股*鈎 + 2*r*鈎*sqrt(股^2 + 鈎^2) - 股*鈎^2)/(鈎*(2*r - 鈎)), (-r*sqrt(股^2 + 鈎^2) - 股*(r - 鈎))/鈎, sqrt(鈎^4*(2*r - 鈎)^2 + (2*r^2*股 + 2*r^2*sqrt(股^2 + 鈎^2) - 2*r*股*鈎 - 2*r*鈎*sqrt(股^2 + 鈎^2) + 股*鈎^2)^2)/(鈎*Abs(2*r - 鈎)))

   (6.0, 7.0, 10.0)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, r) = (8, 15, 4/2)
   Abs = abs
   (a, x, length) = ((-2*r^2*股 - 2*r^2*sqrt(股^2 + 鈎^2) + 2*r*股*鈎 + 2*r*鈎*sqrt(股^2 + 鈎^2) - 股*鈎^2)/(鈎*(2*r - 鈎)), (-r*sqrt(股^2 + 鈎^2) - 股*(r - 鈎))/鈎, sqrt(鈎^4*(2*r - 鈎)^2 + (2*r^2*股 + 2*r^2*sqrt(股^2 + 鈎^2) - 2*r*股*鈎 - 2*r*鈎*sqrt(股^2 + 鈎^2) + 股*鈎^2)^2)/(鈎*Abs(2*r - 鈎)))
   @printf("斜線の長さ = %g;  a = %g;  x = %g\n", length, a, x)
   plot([0, 股, 0,  0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   circle(x, r, r)
   segment(0, 鈎, a, 0, :green)
   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, 鈎, " 鈎", :blue, :left, :bottom)
       point(股, 0, " 股", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, "a ", :blue, :right, :bottom, delta=delta/2)
       point(x, r, "r,(x,r)", :red, :center, delta=-delta/2)
   end
end;

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

算額(その718)

2024年02月22日 | Julia

算額(その718)

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

長方形内に菱形と大円 2 個,小円 4 個をいれる。小円は長方形の二辺に内接し菱形と大円に外接している。
長方形の短辺が寸五分のとき,小円の直径はいかほどか。

長方形の中心を原点として,長辺と短辺を 2a, 2b とする。
大円の半径と中心座標を r1, (0, r1)
小円の半径と中心座標を r2, (a - r1, b - r2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     r1::positive, r2::positive
r1 = b/2
eq1 = (a - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = dist(0, b, a, 0, a - r2, b - r2) - r2^2
res = solve([eq1, eq2], (r2, a))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (b*(9 - sqrt(17))/16, b*(5 + 3*sqrt(17))/16)

「問」では「平寸五分」と一文字が欠損しているがこれは「二」であることがわかる。
小円の直径は 2(短辺/2*(9 - √17)/16 = 0.7620147459972405 つまり 「七分六厘二毛有奇」

短辺 = 2.5  # 二寸五分 = 2b
2(短辺/2*(9 - √17)/16)

   0.7620147459972405

「術」は,「から十七の平方根を引いたものに,短辺を十六で割ったものを掛ける」とあるが,五は九の誤記である。

("誤", (5 - sqrt(17))*(短辺/16)) |> println
("正", (9 - sqrt(17))*(短辺/16)) |> println

   ("誤", 0.13701474599724053)
   ("正", 0.7620147459972405)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = 2.5/2
   r1 = b/2
   (r2, a) = (b*(9 - sqrt(17))/16, b*(5 + 3*sqrt(17))/16)
   @printf("小円の直径 = %g\n", 2r2)
   @printf("r2 = %g, a = %g\n", r2, a)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:gray70, lw=0.5)
   plot!([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:green, lw=0.5)
   circle(0, r1, r1)
   circle(0, -r1, r1)
   circle4(a - r2, b - r2, r2, :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(0, r1, " 大円:r1\n (0,r1)", :red, :left, :vcenter)
       point(a - r2, b - r2, " 小円:r2\n (a-r2,b-r2)", :blue, :center, delta=-delta/2)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その717)

2024年02月22日 | Julia

算額(その717)

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

外円内に正三角形と弦を描き,区画された領域に菱形 1 個,甲円 1 個,乙円 3 個,丙円 2 個を入れる。
菱形の長い方の対角線が 1 寸 5 分のとき,丙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, r1 - R)
乙円の半径と中心座標を r2, (x2, y + r2), (0, y + r2)
丙円の半径と中心座標を r3, (x3, 2r1 - R - r3)
菱形の長い方の対角線の長さを 2x
乙円が接している弦と y 軸との交点座標を (0, y) とする。
R = 4r1, r3 = 3r1/4, x3 = 2sqrt(r1*r3) は容易にわかるので,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, x::positive, y::positive,
     r1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive
@syms d
R = 4r1
r3 = 3r1/4
x3 = 2*sqrt(r1*r3)
eq1 = x3^2 + (r3 - r1)^2 - (r1 + r3)^2
eq2 = x3^2 + (2r1 - R - r3)^2 - (R - r3)^2
eq3 = x2^2 + (y + r2)^2 - (R - r2)^2
eq3 = eq3 |> expand
eq4 = dist(sqrt(Sym(3))/2*R, -R/2, 0, R, x2, y + r2) - r2^2
eq4 = eq4 |> expand
eq5 = dist(sqrt(Sym(3))/2*R, -R/2, 0, R, 0, y + r2) - r2^2
eq5 = apart(eq5, d);
eq6 = ((y + 2r1 - R)/2 - 2r1 + R)/(sqrt(Sym(3))/2*R - x) - sqrt(Sym(3))
eq6 = eq6 |> factor;
eq7 = 16R*r3 - 3R^2;

res = solve([eq3, eq4, eq5, eq6], (r1, r2, x2, y))

   1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (25*sqrt(3)*x/111, 8*sqrt(3)*x/37, 32*x/37, 28*sqrt(3)*x/111)

丙円の直径は,甲円の直径の 3/4 倍であり,甲円の直径は 2(25*sqrt(3)*x/111) である
菱形の長い方の対角線の長さが 2x のとき,丙円の直径は 2(25*sqrt(3)*x/111)*(3/4) である。

2x = 1.5 寸のとき 0.4388642248907628 すなわち,「4分3厘8毛有奇」である。

x = 1.5/2
2(25*sqrt(3)*x/111)*(3/4)

   0.4388642248907628

「術」は,「3 の平方根に菱長と 25 を掛けて 148 で割る」sqrt(3)*1.5*25/148 = 0.4388642248907628 である。

sqrt(3)*1.5*25/148

   0.4388642248907628

その他のパラメータは以下のとおりである。

R = 1.1703, r1 = 0.292576;  r2 = 0.280873;  x2 = 0.648649;  r3 = 0.219432;  x3 = 0.506757

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x = 15//20
   (r1, r2, x2, y) = (25*sqrt(3)*x/111, 8*sqrt(3)*x/37, 32*x/37, 28*sqrt(3)*x/111)
   R = 4r1
   r3 = 3r1/4
   x3 = 2*sqrt(r1*r3)
   @printf("丙円の直径 = %g\n", 2r3)
   @printf("R = %g, r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", R, r1, r2, x2, r3, x3)
   plot([√3/2*R, 0, -√3/2*R, √3/2*R], [-R/2, R, -R/2, -R/2], color=:gray70, lw=0.5)
   circle(0, 0, R)
   circle(0, r1 - R, r1)
   circle(x3, 2r1 - R - r3, r3, :blue)
   circle(-x3, 2r1 - R - r3, r3, :blue)
   circle(0, y + r2, r2, :orange)
   circle(x2, y + r2, r2, :orange)
   circle(-x2, y + r2, r2, :orange)
   plot!([x, 0, -x, 0, x], [(y + 2r1 - R)/2, y, (y + 2r1 - R)/2, 2r1 - R, (y + 2r1 - R)/2], color=:green, lw=0.5)
   xa = sqrt(R^2 - y^2)
   segment(-xa, y, xa, y)
   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(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(0, y, " y", :black, :left, :top, delta=-delta)
       point(0, y + r2, " 乙円:r2\n (0,y+r2)", :black, :left, :vcenter)
       point(x2, y + r2, " 乙円:r2\n (x2,y+r2)", :black, :center, delta=-delta/2)
       point(x, (y + 2r1 - R)/2, "(x,(y+2r1-R)/2) ", :green, :right, :vcenter)
       point(0, 2r1 - R, " 2r1-R", :red, :center, :bottom, delta=delta/2)
       point(0, r1 - R, "甲円:r1\n(0,r1-R)", :red, :center, :bottom, delta=delta/2)
       point(x3, 2r1 - R - r3, "丙円:r1\n(x3,2r1-R-r3)", :black, :center, :bottom, delta=delta/2)
   end
end;

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

算額(その716)

2024年02月22日 | Julia

算額(その716)

埼玉県秩父市大宮 秩父神社 明治20年

山口正義(2015): やまぶき, 第27号
https://yamabukiwasan.sakura.ne.jp/ymbk27.pdf

楕円の中に大円4個,小円3個が入っている。小円の直径が 1 寸のとき,楕円の長径はいかほどか。

楕円の長半径,短半径を a, b,中心座標を (0, 0); b = 3r2
大円の半径と中心座標を r1, (x1, 0), (0, r2); r1 = 2r2
楕円と大円の接点座標を (x0, y0)
小円の半径を r2
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, r2::positive,
     a::positive, b::positive,
     x0::positive, y0::positive
r1 = 2r2
b = 3r2
eq1 = x1^2 + r2^2 - 4r1^2
eq2 = (x0 - x1)^2 + y0^2 - r1^2
eq3 = x0^2/a^2 + y0^2/b^2 - 1
eq4 = -b^2*x0/(a^2*y0) - (x1 - x0)/y0
res = solve([eq1, eq2, eq3, eq4], (a, x1, x0, y0))

   1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (6*r2, sqrt(15)*r2, 4*sqrt(15)*r2/3, sqrt(21)*r2/3)

楕円の長半径は小円の半径の 6 倍である。
小円の半径が 1/2 寸のとき,長半径は 3 寸(長径は 6 寸)である。

r2 = 1//2
(6.0*r2, sqrt(15)*r2, 4*sqrt(15)*r2/3, sqrt(21)*r2/3)

   (3.0, 1.9364916731037085, 2.581988897471611, 0.7637626158259733)

楕円の長径 = 6;  a = 3;  x1 = 1.93649;  x0 = 2.58199;  y0 = 0.763763

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   r1 = 2r2
   b = 3r2
   (a, x1, x0, y0) = (6*r2, sqrt(15)*r2, 4*sqrt(15)*r2/3, sqrt(21)*r2/3)
   @printf("楕円の長径 = %g;  a = %g;  x1 = %g;  x0 = %g;  y0 = %g\n", 2a, a, x1, x0, y0)
   plot()
   circle(x1, 0, r1)
   circle(-x1, 0, r1)
   circle.(0, [0, 2r2, -2r2], r2)
   circle(0, r2, r1, :green)
   circle(0, -r2, r1, :green)
   ellipse(0, 0, a, b, color=: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(x0, y0, "(x0,y0)", :blue, :left, :bottom, delta=delta/2)
       point(x1, 0, "x1", :magenta, :center, :bottom, delta=delta/2)
       point(0, r2, " r2", :green, :center, delta=-delta/2)
       point(0, b, " b", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その715)

2024年02月21日 | Julia

算額(その715)

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード:円6個,外円,楕円2個

外円内に長径が互いに直行する 2 つの楕円がある。中央に大円,その周りに 4 個の小円が入っている。
外円と大円の直径がそれぞれ 10 寸,3 寸のとき,小円の直径はいかほどか。

計算を簡単にするために図形を45度回転させたものを考える。
外円の半径と中心座標を R, (0, 0 )
大円の半径と中心座標を r1, (0, 0)
小円の半径と中心座標を r2, (r1 + r2, 0)
とおく。
楕円の長半径と短半径は a = R, b = r1 である。
楕円と小円の接点座標を (x0, y0)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, a::positive, b::positive,
     r1::positive, r2::positive,
     x0::positive, y0::positive
(a, b) = (R, r1)
eq1 = (x0 - r1 - r2)^2 + y0^2 - r2^2
eq2 = -b^2*x0/(a^2*y0) - (r1 + r2 - x0)/y0
eq3 = x0^2/a^2 + y0^2/b^2 - 1
res = solve([eq1, eq2, eq3], (r2, x0, y0))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (r1 - 2*r1^3/R^2, 2*r1, r1*sqrt(R^2 - 4*r1^2)/R)

小円の直径は 2(r1 - 2*r1^3/R^2) である。
外円,大円の半径がそれぞれ 10 寸,3 寸のとき,小円の直径は 2.46 寸である。

(R, r1) = (10, 3) ./ 2
2(r1 - 2*r1^3/R^2)  # r2

   2.46

2*r1, r1*sqrt(R^2 - 4*r1^2)/R  # x0, y0

   (3.0, 1.2)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1) = (10, 3) .// 2
   (a, b) = (R, r1)
   (r2, x0, y0) = (r1 - 2*r1^3/R^2, 2*r1, r1*sqrt(R^2 - 4*r1^2)/R)
   @printf("小円の直径 = %g;  r2 = %g;  x0 = %g;  y0 = %g\n", 2r2, r2, x0, y0)
   plot()
   circle(0, 0, R)
   circle(0, 0, r1)
   circle42(0, r1 + r2, r2, :green)
   ellipse(0, 0, R, r1, color=:blue)
   ellipse(0, 0, r1, R, color=: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(x0, y0, "(x0,y0)", :blue, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(0, r1, " r1", :red, :left, delta=-delta/2)
       point(0, 0, "大円:r1\n(0,0)", :red, :center, :vcenter)
       point(r1 + r2, 0, "小円:r2\n(r1+r2,0)", :green, :center, :vcenter)
   end
end;

 

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

算額(その714)

2024年02月21日 | Julia

算額(その714)

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

外円に内接する正六角形,正三角形があり,甲円,乙円 3 個ずつを入れる。
ただし,甲円は外円に内接し,正三角形,正六角形に外接する。また,甲円は正三角形に外接し,正六角形に内接する。
乙円の直径が 1.5 寸のとき,外円の直径はいかほどか。

外円の直径を 1 として正三角形,正六角形,甲円,乙円を定め,乙円の直径と外円の直径の比に基づいて乙円の直径が 1.5 寸のときの外円の直径を求める。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (x2, y2)
として以下の連立方程式を解く。

まず,半径 R = 1/2 の外円と,外円に内接する正六角形に内接し,正三角形に外接する円(甲円)の半径 (r1) と中心座標 (x1, y1) を求める。

include("julia-source.txt");

using SymPy
@syms x1::positive, y1::positive, r1::positive, R::positive
s3 = sqrt(Sym(3))
R = 1//2
eq0 = s3*R/2 - r1 - x1
eq1 = x1^2 + y1^2 - (R - r1)^2
eq2 = dist(0, R, s3*R//2, -R//2, x1, y1) - r1^2
res1 = solve([eq0, eq1, eq2], (x1, y1, r1))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (-11*sqrt(3)/4 - sqrt(73 - 40*sqrt(3)) + 7, -949*sqrt(219 - 120*sqrt(3))/529 - 520*sqrt(3)*sqrt(219 - 120*sqrt(3))/529 - 5/4 + sqrt(3) + 880*sqrt(3)*sqrt(73 - 40*sqrt(3))/529 + 1606*sqrt(73 - 40*sqrt(3))/529, -7 + sqrt(73 - 40*sqrt(3)) + 3*sqrt(3))

res1[1][1].evalf(), res1[1][2].evalf(), res1[1][3].evalf()

   (0.308657048910078, 0.214101615137755, 0.124355652982141)

半径 R = 1 の外円に内接する正六角形に内接し,正三角形に外接する円(乙円)の半径 (r2) と中心座標 (x2, y2) を求める。

using SymPy
@syms x2, y2, r2, R
s3 = sqrt(Sym(3))
R = 1//2
eq0 = s3*R//2 - r2 - x2
eq1 = (R//2 + r2)*cosd(Sym(30)) - x2
eq2 = (R//2 + r2)*sind(Sym(30)) - y2
res2 = solve([eq0, eq1, eq2], (x2, y2, r2))

   Dict{Any, Any} with 3 entries:
     y2 => -1/4 + sqrt(3)/4
     x2 => 3/4 - sqrt(3)/4
     r2 => -3/4 + sqrt(3)/2

(3/4 - sqrt(3)/4, -1/4 + sqrt(3)/4, -3/4 + sqrt(3)/2)

   (0.3169872981077807, 0.1830127018922193, 0.1160254037844386)

外円の半径が R = 1/2 なら乙円の半径は 0.116025403784439 になる。
乙円の半径が 1.5/2 なら外円の直径は 6.464101615137762 寸である。

「答」では 6.5 寸となっている。
「術」も,「乙円の直径を 13 倍して,3 で割る」ということなので 1.5✕13/3 = 6.5 ということになる。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 1//2
   plot([√3, 0, -√3, √3].*R/2, [-1, 2, -1, -1].*R/2, color=:green, lw=0.5)
   plot!([√3, √3, 0, -√3, -√3, 0, √3].*R/2, [-1, 1, 2, 1, -1, -2, -1].*R/2, color=:blue, lw=0.5)
   circle(0, 0, R)
   (x1, y1, r1) = (0.308657048910078, 0.214101615137755, 0.124355652982141)
   rotate(x1, y1, r1, :magenta, angle=60)
   (x2, y2, r2) = (0.316987298107781, 0.183012701892219, 0.116025403784439)
   rotate(x2, y2, r2, :brown, angle=60)
   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, y1, " 甲円:r1,(x1,y1)", :magenta, :left, :vcenter)
       point(x2, y2, " 乙円:r2,(x2,y2)", :brown, :left, :vcenter)
       point(√3R/2, R/2, " (√3R/2,R/2)", :red, :left, :vcenter)
       # plot!(xlims=(0.15, 0.5), ylims=(0.05, 0.35))
   end
end;

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

算額(その713)

2024年02月20日 | Julia

算額(その713)

八六 加須市多聞寺 愛宕神社 明治13年(1880)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

キーワード:円1個,楕円3個,正三角形,正六角形

正六角形の中に正三角形と楕円 3 個が入っている。正三角形の中には直径 1 寸の円が入っている。楕円の短径が 0.44 寸のとき,楕円の長径はいかほどか。

円の半径と中心座標を r, (0, 0)
下にある楕円の長半径と短半径を a, b,中心座標を (0, -r - b)
楕円と正六角形の辺の接点の座標を (x0, y0)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r::positive, a::positive, b::positive,
     x0::positive, y0::negative
eq1 = x0^2/a^2 + (y0 + r + b)^2/b^2 - 1
eq2 = -b^2*x0/(a^2*(y0 + r + b)) - 1/√Sym(3)
eq3 = (2r + y0)/x0 - 1/√Sym(3);
solve([eq1, eq2, eq3], (a, x0, y0))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (sqrt(3)*sqrt(r)*sqrt(-2*b + r), sqrt(3)*r*(2*b - r)/(b - r), r^2/(b - r))

楕円の長半径 a は,r と b の関数である。
sqrt(3r*(r - 2b))
r = 1/2 寸, b = 0.44/2 寸 のとき a = 0.3 である(長径は 0.6 寸)

r = 1/2
b = 0.44/2
sqrt(3r*(r - 2b))

   0.3

接点の座標は (x0, y0) = (0.185577, -0.892857) である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, b) = (1/2, 0.44/2)
   (a, x0, y0) = (sqrt(3)*sqrt(r)*sqrt(-2*b + r), sqrt(3)*r*(2*b - r)/(b - r), r^2/(b - r))
   @printf("長径 = %g;  x0 = %g;  y0 = %g\n", 2a, x0, y0)
   plot([√3, 0, -√3, √3].*r, [-1, 2, -1, -1].*r, color=:green, lw=0.5)
   plot!([√3, √3, 0, -√3, -√3, 0, √3].*r, [-1, 1, 2, 1, -1, -2, -1].*r, color=:blue, lw=0.5)
   circle(0, 0, r)
   l = r + b

   ellipse(0, -l, a, b)
   ellipse(l*cosd(30), l*sind(30), a, b, φ=120)
   ellipse(l*cosd(150), l*sind(150), a, b, φ=240)
   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, -r, "-r", :red, :center, :bottom, delta=delta/2)
       point(0, -l, "-r-b")
       point(x0, y0, "(x0,y0)")
   end
end;

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

算額(その712)

2024年02月20日 | Julia

算額(その712)

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

等脚台形の中に正方形と 3 個の等円が入っている。上底の長さが 3 寸のとき,下底の長さはいかほどか。

上底,下底の長さをそれぞれ 2b, 2a とする。
正方形の一辺の長さを 2c
等円の半径と中心座標を r, (0, 2c + r), (c + r, r)
台形の高さは 2c + 2r
とおいて,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     c::positive, r::positive
h = 2c + 2r
eq1 = (c + r)^2 + 4c^2 - 4r^2
eq2 = dist(a, 0, b, h, 0, 2c + r) - r^2
eq3 = dist(a, 0, b, h, c + r, r) - r^2
res = solve([eq1, eq2, eq3], (a, c, r)) 

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (69*b/5, 9*b/5, 3*b)

下底の長さは上底の長さの 69/5 倍,
正方形の一辺の長さは上底の長さの 9/5 倍,
等円の直径は上底の長さの 3 倍である。

b = 3/2
2b .*(69/5, 9/5, 3)

   (41.400000000000006, 5.4, 9.0)

上底の長さが 3 寸のとき,下底の長さは 41.4 寸である。
「答」では 27.6 寸となっている。
「術」では,「上底の 13.8 倍」と書いている。それに従うと 3✕13.8 = 41.4 寸になる。答えは 2✕13.8 = 27.6 寸 と書いてしまったのであろう。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = 3//2
   (a, c, r) = (69*b/5, 9*b/5, 3*b)
   @printf("上底の長さ = %g;  下底の長さ = %g;  正方形の一辺の長さ = %g;  等円の直径 = %g\n", 2b, 2a, 2c, 2r)
   h = 2c + 2r
   #@printf("r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g\n", r1, x1, r2, x2, r3)
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:green, lw=0.5)
   circle.([1, -1].*(c + r), r, r)
   circle(0, 2c + r, r)
   rect(-c, 0, c, 2c, :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(c, 0, " c", :blue, :left, :bottom, delta=delta/2)
       point(0, 2c, " 2c", :blue, :left, delta=-delta/2)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
       point(b, h, " (b,h)", :green, :center, :bottom, delta=delta)
       point(c + r, r, "(c+r,r)", :red, :center, delta=-delta/2)
       point(0, 2c + r, " 2c+r", :red, :left, :vcenter)
   end
end;

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

算額(その711)

2024年02月20日 | Julia

算額(その711)

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年(部分拡大図)(加須市)
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

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

直線の上に甲円,乙円が 2 個ずつ,丙円が 1 個載っている。丙円の上に正方形が 1 個 載っている。丙円の直径が 1 寸のとき,甲円の直径はいかほどか。

甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (0, r3)
とおき,以下の連立方程式を解く。
甲円と乙円の中心を結ぶ直線は x 軸と角度 22.5° で交わる。
tan(22.5°) =  tan(PI/8) は,思いの外きれいな数 √2 - 1 である。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive,
     r2::positive, x2::positive,
     r3::positive
eq1 = r2/(x2 + 2r3) - tan(PI/8)
eq2 = r1/(x1 + 2r3) - tan(PI/8)
eq3 = x1 - x2 - 2sqrt(r1*r2)
eq4 = x2 - 2sqrt(r2*r3)
res = solve([eq1, eq2, eq3, eq4], (r1, x1, r2, x2)) 

   1-element Vector{NTuple{4, Sym}}:
    (2*r3*(-4*sqrt(2) - 2*sqrt(4 - 2*sqrt(2)) + 4*sqrt(2 - sqrt(2)) + 7), 2*r3*(-2 + 2*sqrt(4 - 2*sqrt(2)) + 3*sqrt(2)), 2*r3, 2*sqrt(2)*r3)

甲円の半径は 乙円の半径の 2*(-4*sqrt(2) - 2*sqrt(4 - 2*sqrt(2)) + 4*sqrt(2 - sqrt(2)) + 7) = 14 - 8√2 - 4sqrt(4 - 2√2) + 8sqrt(2 - √2) 倍である。

乙円の直径が 1 寸のとき,甲円の直径は 4.479657617687099,4寸4分7厘9毛有奇である。

14 - 8√2 - 4sqrt(4 - 2√2) + 8sqrt(2 - √2)

   4.479657617687099

その他のパラメータは以下のとおりである。正方形の大きさは任意。

r1 = 2.23983;  x1 = 4.40743;  r2 = 1;  x2 = 1.41421;  x3 = 0.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1//2
   a = r3*12 / √2
   (r1, x1, r2, x2) = (2*r3*(-4*sqrt(2) - 2*sqrt(4 - 2*sqrt(2)) + 4*sqrt(2 - sqrt(2)) + 7), 2*r3*(-2 + 2*sqrt(4 - 2*sqrt(2)) + 3*sqrt(2)), 2*r3, 2*sqrt(2)*r3)
   @printf("r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g\n", r1, x1, r2, x2, r3)
   plot([0, 1, 0, -1, 0] .* a, [0, 1, 2, 1, 0] .* a .+ 2r3)
   circle(x1, r1, r1)
   circle(-x1, r1, r1)
   circle(x2, r2, r2, :orange)
   circle(-x2, r2, r2, :orange)
   circle(0, r3, r3, :blue)
   abline(-2r3, 0, tand(45/2), -2.5r3, x1+r1) 
   abline(-2r3, 0, 1, -2.5r3, 1.2x1) 
   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(-2r3, 0, "-2r3", :black, :left, delta=-delta/2)
       point(0, r3, " r3", :blue, :left, delta=-delta/2)
       point(x2, r2, " 乙円:r2\n (x2,r2)", :black, :left, delta=-delta/2)
       point(x1, r1, " 甲円:r1\n (x1,r1)", :red, :left, delta=-delta/2)
   end
end;

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

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

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