裏 RjpWiki

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

算額(その778)

2024年03月13日 | Julia

算額(その778)

川田保知:『算法極数小補解義』 文化15年戊寅正月
山口正義:やまぶき,第26号

https://yamabukiwasan.sakura.ne.jp/ymbk26.pdf

直角三角形内に斜線(文斜と命名)を入れ,股の一部分を武斜と命名する。
文斜と武斜を 2 寸,7 寸に固定したとき,股は武斜 + x になる。面積が最大になるときの股はいかほどか。

直角三角形の面積は x の関数になる。面積を求める式を微分し,接線の傾きが 0 になるときの x を求める。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 文斜::positive, 武斜::positive, x::positive
S = sqrt(文斜^2 - x^2) * (武斜 + x)//2  # 鈎*股/2
S |> println

   (x/2 + 武斜/2)*sqrt(-x^2 + 文斜^2)

直角三角形の面積 S は x の関数である。

たとえば,文斜,武斜がそれぞれ 2, 7 のとき,面積は x が 0.5 前後で最大になる。

S2 = S(文斜 => 2, 武斜 => 7)
using Plots
pyplot(size=(500, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
p1 = plot(S2, xlims=(0, 2), xlabel="x", ylabel="面積 S")
p2 = plot(S2, xlims=(0.4, 0.6), xlabel="x", ylabel="面積 S", title="左図の拡大図", titlefont=11)
plot!(p1, p2)

接線の傾きが 0 になるときの x は,(sqrt(8*文斜^2 + 武斜^2) - 武斜)/4 である。

山口はまだ術を解読できていないとしているが「術曰置併文斜▢八段武斜▢開平方」の部分は sqrt(8*文斜^2 + 武斜^2) に相当するであろう。二箇所の▢は「乗」と思われる。引き続く「加武斜三段皈」はわからない。特に「皈」が。

f = diff(S, x)
res2 = solve(f, x)[1]
res2 |> println
res2.evalf() |> println

   -武斜/4 + sqrt(8*文斜^2 + 武斜^2)/4
   -0.25*武斜 + 0.707106781186548*(文斜^2 + 0.125*武斜^2)^0.5

文斜,武斜がそれぞれ 2, 7 のとき,面積は x が 0.5 のとき最大になる。
また,そのときの股は,股 = 武斜 + 0.5 = 7.5 である。

x = res2(文斜 => 2, 武斜 => 7)
x |> println

   1/2

文斜,武斜がそれぞれ 2, 7 のとき,x が 0.5 のとき,面積は最大値 7.26184377413891 になる。

文斜が 2,武斜が 7 のとき,股(= 武斜 + x)が 7.5 のとき,直角三角形の面積が最大になる(x = 0.5)。

S = sqrt(文斜^2 - x^2) * (武斜 + x)//2
S |> println
S(x => 0.5, 文斜 => 2, 武斜 => 7).evalf() |> println

   sqrt(文斜^2 - 1/4)*(武斜/2 + 1/4)
   7.26184377413891

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (文斜, 武斜) = (2, 7)
   x = 0.5
   鈎 = sqrt(文斜^2 - x^2)
   股 = 武斜 + x
   @printf("文斜が %g,武斜が %g のとき,股(= 武斜 + x)が %g のとき,直角三角形の面積が最大になる(x = %g)。\n", 文斜, 武斜, 武斜 + x, x)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:red, lw=0.5)
   segment(0, 鈎, x, 0, :blue)
   segment(x, 0, 股, 0, :green, lw=1)
   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, "x+武斜 = 股", :red, :right, delta=-3delta)
       point(0, 鈎, " 鈎", :red, :left, :bottom, delta=delta)
       point(x, 0, " x", :blue, :left, :bottom, delta=delta)
       point(x, 鈎/2, "文斜", :blue, :left, mark=false)
       point(股/2, 0, "武斜", :green, :left, :bottom, delta=2delta, mark=false)
       plot!(ylims=(-0.5, 2.1))
   end
end;

 

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

算額(その777)

2024年03月13日 | Julia

算額(その777)

川田保知:『算法極数小補解義』 文化15年戊寅正月
山口正義:やまぶき,第26号

https://yamabukiwasan.sakura.ne.jp/ymbk26.pdf

2 個の甲円が交差してできる区画に,乙円と丙円を入れる。丙円の直径が最大になるときの乙円の直径を求めよ。

甲円の半径と中心座標を r1, (r1, 0), (r1 + 2r2, 0)
乙円の半径と中心座標を r2, (r2, 0)
丙円の半径と中心座標を r3, (x3, y3)
とおき,r2 を変数のまま r3, x3, y3 を求める。
r3 は r2 の関数になるので,r3 を r2 で微分し,接線の傾きが 0 になるときの r2 を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, x3::positive, y3::positive
r1 = 35//2
eq1 = (r1 - x3)^2 + y3^2 - (r1 - r3)^2
eq2 = (x3 - r2)^2 + y3^2 - (r2 + r3)^2
eq3 = (r1 + 2r2 - x3)^2 + y3^2 - (r1 + r3)^2
res1 = solve([eq1, eq2, eq3], (r3, x3, y3));

丙円の半径は,乙円の半径の関数である。

r3 = res1[1][1] |> simplify
r3 |> println

   r2*(1225 - 4*r2^2)/(4*r2^2 + 1225)

res1[1][2] |> simplify |> println
res1[1][3] |> simplify |> println

   r2*(2*r2 + 35)^2/(4*r2^2 + 1225)
   70*r2*sqrt(1225 - 4*r2^2)/(4*r2^2 + 1225)

乙円の半径 r2 が 7.5〜10.0 の範囲内で丙円の半径 r3 が最大になる。

using Plots
pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(r3, xlims=(0, 17.5), xlabel="乙円の半径 r2", ylabel="丙円の半径 r3")

r3 を r2 で微分して,接線の傾きが 0 になるときの r2 を求めると 8.50269475574130 である。つまり,乙円の直径が 17.0053895114826 のとき丙円の直径が最も大きくなる。

f = diff(r3, r2)
res2 = solve(f, r2)[1]
res2 |> println
res2.evalf() |> println

   35*sqrt(-1/2 + sqrt(5)/4)
   8.50269475574130

もっとも大きいときの丙円の半径は 5.25495435501361(直径は 10.5099087100272 である)。

2r3(r2 => res2[1]).evalf() |> println

   10.5099087100272

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

   r1 = 17.5;  r2 = 8.50269;  r3 = 5.25495;  x3 = 15.1871;  y3 = 12.0246

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 35//2
   r2 = 35sqrt(√5 - 2)/2
   t = r2/(4r2^2 + 1225)
   (r3, x3, y3) = t .* (
       1225 - 4r2^2,
       (2*r2 + 35)^2,
       70sqrt(1225 - 4r2^2))
   @printf("乙円の直径が %g のとき,丙円の直径は最大値 %g になる。\n", 2r2, 2r3)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n", r1, r2, r3, x3, y3)
   plot()
   circle(r1, 0, r1)
   circle(r1 + 2r2, 0, r1)
   circle(r2, 0, r2, :blue)
   circle22(x3, y3, r3, :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(r1, 0, " 甲円:r1,(r1,0)", :red, :left, delta=-delta)
       point(r1 + 2r2, 0, " 甲円:r1,(r1+2r2,0)", :red, :left, delta=-delta)
       point(r2, 0, "乙円:r2,(r2,0)", :blue, :center, delta=-delta)
       point(x3, y3, "丙円:r3\n(x3,y3)", :green, :center, :bottom, delta=delta)
   end
end;

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

算額(その776)

2024年03月13日 | Julia

算額(その776)

埼玉県東松山市 岩殿観音(正法寺) 文政6年(1823)
山口正義:やまぶき,第30号

https://yamabukiwasan.sakura.ne.jp/ymbk30.pdf

二四 武州比企郡紫竹村 観世音堂 文政6年(1823)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.
キーワード:円6個,半円,正三角形

外円(半円)内に大円,小円,正三角形,および正三角形内に全名円,等円が入っている。正三角形の一辺の長さは外円の半径に等しい。等円の直径が与えられたとき,外円,大円,小円の直径を求めよ。

外円の半径と中心座標を R,(0, 0)
大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (x2, r2)
正三角形の一辺の長さは R
全名円の半径と中心座標を r3, (-R/2, r3)
等円の半径と中心座標を r4, (x4, r4)
とおく。
大円,小円の配置と正三角形及びその内部の円の配置は独立である。
既知の変数は r4 であるが,まずは R を基準として正三角形の内部の全名円,等円のパラメータを求める。
全名円の半径は r3 = 2R/√3 であることはすぐわかる。r4, x4 については以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, x1::positive, r2::positive, x2::positive,
     r3::positive, r4::positive, x4::negative
s3 = sqrt(Sym(3))
r3 = (R/2)/s3
eq1 = r4*(R/2) + r3*x4
eq2 = (R/2 + x4)^2 + (r3 - r4)^2 - (r3 + r4)^2;
res1 = solve([eq1, eq2], (r4, x4))

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (sqrt(3)*R/18, -R/6)
    (sqrt(3)*R/2, -3*R/2)

2 組の解が得られるが,最初のものが適解である。
等円の半径 r4 は r4 = R*√3/18 である。r4 が既知であるならば,R = 18r4/√3 である。
ついで,R = 18r4/√3 として,大円と小円のパラメータを求める。

R = 18r4/s3
eq3 = x1^2 + r1^2 - (R - r1)^2
eq4 = x2^2 + r2^2 - (R - r2)^2
eq5 = (x2 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq6 = dist(-R/2, s3*R/2, 0, 0, x1, r1) - r1^2;
res = solve([eq3, eq4, eq5, eq6], (r1, x1, r2, x2));

それぞれは若干長い式になるが,二重根号を外すなどの簡約化を行うと以下のようになる。

res[1][1] |> sympy.sqrtdenest |> simplify |> println
res[1][2] |> sympy.sqrtdenest |> simplify |> println
res[1][3] |> sympy.sqrtdenest |> simplify |> println
res[1][4] |> sympy.sqrtdenest |> simplify |> println

   18*r4*(2 - sqrt(3))
   6*r4*(-3 + 2*sqrt(3))
   18*r4*(-62*sqrt(2) - 41*sqrt(3) + 24*sqrt(6) + 150)/529
   6*r4*(-12*sqrt(2) - 2*sqrt(3) + 9 + 18*sqrt(6))/23

大円の直径は等円の直径の 18(2 - √3) 倍,小円の直径は等円の直径の 18(24√6 - 62√2 - 41√3 + 150)/529 倍である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   s3 = √3
   r4 = 1
   R = 18r4/s3
   x4 = -R/6
   r3 = R/2s3
   (r1, x1, r2, x2) = r4 .* (18(2 - √3), 6(2√3 - 3), 18(24√6 - 62√2 - 41√3 + 150)/529, 6(18√6 - 12√2 - 2√3 + 9)/23)
   @printf("R = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  r4 = %g;  x4 = %g\n", R, r1, x1, r2, x2, r3, r4, x4)
   plot([-R, 0, -R/2, -R], [0, 0, R√3/2, 0], color=:black, lw=0.5)
   circle(0, 0, R, beginangle=0, endangle=180)
   circle(x1, r1, r1, :blue)
   circle(x2, r2, r2, :orange)
   circle(-R/2, r3, r3, :green)
   circle(x4, r4, r4, :magenta)
   circle(-R - x4, r4, r4, :magenta)

   circle(-R/2, 2r3 + r4, r4, :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(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(x1, r1, "大円:r1,(x1,r1)", :blue, :center, delta=-delta)
       point(x2, r2, "小円:r2\n(x2,r2)", :black, :center, delta=-delta)
       point(-R/2, r3, "全名円:r3\n(-R/2,r3)", :green, :center, delta=-delta)
       point(x4, r4, "等円:r4,(x4,r4)", :magenta, :left, delta=-delta)
   end
end;

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

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

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