裏 RjpWiki

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

算額(その390)

2023年08月15日 | Julia

算額(その390)

埼玉県あきる野市 二宮神社 寛政6年(1794)
https://yamabukiwasan.sakura.ne.jp/ymbk17.pdf

鈎股弦(直角三角形)内に円と正方形が入っている。
股と長弦の和が 201寸6分,鈎と円径(円の直径)と方面(正方形の一辺)の和が 188 寸である。それぞれの長さはいくつか。

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

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive, 中鈎::positive,
     短弦::positive, 長弦::positive, 円径::positive, 方面::positive,
     x::positive, y::positive;
eq1 = 股 + 長弦 - 2016//10
eq2 = 鈎 + 円径 + 方面 - 188
eq3 = 中鈎^2 + 短弦^2 - 鈎^2
eq4 = 中鈎^2 + 長弦^2 - 股^2
eq5 = (鈎 - 方面)*(股 - 方面) - 方面^2
eq6 = 長弦 + 短弦 - 弦
eq7 = 鈎 + 股 - 弦 - 円径
eq8 = (鈎 + 股 + 弦)円径/2 - 鈎*股
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8],
   (鈎, 股, 弦, 中鈎, 短弦, 長弦, 円径, 方面))

   1-element Vector{NTuple{8, Sym}}:
    (84, 112, 140, 336/5, 252/5, 448/5, 56, 48)

股の長さは 112 寸である。

中鈎と弦の交点の座標 (x, y) を求める。

(鈎, 股, 弦, 中鈎, 短弦, 長弦, 円径, 方面) = (84, 112, 140, 336//5, 252//5, 448//5, 56, 48)
eq9 = x^2 + (鈎 - y)^2 - 短弦^2
eq10 = (股 - x)^2 + y^2 - 長弦^2
solve([eq9, eq10], (x, y))

   1-element Vector{Tuple{Sym, Sym}}:
    (1008/25, 1344/25)

   鈎 = 84;  股 = 112;  弦 = 140;  中鈎 = 67.2;  短弦 = 50.4;  長弦 = 89.6;  円径 = 56;  方面 = 48
   x = 40.32;  y = 53.76

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, 弦, 中鈎, 短弦, 長弦, 円径, 方面) = (84, 112, 140, 336/5, 252/5, 448/5, 56, 48)
   (x, y) = (1008/25, 1344/25)
   @printf("鈎 = %g;  股 = %g;  弦 = %g;  中鈎 = %g;  短弦 = %g;  長弦 = %g;  円径 = %g;  方面 = %g\n", 鈎, 股, 弦, 中鈎, 短弦, 長弦, 円径, 方面)
   @printf("x = %g;  y = %g\n", x, y)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(円径/2, 円径/2, 円径/2)
   rect(0, 0, 方面, 方面, :blue)
   segment(0, 0, x, y, :green)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, 鈎, "鈎 ", :black, :right, :vcenter)
       point(股, 0, "股", :black, :center, delta=-delta)
       point(方面, 0, "方面", :blue, :center, delta=-delta)
       point(円径/2, 円径/2, "(円径/2,円径/2)", :red, :center, delta=-delta)
       point(x, y, " (x,y)", :green, :left, :bottom)
       point(2x/3, 2y/3, "中鉤", :green, mark=false)
       point(x/2, (鈎 + y)/2, "   短弦 = (0,鈎)〜(x,y)", :black, mark=false)
       point((股 + x)/2, y/2, "   長弦 = (x,y)〜(股,0)", :black, mark=false)
       point(股/2, 鈎/2, "   弦 = 短弦 + 長弦", :black, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       plot!(xlims=(-8, 120), ylims=(-8, 90))
   else
       plot!(showaxis=false)
   end
end;

 

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

算額(その389)

2023年08月15日 | Julia

算額(その389)

埼玉県東松山市正代 世明寿寺 明治10年(1877)
山口正義 毛呂周辺の算額
https://yamabukiwasan.sakura.ne.jp/22moroshuuhen.pdf

鉤股弦内に,甲円,乙円,菱形を入れる。甲円,乙円の直径がそれぞれ 40 寸,28 寸のとき,菱形の一辺の長さはいかほどか。

x, y, z を以下のように定め,連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms 股::positive, x::positive, y::positive, z::positive;

(弦, r1, r2) = (180, 40//2, 32//2)
eq1 = 股^2 + x^2 - 弦^2
eq2 = y + z - (股 - z) - 2r2
eq3 = (x - y) + (股 - z) - (弦 - (股 - z)) - 2r1
eq4 = (x - y)z - (股 - z)y
res = solve([eq1, eq2, eq3, eq4], (股, x, y, z))

   2-element Vector{NTuple{4, Sym}}:
    (108, 144, 56, 42)
    (144, 108, 48, 64)

股 > x なので,2番目のものが適解である。
菱形の一辺は,股 - z = 80 である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (40, 32) .// 2
   (股, x, y, z) = (144, 108, 48, 64)
   plot([0, 股, 0, 0], [0, 0, x, 0], color=:black, lw=0.5)
   circle(r1, y + r1, r1)
   circle(r2, r2, r2, :blue)
   segment(0, y, 股 - z, y)
   segment(0, y, z, 0)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(r1, y + r1, "甲:r1\n(r1,y+r1)", :red, :center, delta=-delta)
       point(r2, r2, "乙:r2\n(r2,r2)", :blue, :center, delta=-delta)
       point(0, x, "x ", :black, :right, :vcenter)
       point(0, y, "y ", :black, :right, :vcenter)
       point(z, 0, " z", :black, :center, delta=-delta)
       point(股, 0, " 股", :black, :center, delta=-delta)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       plot!(xlims=(-7, 150), ylims=(-7, 110))
   else
       plot!(showaxis=false)
   end
end;

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

算額(その388)

2023年08月15日 | Julia

算額(その388)

埼玉県本庄市都島 正観寺 戸塚盛政の算額
山口正義 やまぶき 第40号
https://yamabukiwasan.sakura.ne.jp/ymbk40.pdf

甲,乙,丙,丁,戊の 5 個の正方形がある。甲,乙の正方形の一辺の 3 乗和が 189,丙,丁,戊の正方形の一辺の 3 乗和が 36 である。それぞれの正方形の一辺の長さは等差数列になっている。正方形の一辺の長さは,甲 > 乙 > 丙 > 丁 > 戊。
戊の正方形の一辺の長さと公差を求めよ。

この方程式系は,「9 次式を解くことになりかなり難しい」と書かれているが,SymPy では何の苦もなく解ける。
その9次式は,戊の正方形の一辺の長さを x としたとき,次のようなものである。
114604875x^9 - 9482437125x^6 + 200810066625x^3 - 191442234375 = 0
solve() で解くと,実数解は x = 1 だけで,それ以外は複素数解である。
ただ,この 9 次式は,定数 189,36 のときのものであり,その他の場合(例えば 700,500の場合には別の式になる)。

using SymPy

@syms x
eq = 114604875x^9 - 9482437125x^6 + 200810066625x^3 - 191442234375
res = solve(eq)

eq を因数分解すると以下のようになり,実数解は x = 1 のみであることが確認できる。

factor(eq) |> println

   1488375*(x - 1)*(x^2 + x + 1)*(77*x^6 - 6294*x^3 + 128625)

また,SymPy では,以下の 5 次元連立方程式を解くことで解を求めることができる。

using SymPy

@syms 甲::positive, 乙::positive, 丙::positive, 丁::positive, 戊::positive;

eq1 = 甲^3 + 乙^3 - 189
eq2 = 丙^3 + 丁^3 + 戊^3 - 36
eq3 = (甲 - 乙) - (乙 - 丙)
eq4 = (甲 - 乙) - (丙 - 丁)
eq5 = (甲 - 乙) - (丁 - 戊)
res = solve([eq1, eq2, eq3, eq4, eq5], (甲, 乙, 丙, 丁, 戊))
res |> println

   NTuple{5, Sym}[(5, 4, 3, 2, 1)]

一辺の長さは,甲 = 5,乙 = 4,丙 = 3,丁 = 2,戊 = 1 で,公差は 1 である。

この方法の場合も
eq1 = 甲^3 + 乙^3 - 700
eq2 = 丙^3 + 丁^3 + 戊^3 - 500
のような場合には解は一定の時間内では求まりそうにない。

 

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

算額(その387)

2023年08月15日 | Julia

算額(その387)

川田保知『淇澳集』
山口正義 やまぶき 第26号
https://yamabukiwasan.sakura.ne.jp/ymbk26.pdf

正方形の内外に大円,小円が配置されている。大円の直径が 2 寸のとき,小円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy

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

r1 = 2//2
x2 = y2 = r1 + r2/sqrt(Sym(2))
eq1 = (x2 - 2r1)^2 + y2^2 - (r1 + r2)^2
res = solve(eq1, r2)[1]
res |> println

   1/2

小円の直径は大円の直径の 1/2 である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (2, 1) .// 2
   @printf("r1 = %g;  r2 = %g\n", r1, r2)
   plot()
   rect(-r1, -r1, r1, r1, :black)
   circle(0, 0, r1)
   circle42(0, 2r1, r1)
   circle4(r1 + r2/√2, r1 + r2/√2, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, 2r1, " 2r1  大円", :red)
       point(r1 + r2/√2, r1 + r2/√2, " 小円:r2\n(r1+r2/√2,r1+r2/√2)", :blue, :left, :bottom, delta=delta)
       point(r1, 0, " r1", :red, :left, :bottom, delta=delta/2)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その386)

2023年08月15日 | Julia

算額(その386)

川田保知『算法極数小補解義』
山口正義 やまぶき 第26号
https://yamabukiwasan.sakura.ne.jp/ymbk26.pdf

鈎股弦に斜線を引き文斜と名付ける,文斜と股の交点の右側を武斜としたとき,鉤股弦の面積が最も大きくなるのは股がいくつのときか。

股 = x + 武斜 として,以下の二元連立方程式を解き,面積と鈎を求める。それぞれは x, 文斜,武斜 の関数になる。

include("julia-source.txt");

using SymPy

@syms 鈎::positive,  股::positive,
     弦::positive, 文斜::positive, 武斜::positive,
     x::positive, 面積::positive;

eq1 = 鈎^2 + x^2 - 文斜^2
eq2 = 鈎*(武斜 + x)/2 - 面積
res = solve([eq1, eq2], (面積, 鈎))

   1-element Vector{Tuple{Sym, Sym}}:
    (sqrt(-x + 文斜)*sqrt(x + 文斜)*(x + 武斜)/2, sqrt(-x + 文斜)*sqrt(x + 文斜))

面積は sqrt(-x + 文斜)*sqrt(x + 文斜)*(x + 武斜)/2 である。
文斜,武斜が特定の値 2, 7 をとるときに,面積が最大になる x と,そのときの面積を求める。

(文斜, 武斜) = (2, 7)
f = sqrt(2 - x)*sqrt(x + 2)*(x + 7)/2
f |> println

   sqrt(2 - x)*sqrt(x + 2)*(x + 7)/2

using Plots
pyplot(size=(400, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(f, xlims=(0, 2), xlabel="x", ylabel="面積")

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

   1/2
   0.500000000000000
   7.26184377413891

x が 0.5,すなわち 股 = 武斜 + 0.5 = 7.5 のとき,面積の最大値が 7.26184377413891 になる。

なお,このとき鈎は sqrt(2 - x)*sqrt(x + 2) = 1.9364916731037083 である。
1.9364916731037083 * 7.5 / 2 = 7.261843774138906 である。

   x = 0.5;  鈎 = 1.93649;  股 = 7.5;  面積 = 7.26184

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (文斜, 武斜) = (2, 7)
   x = 0.5
   股 = 武斜 + x
   (面積, 鈎) = (sqrt(-x + 文斜)*sqrt(x + 文斜)*(x + 武斜)/2, sqrt(-x + 文斜)*sqrt(x + 文斜))
   @printf("x = %g;  鈎 = %g;  股 = %g;  面積 = %g\n", x, 鈎, 股, 面積)
   plot()
   segment(0, 0, 0, 鈎, :red)
   segment(0, 0, x, 0, :blue)
   segment(x, 0, x + 武斜, 0, :magenta)
   segment(0, 鈎, x, 0, :green)
   segment(0, 鈎, 股, 0, :black)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, 鈎/2, "鈎 ", :red, :right, mark=false)
       point(x/2, 鈎/2, " 文斜", :green, mark=false)
       point(x + 武斜/2, 0, "武斜", :magenta, :left, :bottom, mark=false, delta=delta)
       point(x/2, 0, "x", :black, :left, :bottom, mark=false, delta=delta)
       dimension_line(0, -0.2, x + 武斜, -0.2, "\n股")
       plot!(xlims=(-0.5, 8), ylims=(-0.7, 2.3))
   else
       plot!(showaxis=false)
   end
end;

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

算額(その385)

2023年08月15日 | Julia

算額(その385)

川田保知『算法極数小補解義』
山口正義 やまぶき 第26号
https://yamabukiwasan.sakura.ne.jp/ymbk26.pdf


直径 35 寸の 2 個の甲円が交わっている。交わっていない部分に乙円 1 個と丙円 2 個を入れる。丙円の直径が最大になる時の乙円の直径を求めよ。

丙円の半径と中心座標を r3, (x3, y3)
乙円の半径と中心座標を r2, (r2, 0)
甲円の半径と中心座標を r1, (r1, 0) および (2r2 + r1, 0, r1)
とおき,以下の連立方程式をとき,r3, x3, y3 を求める。それぞれは r1 と r2 の関数である。

include("julia-source.txt");

using SymPy

@syms r1::positive,  r2::positive,
     r3::positive, x3::positive, y3::positive;

eq1 = (2r2 + r1 - x3)^2 + y3^2 - (r1 + r3)^2
eq2 = (x3 - r2)^2 + y3^2 - (r2 + r3)^2
eq3 = (x3 - r1)^2 + y3^2 - (r1 - r3)^2
res = solve([eq1, eq2, eq3], (r3, x3, y3))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (r2*(r1 - r2)*(r1 + r2)/(r1^2 + r2^2), r2*(r1 + r2)^2/(r1^2 + r2^2), 2*r1*r2*sqrt(r1 - r2)*sqrt(r1 + r2)/(r1^2 + r2^2))

r1 = 35//2
f = r2*(r1 - r2)*(r1 + r2)/(r1^2 + r2^2)
f |> println

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

pyplot(size=(400, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(f, xlims=(0, r1))

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

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

乙円の半径が 8.50269475574130 のときに,丙円の半径が最大値 5.25495435501361 になる。
すなわち,
乙円の直径が 17.0053895114826 のときに,丙円の直径が最大値 10.5099087100272 になる。

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

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 35//2
   r2 = res2
   (r3, x3, y3) = (r2*(r1 - r2)*(r1 + r2)/(r1^2 + r2^2), r2*(r1 + r2)^2/(r1^2 + r2^2), 2*r1*r2*sqrt(r1 - r2)*sqrt(r1 + r2)/(r1^2 + r2^2))
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n", r1, r2, r3, x3, y3)
   plot()
   circle(r1, 0, r1, :blue)
   circle(2r2 + r1, 0, r1, :blue)
   circle(r2, 0, r2)
   circle(x3, y3, r3, :green)
   circle(x3, -y3, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1, 0, " r1,甲円", :blue, delta=-delta)
       point(2r2 + r1, 0, " r2r2+1,甲円", :blue, delta=-delta)
       point(r2, 0, " r2,乙円", :red, delta=-delta)
       point(x3, y3, "丙円:r3\n(x3,y3)", :green, :center, delta=-delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

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

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