裏 RjpWiki

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

算額(その1293)

2024年09月14日 | Julia

算額(その1293)

十七 群馬県高崎市八幡町 八幡宮 文化7年(1810)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円5個,菱形

菱形の中に,大円 2 個,中円 2 個,小円 1 個を容れる。菱形の対角線の長い方が 12 寸,短い方が 5 寸のとき,小円の直径はいかほどか。

菱形の対角線を 2a, 2b; a > b
大円の半径と中心座標を r1, (r3 + r1, 0)
中円の半径と中心座標を r2, (0, r3 + r2)
小円の半径と中心座標を r3, (0, 0)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     r1::positive, r2::positive, r3::positive;
sinθ = b/sqrt(a^2 + b^2)
cosθ = a/sqrt(a^2 + b^2)
eq1 = (r3 + r1)^2 + (r3 + r2)^2 - (r1 + r2)^2 |> expand
eq2 = (a - r3 - r1)*sinθ - r1
eq3 = (b - r3 - r2)*cosθ - r2;

SymPy では能力的に一度に解けないので,逐次解いていく。
まず,eq1 を解いて r1 を求める。

ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println

   r3*(r2 + r3)/(r2 - r3)

eq2 の r1 に,上で求めた r1 を代入し,r2 を求める。

eq12 = eq2(r1 => ans_r1)
ans_r2 = solve(eq12, r2)[1]
ans_r2 |> println

   r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))

eq3 に,上で求めた r1, r2 を代入し,r3 を求める。

eq13 = eq3(r1 => ans_r1, r2 => ans_r2) |> simplify |> numerator
ans_r3 = solve(eq13, r3)[2] |> factor
ans_r3 |> println

   -a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2

解として得られた r3 を 代入して r2,r1 を確定する。

a = 12/2
b = 5/2
r3 = -a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2
r2 = r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))
r1 = r3*(r2 + r3)/(r2 - r3)
(r1, r2, r3) |> println

   (1.5296184351192643, 0.9631806558860881, 0.49337363357064845)

菱形の対角線の長い方が 12 寸,短い方が 5 寸のとき,小円の直径は 2*0.49337363357064845 = 0.9867472671412969 寸である。「答」では「小円の径一寸三分二厘半微弱」としているが,正しいとは思えない。

すべてのパラメータは以下のとおりである。

   a = 6;  b = 2.5;  r1 = 1.52962;  r2 = 0.963181;  r3 = 0.493374

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = -a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2
   r2 = r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))
   r1 = r3*(r2 + r3)/(r2 - r3)
   @printf("菱形の対角線の長い方が %g,短い方が %g のとき,小円の直径は %g である。\n", 2a, 2b, 2r3)
   @printf("a = %g;  b = %g;  r1 = %g;  r2 = %g;  r3 = %g\n", a, b, r1, r2, r3)
   plot([a , 0, -a, 0, a], [0, b, 0, -b, 0], color=:green, lw=0.5)
   circle2(r3 + r1, 0, r1)
   circle22(0, r3 + r2, r2, :blue)
   circle(0, 0, r3, :magenta)
   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, 0, " a", :green, :left, :vcenter)
       point(0, b, "b", :green, :center, :bottom, delta=delta)
       point(r3 + r1, 0, "大円:r1\n(r3+r1,0)", :red, :center, delta=-delta)
       point(0, r3 + r2, "中円:r2\n(0,r3+r2)", :blue, :center, :bottom, delta=delta)
       point(0, 0, "小円:r3,(0,0)", :magenta, :right, :vcenter, deltax=-7delta)
   end  
end;

draw(12/2, 5/2, true)

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

算額(その1292)

2024年09月14日 | Julia

算額(その1292)

百四十四 群馬県勢多郡宮城村柏倉 諏訪神社 (1914)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円6個,直線上

直線上にある大円 3 個が交わっている。隙間に 3 個の小円を容れる。大円の直径が与えられたとき,小円の直径を得る術を述べよ。

大円の半径と中心座標を r1, (0, 0), (r1 - r2, 0)
小円の半径と中心座標を r2, (0, 0), (0, r1 - r2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

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

eq = 2(r1 - r2)^2 - (r1 + r2)^2
res = solve(eq, r2)[1]  # 1 of 2
res |> println

   r1*(3 - 2*sqrt(2))

小円の半径 r2 は,大円の半径 r1 の (3 - 2√2) 倍である。
大円の直径が 1 寸のとき,小円の直径は (3 - 2√2) = 0.1715728752538097 寸である。

function draw(r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = r1*(3 - 2√2)
   @printf("大円の直径が %g のとき,小円の直径は %g である。\n", 2r1, 2r2)
   plot()
   circle(0, 0, r1)
   circle2(r1 - r2, 0, r1)
   circle(0, 0, r2, :blue)
   circle22(0, r1 - 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(r1 - r2, 0, "大円:r1\n(r1-r2,0)", :red, :right, :vcenter, deltax=-delta)
       point(0, r1 - r2, "小円:r2,(0,r1-r2)", :blue, :center, :bottom, delta=7delta)
       point(r2, 0, "r2 ", :blue, :right, :vcenter)
   end  
end;

draw(1/2, true)

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

算額(その1291)

2024年09月14日 | Julia

算額(その1291)

百四十一 群馬県藤岡市鮎川 北野神社 明治24年(1891)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円1個,直角三角形,正方形

直角三角形と,その鈎(直角を挟む 2 辺のうち,短い方)を共有する正方形で区分される領域に黄円を描く。弦(直角三角形の斜辺)の長さが 6.5 寸のとき,黄円が最大になるのは鈎がどれほどのときか。

鈎が 0 に近づく場合,黄円の直径も 0 に近づく。鈎が「弦/√2」 に近づく場合,黄円の直径も 0 に近づく。両者の中ほど(真ん中ではない)のときに黄円の直径が最も大きくなる。

計算の都合上,左右反転させて図を描く。
鈎(正方形の一辺のながさ)を「鈎」
黄円の半径と中心座標を r, (a + r, 0)
股(直角を挟む 2 辺のうち,長い方)を「股」
弦を「弦」
正方形の辺と弦の交点座標を (a, c)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     a::positive, c::positive, r::positive;
a = 鈎
eq1 = 鈎^2 + 股^2 - 弦^2
eq2 = c*股 - 鈎*(股 - a)
eq3 = (股 - a) + c - 弦*c/a - r
eq4 = dist2(股, 0, 0, 鈎, a + r, r, r)
res = solve([eq1, eq2, eq4], (r, 股, c))[1]

   ((-弦*sqrt(弦^2 - 鈎^2)*sqrt(弦^4 - 弦^2*鈎^2 - 2*弦^2*鈎*sqrt(弦^2 - 鈎^2) + 2*鈎^3*sqrt(弦^2 - 鈎^2)) + (弦^2 - 2*鈎^2)*(弦^2 - 鈎^2))/(2*(弦^2 - 鈎^2)^(3/2)), sqrt(弦^2 - 鈎^2), -鈎^2/sqrt(弦^2 - 鈎^2) + 鈎)

黄円の半径 r は,r = ((-弦*sqrt(弦^2 - 鈎^2)*sqrt(弦^4 - 弦^2*鈎^2 - 2*弦^2*鈎*sqrt(弦^2 - 鈎^2) + 2*鈎^3*sqrt(弦^2 - 鈎^2)) + (弦^2 - 2*鈎^2)*(弦^2 - 鈎^2))/(2*(弦^2 - 鈎^2)^(3/2)), sqrt(弦^2 - 鈎^2), -鈎^2/sqrt(弦^2 - 鈎^2) + 鈎) と表すことができる。

弦が 6.5 のとき,鈎と黄円の半径の関係は下図のようである。鈎が 2 〜 3 の範囲で黄円の半径 r が最大になる。

using Plots
pyplot(size=(300, 200), grid=false, aspectratio=:none, label="")
plot(res[1](弦 => 6.5), xlims=(0, 6.5/√2), xlabel="鈎", ylabel="r")
hline!([0])

黄円の半径が最大になるときの鈎の値を求めるには,r の式において,鈎で微分した導関数の値が 0 になるときの鈎の値を求めればよい。

solve(diff(res[1], 鈎), 鈎)[1] |> println

   弦*(-1/2 + sqrt(3)/2)

鈎が 弦*(√3 - 1)/2 のとき,黄円の半径(直径)は最大値を取る。
弦が 6.5 寸の場合には,鈎 = 6.5*(√3 - 1)/2 = 2.3791651245988508 のとき,半径 r は最大値 0.584868958605704 になる。
直径は 1.16973791721141 である。

res[1](弦 => 6.5, 鈎 => 2.3791651245988508) |> println
2res[1](弦 => 6.5, 鈎 => 2.3791651245988508) |> println

   0.584868958605704
   1.16973791721141

function draw(弦, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 鈎 = 弦*(√3 - 1)/2
   (r, 股, c) = ((-弦*sqrt(弦^2 - 鈎^2)*sqrt(弦^4 - 弦^2*鈎^2 - 2*弦^2*鈎*sqrt(弦^2 - 鈎^2) + 2*鈎^3*sqrt(弦^2 - 鈎^2)) + (弦^2 - 2*鈎^2)*(弦^2 - 鈎^2))/(2*(弦^2 - 鈎^2)^(3/2)), sqrt(弦^2 - 鈎^2), -鈎^2/sqrt(弦^2 - 鈎^2) + 鈎)
   @printf("弦(斜線)の長さが %g のとき,鈎(正方形の一辺の長さ)が %g のときに黄円の直径が最大値 %g になる。\n", 弦, 鈎, 2r)
   plot([0, a, a, 0, 0], [0, 0, a, a], color=:green, lw=0.5)
   circle(a + r, r, r)
   segment(0, 鈎, 股, 0, :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, a, "a ", :green, :right, :bottom, delta=3delta)
       point(a, a, "(a,a)", :green, :right, :bottom, delta=delta)
       point(a, 0, "a", :green, :center, delta=-2delta)
       point(股, 0, "股", :blue, :center, delta=-2delta)
       point(a + r, r, "黄円:r\n(a+r,r)", :red, :center, delta=-2delta)
       δ = 3delta
       dimension_line(2delta, 鈎 + 3delta, 股 + 2delta, 3delta, "弦", delta=3delta, deltax=3delta)
       dimension_line(-4delta, 0, -4delta, 鈎, "鈎", deltax=-4delta)
       plot!(xlims=(-15delta, 股 + 3delta), ylims=(-9delta, 鈎 + 3delta))
   end  
end;

draw(6.5, true)

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

算額(その1290)

2024年09月14日 | Julia

算額(その1290)

百四十一 群馬県藤岡市鮎川 北野神社 明治24年(1891)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円4個,楕円

楕円の中に小円 2 個を容れる。小円はそれぞれが楕円と 2 点で接する。さらに,小円と外接する大円 2 個を描く。大円は互いに外接し,2 個の小円とも外接する。楕円の長径と短径,小円の直径が与えられた解き,大円の直径を求める術を述べよ。

楕円の長半径,短半径を a, b
大円の半径と中心座標を r1, (0, r1)
小円の半径と中心座標を r2, (x2, 0)
とおき,以下の連立方程式を解く。

eq1 は「算法助術の公式84」による。

include("julia-source.txt");

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

   (-(-a^2*b^2 + a^2*r2^2 + b^4)/(2*b^2*r2), sqrt(-(a - b)*(-b + r2))*sqrt(a + b)*sqrt(b + r2)/b)

res[1] |>  simplify |> println
res[2] |>  simplify |> println

   (a^2*b^2 - a^2*r2^2 - b^4)/(2*b^2*r2)
   sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2)/b

たとえば,長径 = 6,短径 = 3,小円の直径 = 2 のとき,大円の半径は 1.375 である(直径は 2.75)。

a = 6/2; b = 3/2; r2 = 2/2
(a^2*b^2 - a^2*r2^2 - b^4)/(2*b^2*r2)

   1.375

function draw(a, b, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = (a^2*b^2 - a^2*r2^2 - b^4)/(2*b^2*r2)
   x2 = sqrt(a^2*b^2 - a^2*r2^2 - b^4 + b^2*r2^2)/b
   @printf("長径が %g,短径が %g,小円の直径が %g のとき,大円の直径は %g である。\n", 2a, 2b, 2r2, 2r1)
   plot()
   ellipse(0, 0, a, b, color=:blue)
   circle22(0, r1, r1, :green)
   circle2(x2, 0, r2)
   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,(0,r1)", :green, :center, delta=-delta)
       point(x2, 0, "小円:r2,(x,0)", :red, :center, delta=-delta)
       point(0, b, "b", :blue, :center, :bottom, delta=delta)
       point(a, 0, " a", :blue, :left, :vcenter)
   end  
end;

draw(6/2, 3/2, 2/2, true)

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

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

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