裏 RjpWiki

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

算額(その394)

2023年08月17日 | Julia

算額(その394)

東京都府中市 大國魂神社 明治18年(1885)
佐藤健一『多摩の算額』
山口正義『やまぶき』第8号
https://yamabukiwasan.sakura.ne.jp/ymbk8.pdf

もとは以下のものか。

東都愛宕山 天明8年戊申5月(1788)
東都大久保住 小林勝之助佐政

藤田貞資(1789):神壁算法
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

直線の上に,大円,小円,甲円,乙円,丙円が載っている。甲円,乙円の直径がそれぞれ 45 寸,40 寸のとき,丙円の直径を求めよ。

大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (0, r2)
甲円の半径と中心座標を r3, (x3, r3)
乙円の半径と中心座標を r4, (x4, r4)
丙円の半径と中心座標を r5, (x5, y5)
として,以下の連立方程式を立て,nlsove() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, x1::positive, r2::positive, r3::positive,
     x3::positive, r4::positive, x4::positive, r5::positive,
     x5::positive, y5::positive;

(r3, r4) = (45, 40) .// 2
eq1 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x1 - x5)^2 + (r1 - y5)^2 - (r1 + r5)^2
eq4 = x4^2 + (r2 - r4)^2 - (r2 + r4)^2
eq5 = x5^2 + (r2 - y5)^2 - (r2 + r5)^2
eq6 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
eq7 = (x3 - x5)^2 + (y5 - r3)^2 - (r3 + r5)^2
eq8 = (x5 - x4)^2 + (y5 - r4)^2 - (r4 + r5)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (r1, x1, r2, x3, x4, r5, x5, y5))

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)
   (r1, x1, r2, x3, x4, r5, x5, y5) = u
   return [
       x1^2 + (r1 - r2)^2 - (r1 + r2)^2,  # eq1
       (r1 - 45/2)^2 - (r1 + 45/2)^2 + (x1 - x3)^2,  # eq2
       -(r1 + r5)^2 + (r1 - y5)^2 + (x1 - x5)^2,  # eq3
       x4^2 + (r2 - 20)^2 - (r2 + 20)^2,  # eq4
       x5^2 - (r2 + r5)^2 + (r2 - y5)^2,  # eq5
       (x3 - x4)^2 - 1800,  # eq6
       -(r5 + 45/2)^2 + (x3 - x5)^2 + (y5 - 45/2)^2,  # eq7
       -(r5 + 20)^2 + (-x4 + x5)^2 + (y5 - 20)^2,  # eq8
   ]
end;

iniv = [big"120.0", 200, 70, 30, 20, 17, 17, 17]
res = nls(H, ini=iniv);
println(res);


   (BigFloat[180.0000000000000000315639921488776459575236334065858559417549014857970348855221, 254.5584412271571087953793012664368158326051549772724874493636018118273440778069, 89.99999999999999999414137791862049647293969079952184651651057986232141304475634, 127.2792206135785543912891769702404571933059905375609423589882744147464141850234, 84.85281374238570292723843621648770921579775164950411289786217536675179950468641, 17.99999999999999999772191395059268200439910616460700693283485445518118759929732, 101.8233764908628435131618139128282843444860986198208175883681931129179867092681, 53.99999999999999999717741648049545461928635214270576638615330713361282330565839], true)

r1 = 180;  x1 = 254.558;  r2 = 90;  x3 = 127.279;  x4 = 84.8528;  r5 = 18;  x5 = 101.823;  y5 = 54
丙円の直径 = 36

x1, x3, x4, x5 はそれぞれ 180√2, 90√2, 60√2, 72√2 である。

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x1, r2, x3, x4, r5, x5, y5) = res[1]
   (r3, r4) = (45, 40) .// 2
   @printf("r1 = %g;  x1 = %g;  r2 = %g;  x3 = %g;  x4 = %g;  r5 = %g;  x5 = %g;  y5 = %g\n", r1, x1, r2, x3, x4, r5, x5, y5)
   @printf("丙円の直径 = %g\n", 2r5)
   plot()
   circle(x1, r1, r1)
   circle(0, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   circle(x4, r4, r4, :magenta)
   circle(x5, y5, r5, :orange)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, r1, " 大:r1,(x1,r1)", :red, :left, :vcenter)
       point(0, r2, "小:r2,  (0,r2)", :blue, :center, :vcenter)
       point(x3, r3, "      甲:r3,(x3,r3)", :green, :left, :vcenter)
       point(x4, r4, "乙:r4,(x4,r4)   ", :magenta, :right, :bottom)
       point(x5, y5, "   丙:r5,(x5,y5) ", :orange, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その393)

2023年08月17日 | Julia

算額(その393)

埼玉県東秩父村御堂 浄蓮寺 明治34年(1901)
https://yamabukiwasan.sakura.ne.jp/ymbk7.pdf

代数の問題なので,3 問一括で記事にする。


第1問 大円内に小円がある。小円の直径の平方根と大円の直径の立方根の和が 4 寸,また小円の直径は大円の直径より 4 寸小さいとき,大円と小円の直径を求めよ。

小円,大円の半径を R1, R2 とし,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R1::positive, R2::positive;
eq1 = R1^(1//2) + R2^(1//3) - 4
eq2 = R2 - R1 - 4
solve([eq1, eq2], (R1, R2))

   1-element Vector{Tuple{Sym, Sym}}:
    (4, 8)

小円の直径は 4,大円の直径は 8 である。

using Plots

function draw()
    pyplot(size=(200, 200), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R1, R2) = (4, 8)
   plot(showaxis=false)
   circle(0, 0, R2/2)
   circle(0, -2, R1/2, :blue)

end;


第2問 鈎股弦において,弦が 997 寸,面積が 172050 平方寸のとき鈎はいくつか。

鈎,股を 「鉤」,「股」として,以下の連立方程式を解く。

@syms 鈎::positive, 股::positive, 弦::positive
(弦, 面積) = (997, 172050)
eq1 = 鈎^2 + 股^2 - 弦^2
eq2 = 鈎*股/2 - 面積
solve([eq1, eq2], (鈎, 股))

   2-element Vector{Tuple{Sym, Sym}}:
    (372, 925)
    (925, 372)

鈎 < 股 なので,鈎 = 372 寸,股 = 925 寸である。


第3問 方錐(四角錐)の体積が 32 立方寸,底面の正方形の方面(辺)の平方根と高さの和が 8(無名数)のとき,方面と高さはいくつか。

底面の正方形の方面を x,高さを h として,以下の連立方程式を解く。

@syms x::positive, h::positive
eq1 = x^2*h/3 - 32
eq2 = sqrt(x) + h - 8
res = solve([eq1, eq2], (x, h))

   2-element Vector{Tuple{Sym, Sym}}:
    (4, 6)
    (15 + 2*sqrt(41/(9*sqrt(2089)/16 + 1541/16)^(1/3) + 2*(9*sqrt(2089)/16 + 1541/16)^(1/3) + 265/4) + 2*sqrt(-2*(9*sqrt(2089)/16 + 1541/16)^(1/3) - 41/(9*sqrt(2089)/16 + 1541/16)^(1/3) + 4203/(4*sqrt(41/(9*sqrt(2089)/16 + 1541/16)^(1/3) + 2*(9*sqrt(2089)/16 + 1541/16)^(1/3) + 265/4)) + 265/2), 96/(15 + 2*sqrt(41/(9*sqrt(2089)/16 + 1541/16)^(1/3) + 2*(9*sqrt(2089)/16 + 1541/16)^(1/3) + 265/4) + 2*sqrt(-2*(9*sqrt(2089)/16 + 1541/16)^(1/3) - 41/(9*sqrt(2089)/16 + 1541/16)^(1/3) + 4203/(4*sqrt(41/(9*sqrt(2089)/16 + 1541/16)^(1/3) + 2*(9*sqrt(2089)/16 + 1541/16)^(1/3) + 265/4)) + 265/2))^2)

res[1][1].evalf(), res[1][2].evalf()

   (4.00000000000000, 6.00000000000000)

res[2][1].evalf(), res[2][2].evalf()

   (63.6210823302620, 0.0237175118817425)

算額の答えは,「方面 4 寸,高さ 6 寸」ただ一つであるが,答えはもう一組ある。
方面 63.6210823302620 寸, 高さ 0.0237175118817425 寸」という,すごく扁平な四角錐も条件を満たす。

このような,整数解と実数解が併存する問題はときどきある。和算のやり方では実数解の存在に気づかないのかもしれない。

(x, h) = res[2]
x^2 * h / 3 |> println
(sqrt(x) + h).evalf() |> println

   32
   8.00000000000000

eq2 から h = 8 - sqrt(x) を,eq1 に代入する

@syms x
eq = x^2*(8 - sqrt(x)) - 96
res = solve(eq, x)
res[1].evalf(), res[2].evalf()

   (4.00000000000000, 63.6210823302620)

eq を図に描くと以下のようになる

include("julia-source.txt");
using Plots
pyplot(size=(400, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(eq, xlims=(0, 65), xlabel="x", ylabel="体積")
point(res[1], 0, res[1].evalf())
point(res[2], 0, res[2].evalf(), :green, :right)
hline!([0])

 

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

算額(その392)

2023年08月16日 | Julia

算額(その392)

埼玉県鴻巣市 新井稲荷神社 明治25年(1892)
https://yamabukiwasan.sakura.ne.jp/ymbk351.pdf

外円の中に弦を入れ,その上下に大円,小円を入れる。弦の長さが 4 寸,小円の直径が 1 寸のとき,大円の直径はいかほどか。

外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (x1, r1 - r0 + 2r2)
小円の半径と中心座標を r2, (x2, 3r2 - r0)
弦の長さを「弦」
とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, x1::negative, r2::positive, x2::positive, 弦::positive;

eq1 = (x2 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x2^2 + (3r2 - r0)^2 - (r0 - r2)^2
eq3 = x1^2 + (r1 -r0 +2r2)^2 - (r0 - r1)^2
eq4 = r0^2 - (2r2 - r0)^2 - (弦//2)^2

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

   2-element Vector{NTuple{4, Sym}}:
    (r2 + 弦^2/(16*r2), (-16*r2^2 + 2*弦^2 - sqrt(-4*r2 + 弦)*sqrt(4*r2 + 弦)*sqrt(-16*r2^2 + 弦^2) - sqrt(-256*r2^4 + 弦^4))/(32*r2), sqrt(-16*r2^2 + 弦^2)/4 + sqrt(16*r2^2 + 弦^2)/4, sqrt(-16*r2^2 + 弦^2)/2)
    (r2 + 弦^2/(16*r2), (-16*r2^2 + 2*弦^2 - sqrt(-4*r2 + 弦)*sqrt(4*r2 + 弦)*sqrt(-16*r2^2 + 弦^2) + sqrt(-256*r2^4 + 弦^4))/(32*r2), sqrt(-16*r2^2 + 弦^2)/4 - sqrt(16*r2^2 + 弦^2)/4, sqrt(-16*r2^2 + 弦^2)/2)

2 組目の解が適解である。

大円の直径は 3.93649167310371 寸である。

2*res[2][2](r2 => 1//2, 弦 => 4).evalf() |> println

    3.93649167310371

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, 弦) = (1//2, 4)
   (r0, r1, x1, x2) = (r2 + 弦^2/(16*r2), (-16*r2^2 + 2*弦^2 - sqrt(-4*r2 + 弦)*sqrt(4*r2 + 弦)*sqrt(-16*r2^2 + 弦^2) + sqrt(-256*r2^4 + 弦^4))/(32*r2), sqrt(-16*r2^2 + 弦^2)/4 - sqrt(16*r2^2 + 弦^2)/4, sqrt(-16*r2^2 + 弦^2)/2)
   b = sqrt(r0^2 - (2r2 - r0)^2)
   @printf("r0 = %g;  r1 = %g;  x1 = %g;  x2= %g\n大円の直径 = %g;  弦 = %g;  小円の直径 = %g\n", r0, r1, x1, x2, 2r1, 2b, 2r2)
   plot()
   circle(0, 0, r0, :black)
   circle(x1, r1 - r0 +2r2, r1)
   circle(x2, 3r2 - r0, r2, :blue)
   circle(0, r2 - r0, r2, :blue)
   segment(-b, 2r2 - r0, b, 2r2 - r0)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(x1, r1 - r0 + 2r2, " 大円:r1,(x1,r1-r0+2r2) ", :red, :right, :vcenter)
       point(x2, 3r2 - r0, " 小円:r1,(x1,3r2−r0) ", :blue, :right, :vcenter)
       point(0, r2 - r0, " r2-r0", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

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

算額(その391)

2023年08月16日 | Julia

算額(その391)

埼玉県鴻巣市 新井稲荷神社 明治25年(1892)https://yamabukiwasan.sakura.ne.jp/ymbk351.pdf

大円,中円,小円,乙円,丙円が図のように配置されている。小円,乙円,丙円の直径がそれぞれ 5 寸,4 寸,3 寸のとき,大円の直径はいかほどか。

大円の半径と中心座標を r1, (0, r3 + 2r5 + 2r4 + r1)
中円の半径と中心座標を r2, (x2, y2)
小円の半径と中心座標を r3, (0, 0)
乙円の半径と中心座標を r4, (0, r3 + 2r5 + r4) 
丙円の半径と中心座標を r5, (0, r3 + r5)
とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::positive, y2::positive,
     r3::positive, r4::positive, r5::positive;

eq1 = x2^2 + (r3 + 2r5 + 2r4 + r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + y2^2 - (r2 + r3)^2
eq3 = x2^2 + (r3 + 2r5 + r4 - y2)^2 - (r2 + r4)^2
eq4 = x2^2 + (r3 + r5 - y2)^2 - (r2 + r5)^2

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

   2-element Vector{NTuple{4, Sym}}:
    (-r3*r4^2/(r3*r4 - r3*r5 - r5^2), r5*(r3 + r5)*(r4 + r5)/(r3*r4 - r5^2), -2*sqrt(r3)*sqrt(r4)*r5*sqrt(r3 + r5)*sqrt(r4 + r5)/(r3*r4 - r5^2), (r3^2*r4 + r3*r4*r5 - r4*r5^2 - r5^3)/(r3*r4 - r5^2))
    (-r3*r4^2/(r3*r4 - r3*r5 - r5^2), r5*(r3 + r5)*(r4 + r5)/(r3*r4 - r5^2), 2*sqrt(r3)*sqrt(r4)*r5*sqrt(r3 + r5)*sqrt(r4 + r5)/(r3*r4 - r5^2), (r3^2*r4 + r3*r4*r5 - r4*r5^2 - r5^3)/(r3*r4 - r5^2))

(r3, r4, r5) = (5, 4, 3) .// 2
(-r3*r4^2/(r3*r4 - r3*r5 - r5^2), r5*(r3 + r5)*(r4 + r5)/(r3*r4 - r5^2), -2*sqrt(r3)*sqrt(r4)*r5*sqrt(r3 + r5)*sqrt(r4 + r5)/(r3*r4 - r5^2), (r3^2*r4 + r3*r4*r5 - r4*r5^2 - r5^3)/(r3*r4 - r5^2))

   (10//1, 84//11, -9.127200289462643, 97//22)

(-r3*r4^2/(r3*r4 - r3*r5 - r5^2), r5*(r3 + r5)*(r4 + r5)/(r3*r4 - r5^2), 2*sqrt(r3)*sqrt(r4)*r5*sqrt(r3 + r5)*sqrt(r4 + r5)/(r3*r4 - r5^2), (r3^2*r4 + r3*r4*r5 - r4*r5^2 - r5^3)/(r3*r4 - r5^2))

   (10//1, 84//11, 9.127200289462643, 97//22)

解は 2 組得られるが,両者の違いは,中円の中心座標が右のものか左のものかの違いだけである。

大円の直径は 2(-r3*r4^2/(r3*r4 - r3*r5 - r5^2)) = 20 寸である。
術では 「3920 を 196 で割ればよい」とあるが,なぜそのような数値が出てくるか不明である。

   r1 = 10;  r2 = 7.63636;  x2 = 9.1272;  y2 = 4.40909
   大円の直径 = 20;  中円の直径 = 15.2727

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, r4, r5) = (5, 4, 3) .// 2
   (r1, r2, x2, y2) = (-r3*r4^2/(r3*r4 - r3*r5 - r5^2), r5*(r3 + r5)*(r4 + r5)/(r3*r4 - r5^2), 2*sqrt(r3)*sqrt(r4)*r5*sqrt(r3 + r5)*sqrt(r4 + r5)/(r3*r4 - r5^2), (r3^2*r4 + r3*r4*r5 - r4*r5^2 - r5^3)/(r3*r4 - r5^2))
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", r1, r2, x2, y2)
   @printf("大円の直径 = %g;  中円の直径 = %g\n", 2r1, 2r2)
   plot()
   circle(0, r3 + 2r5 + 2r4 + r1, r1)
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   circle(0, 0, r3, :orange)
   circle(0, r3 + 2r5 + r4, r4, :magenta)
   circle(0, r3 + r5, r5, :green)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, r3 + 2r5 + 2r4 + r1, " r3+2r5+2r4+r1\n\n      大円:r1", :red)
       point(x2, y2, "中円:r2,(x2,y2)", :blue, :center, delta=-delta)
       point(0, r3 + 2r5 + r4, "乙円:r4  r3+2r5+r4 ", :magenta, :right)
       point(0, r3 + r5, "丙円:r5  r3+r5 ", :green, :right)
       point(0, 0, "小円:r5 ", :orange, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その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でシェアする

算額(その384)

2023年08月12日 | Julia

算額(その384)

埼玉県北埼玉郡根古屋 氷川神社
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html

外円の中に甲円 2 個,乙円 6 個を入れる。乙円の直径が 1 寸のとき,甲円の直径はいくらか。

外円の半径と中心座標を r0,(0, 0)
甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (x1, 0) および (x2, y2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms x::positive,  y::positive, x1::positive, x2::positive,
   y2::positive, r1::positive, r2::positive, r0::positive;

r0 = 2r1
eq1 = x1^2 + r1^2 - (r1 + r2)^2
eq2 = x2^2 + (r1 - y2)^2 - (r1 + r2)^2
eq3 = x2^2 + y2^2 - (r0 - r2)^2
eq4 = (x2 - x1)^2 + y2^2 - 4r2^2
res = solve([eq1, eq2, eq3, eq4], (r1, x1, x2, y2))

   1-element Vector{NTuple{4, Sym}}:
    (r2*(1 + sqrt(73))/4, 4*sqrt(2)*r2/sqrt(-3 + sqrt(73)), r2*sqrt(-6 + 2*sqrt(73)), r2*(-5/2 + sqrt(73)/2))

  r2 = 1 のとき r1 = 2.386;  x1 = 2.4025;  x2 = 3.32987;  y2 = 1.772;  r0 = 4.772

甲円の直径は 乙円の直径の (1 + √73)/4 = 2.3860009363293826 倍である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1
   (r1, x1, x2, y2) = (r2*(1 + sqrt(73))/4, 4*sqrt(2)*r2/sqrt(-3 + sqrt(73)), r2*sqrt(-6 + 2*sqrt(73)), r2*(-5/2 + sqrt(73)/2))
   r0 = 2r1
   @printf("r1 = %g;  x1 = %g;  x2 = %g;  y2 = %g;  r0 = %g\n", r1, x1, x2, y2, r0)
   plot()
   circle(0, 0, r0, :black)
   circle(0, r1, r1)
   circle(0, -r1, r1)
   circle(x1, 0, r2, :blue)
   circle(-x1, 0, r2, :blue)
   circle4(x2, y2, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, 2r1, " r0", delta=-delta/2)
       point(0, r1, " r1  甲円")
       point(x1, 0, " x1\n 乙円", :green, :center, delta=-delta/2)
       point(x2, y2, " (x2,y2)\n 乙円", :green, :center, 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でシェアする

算額(その383)

2023年08月12日 | Julia

算額(その383)

番外編 2019 青森公立大学 前期 経営経済学部 第 3 問
https://mathexamtest.web.fc2.com/2019/201911051/2019110510100mj.html#top-0105

BC を x 軸上に置き,点 B を原点とする。
A の座標 (ax, ay)
B の座標 (0, 0)
C の座標 (cx, 0)
D の座標 (dx, 0)
E の座標 (ex, ey)
外接円の半径と中心座標 R, (ox, oy)
AD*cos(∠BAC) を ADcosθ とおく(使わない)。
AB = 7, CA = 13, R = 13√3/3
以下の 9元連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms ax::positive, ay::positive, cx::positive,
   dx::positive, ex::positive, ey::negative,
   ox::positive, oy::negative, R::positive,
   ADcosθ::poisitive, AD::positive,
   AB::positive, CA::positive;

(AB, CA, R) = (7, 13, 13*sqrt(Sym(3))/3)
AD = sqrt((ax - dx)^2 + ay^2)
eq1 = ax^2 + ay^2 - 7^2
eq2 = (cx - ax)^2 + ay^2 - 13^2
eq3 = (ax - ox)^2 + (ay - oy)^2 - R^2
eq4 = ox^2 + oy^2 - R^2
eq5 = (cx- ox)^2 + oy^2 - R^2
eq6 = (ex - ox)^2 + (ey - oy)^2 - R^2
eq7 = dx*13 - AD*sqrt(ex^2 + ey^2)
eq8 = AB^2 + AD^2 - 2*AB*ADcosθ - dx^2
eq9 = CA^2 + AD^2 - 2*CA*ADcosθ - (cx - dx)^2

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (ax, ay, cx, dx, ex, ey, ox, oy, ADcosθ))

   2-element Vector{NTuple{9, Sym}}:
    (7/2, 7*sqrt(3)/2, 15, 21/4, 15/2, -9*sqrt(3)/2, 15/2, -sqrt(3)/6, 35/8)
    (7/2, 7*sqrt(3)/2, 15, 21/4, 15/2, -9*sqrt(3)/2, 15/2, -sqrt(3)/6, 35/8)

2 組求まったが,全く同じだ。

   ax = 3.5;  ay = 6.06218;  cx = 15;  dx = 5.25;
   ex = 7.5;  ey = -7.79423;  ox = 7.5;  oy = -0.288675;  AD*cos(θ) = 4.375
   辺 BC の長さ = 15
   線分 AD の長さ = 7*sqrt(13)/4
   線分 BE の長さ = 3*sqrt(13)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (ax, ay, cx, dx, ex, ey, ox, oy, ADcosθ) = res[2]
   R = 13*sqrt(Sym(3))/3
   @printf("ax = %g;  ay = %g;  cx = %g;  dx = %g;\nex = %g;  ey = %g;  ox = %g;  oy = %g;  AD*cos(θ) = %g\n",
           ax, ay, cx, dx, ex, ey, ox, oy, ADcosθ)
   println("辺 BC の長さ = $cx")
   println("線分 AD の長さ = $(sqrt((ax - dx)^2 + ay^2))")
   println("線分 BE の長さ = $(sqrt(ex^2 + ey^2))")
   plot([0, ax, ex, 0, cx, ax], [0, ay, ey, 0, 0, ay], color=:black, lw=0.5)
   circle(ox, oy, R)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(ax, ay, "A(ax,ay) ", :green, :right, :bottom)
       point(0, 0, "B(0,0) ", :green, :right, :bottom)
       point(cx, 0, " C(cx,0)", :green, :left, :bottom)
       point(dx, 0, " D(dx,0)", :green, :left, :bottom)
       point(ex, ey, " E(ex,ey)", :green, :left, :bottom, delta=1.5delta)
       point(ox, oy, " O(ox,oy)", :red)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       plot!(xlims=(-2, 16))
   else
       plot!(showaxis=false)
   end
end;

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

算額(その382)

2023年08月11日 | Julia

算額(その382)

東京都渋谷区 金王八幡宮(金王八幡神社) 嘉永3年(1850)
https://www.naruto-u.ac.jp/journal/info-edu/j05006.pdf

二十八宿のうち15宿の名前をつけた円がある。角と亢の円周の和が16寸,心,尾,箕の円周の和が30寸,虚,危,室,壁,奎 の周の和が63寸であるとき,角の円周を求めよ。


答は 7 寸 7 分 6 厘 3 毛 3 糸 2 忽 1 微と残りは微少数(7.763321...)である。

include("julia-source.txt");

using Plots

function draw()
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   f(n) = -69//10960*n^2 + 1079//2192*n + 997//137
   names = ["角", "亢", "氐", "房", "心", "尾", "箕", "斗", "牛", "女", "虚", "危", "室", "壁", "奎"]
   plot()
   for (i, name) in enumerate(names)
       point(i, 1, name, :blue, :center, mark=false)
       circle(i, 0, f(i)/30)
   end
   plot!(showaxis=false)
end;

金王八幡宮といえば,冲方丁の「天地明察」にも出てくる。その中に,この問題を参考にしたと思われる算額問題がある。これについては,末尾に記す。

15宿の名前を順番に書いておく。角が一番小さく,奎が一番大きい。
角,亢,氐,房,心,尾,箕,斗,牛,女,虚,危,室,壁,奎

天文(学)に,しかも昔の天文(学)に興味がなければ,読み方も分かりづらいので,この算額の題意を,現代数学風に書き下す。

数列 a[n], n = 1, 2, ..., 15 がある。
a[i] < a[i+1], i = 1, 2, ..., 14 である。
a[1] + a[2] = 16
a[5] + a[6] + a[7] = 30
a[11] + a[12] + a[13] + a[14] + a[15] = 63
がわかっているとき,a[1] を求めよ。

等差数列ではなさそうだが,次に簡単な形式の数列といえば a[n] = p×n^2 + q×n + r であろうか。

第1の条件 a[1] + a[2] = 16 は
(p×1^2 + q×1 + r) + (p×2^2 + q×2 + r)
= 5p + 3q + 2r
のようになる。
2, 3 番目の条件における p, q, r の係数を求める。

using SymPy

function coef(s, e)
   coef2 = coef1 = coef0 = 0
   for n in s:e
       coef2 += n^2
       coef1 += n
       coef0 += 1
   end
   return (coef2, coef1, coef0)
end

coef(1, 2) |> println
coef(5, 7) |> println
coef(11, 15) |> println

   (5, 3, 2)
   (110, 18, 3)
   (855, 65, 5)

以下の3元連立方程式を解けばよい。

@syms p, q, r
eq1 = 5p + 3q + 2r - 16
eq2 = 110p + 18q + 3r - 30
eq3 = 855p + 65q + 5r - 63
res = solve([eq1, eq2, eq3], (p, q, r))

   Dict{Any, Any} with 3 entries:
     p => -69/10960
     q => 1079/2192
     r => 997/137

分数形式を小数点付きで表す。

res[p].evalf(), res[q].evalf(), res[r].evalf()

   (-0.00629562043795620, 0.492244525547445, 7.27737226277372)

a[1] は 1^2p + 1q + r なので,7.76332116788321 が答えである。

(res[p] + res[q] + res[r]).evalf() |> println

   7.76332116788321

行列演算で行うならば,以下のようにする。

A = Sym[5 3 2; 110 18 3; 855 65 5]
b = [16, 30, 63];

pqr = A \ [16.0, 30, 63]
pqr |> println

   Sym[-0.00629562043795620, 0.492244525547445, 7.27737226277372]

pqr' * [1;1;1] |> println

   7.76332116788321

演算子 \ は Julia 特有のものなので,一般言語的には以下のようにする。

inv(A) * b |> println

   Sym[-69/10960, 1079/2192, 997/137]

---------------------------------

さて,冲方丁の天地明察に出てくる問題である。
冲方 丁. 天地明察(特別合本版) (Japanese Edition) (p.237). Kindle 版. 

『今、図の如く、大小の十四宿の星の名を持つ円が並んでいる。角星と亢星の周の長さを足すと九寸である。また房星と心星と尾星の周の長さを足すと十八寸である。さらに女星、虚星、危星、室星、壁星の五つの星の周の長さを足すと四十五寸である。角星の周の長さは何寸であるか問う』

変更点は以下の通り。

- 15宿から14宿に縮小された
- 条件数は変わらず3個であるが,係数行列と定数ベクトルは異なる
- 2 番目の条件式で,金王八幡宮のものは「心,尾,箕の円周の和」であるが天地明察のものでは「房星と心星と尾星の周の長さを足す」となっている
- p*n^2+q*n+r だが,(結果的に)p = 0 である。つまり,簡単な階差数列である

coef(1, 2) |> println
coef(4, 6) |> println
coef(10, 14) |> println


   (5, 3, 2)
   (77, 15, 3)
   (730, 60, 5)

以下の3元連立方程式を解けばよい。

@syms p, q, r
eq1 = 5p + 3q + 2r - 9
eq2 = 77p + 15q + 3r - 18
eq3 = 730p + 60q + 5r - 45

$730 p + 60 q + 5 r - 45$

res = solve([eq1, eq2, eq3])

   Dict{Any, Any} with 3 entries:
     p => 0
     q => 3/7
     r => 27/7

分数形式を小数点付きで表す。なお,p = 0 なので,一般項は
a[n] = 27//7 + 3/7*n
つまり,階差数列である。これを村瀬は「……招差術か」と言っている。

res[p].evalf(), res[q].evalf(), res[r].evalf()

   (0, 0.428571428571429, 3.85714285714286)

a[1] は 1^2p + 1q + r なので,7.76332116788321 が答えである。

res[r] + res[q]

   30/7

1^2*res[p] + 1*res[q] + 1*res[r] |> println

   30/7

因縁の "30/7" が出てきた。

一般項 a(n) を求める関数を定義しておこう。

func_a(n) = (3n + 27)//7;

for i in 1:14
   println("func_a($i) = $(func_a(i))")
end

   func_a(1) = 30//7
   func_a(2) = 33//7
   func_a(3) = 36//7
   func_a(4) = 39//7
   func_a(5) = 6//1
   func_a(6) = 45//7
   func_a(7) = 48//7
   func_a(8) = 51//7
   func_a(9) = 54//7
   func_a(10) = 57//7
   func_a(11) = 60//7
   func_a(12) = 9//1
   func_a(13) = 66//7
   func_a(14) = 69//7

func_a(1) + func_a(2)

   9//1

func_a(4) + func_a(5) + func_a(6)

   18//1

func_a(10) + func_a(11) + func_a(12) + func_a(13) + func_a(14)

   45//1

角星の周の長さは 30/7

func_a(1) |> println

   30//7

行列演算で行うならば,以下のようにする。

A = Sym[5 3 2; 77 15 3; 730 60 5]
b = [9, 18, 45];

A, b

   (Sym[5 3 2; 77 15 3; 730 60 5], [9, 18, 45])

pqr = A \ [9, 18, 45]
pqr |> println

   Sym[0, 3/7, 27/7]

pqr' * [1;1;1] |> println

   30/7

演算子 \ は Julia 特有のものなので,一般言語的には以下のようにする。

inv(A) * b |> println

   Sym[0, 3/7, 27/7]

---

なお,この数列の問題は,初版本(単行本)では解答不能ではないものの,なぜこのような答えになる問題を出したのかという状態であった。

236ページ
『今、図の如く、大小の十五宿の星の名を持つ円が並んでいる。角星と亢星の周の長さを足すと十寸である。また心星と尾星と箕星の周の長さを足すと二十七寸五分である。さらに虚星、危星、室星、壁星、奎星の五つの星の周の長さを足すと四十寸である。角星の周の長さは何寸であるか問う』

角星 + 亢星 = 10寸
心星 + 尾星 + 箕星 = 27寸5分
虚星 + 危星 + 室星 + 壁星 + 奎星 = 40寸

using SymPy

A = Sym[5 3 2; 110 18 3; 855 65 5]
b = [10, 27.5, 40];
pqr = A\b
pqr |> println

   Sym[-0.0942062043795621, 1.64119525547445, 2.77372262773723]

g(n) = -0.0942062043795621*n^2 + 1.64119525547445*n + 2.77372262773723

using SymPy
solve(diff(g))

   8.71065375302662

g(7), g(8), g(9), g(10)

   (9.645985401459837, 9.874087591240855, 9.91377737226275, 9.765054744525518)

9 番目が最大でそれ以降では減少する。そういう意味では,この問題は「無術」なのであろう。

角星の周の長さは 4.32071167883212 ということになる。

pqr' * [1;1;1]

   4.32071167883212

しかし,この数値も特段の意味を持つものではないようだ。

ということで,改訂版(電子版)では著者の意図したように問題(の数値)が変更された

 

 

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

算額(その381)

2023年08月10日 | Julia

算額(その381)

静岡県伊豆市 江川邸 享和2年(1802)
http://www.wasan.jp/sizuoka/egawa.html

長野県長野市若穂 清水寺観音堂 寛政6年(1794)
中村信弥「改訂増補 長野県の算額」(45)
http://www.wasan.jp/zoho/zoho.html

正三角形内に 3 本の斜線を引き,分割された領域に 4 個の等円を内接させる。正三角形の一辺の長さがが与えられるとき,等円の直径はいかほどか。

正三角形の一辺の長さを a,等円の半径を r とする。

以下の連立方程式を解けばよいが,solve() では無理のようなので,nlsove() で数値解を求める。解析解は後半に記載。

include("julia-source.txt");

using SymPy

@syms a::positive, r::positive, x::positive, y::positive,
     x1::positive,y1::positive;

eq1 = distance(0, √Sym(3)a/2, a/2, 0, x, y) - r^2
eq2 = distance(0, √Sym(3)a/2, x1, y1, x, y) - r^2
eq3 = distance(0, √Sym(3)a/2, x1, y1, 0, a/2√Sym(3)) - r^2
eq4 = distance(x1, y1, a/2, 0, x, y) - r^2
eq5 = distance(x1, y1, a/2, 0, 0, a/2√Sym(3)) - r^2;

# res = solve([eq1, eq2, eq3, eq4, eq5])

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)
   (r, x, y, x1, y1) = u
   return [
-r^2 + (-3*a/8 + 3*x/4 + sqrt(3)*y/4)^2 + (-sqrt(3)*a/8 + sqrt(3)*x/4 + y/4)^2,  # eq1
-r^2 + (x - x1*(3*a^2 - 2*sqrt(3)*a*y - 2*sqrt(3)*a*y1 + 4*x*x1 + 4*y*y1)/(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2))^2 + (y - (3*a^2*y - 2*sqrt(3)*a*x*x1 + 2*sqrt(3)*a*x1^2 - 4*sqrt(3)*a*y*y1 + 4*x*x1*y1 + 4*y*y1^2)/(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2))^2,  # eq2
4*a^2*x1^2*(3*a - 2*sqrt(3)*y1)^2/(9*(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2)^2) - r^2 + (sqrt(3)*a/6 - sqrt(3)*a*(3*a^2 - 4*sqrt(3)*a*y1 + 12*x1^2 + 4*y1^2)/(6*(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2)))^2,  # eq3
-r^2 + (x - (a^2*x - 4*a*x*x1 - 2*a*y*y1 + 2*a*y1^2 + 4*x*x1^2 + 4*x1*y*y1)/(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2))^2 + (y - y1*(a^2 - 2*a*x - 2*a*x1 + 4*x*x1 + 4*y*y1)/(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2))^2,  # eq4
a^2*y1^2*(-sqrt(3)*a + 2*sqrt(3)*x1 + 6*y1)^2/(9*(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2)^2) - r^2 + (-a*y1*(3*a - 6*x1 + 2*sqrt(3)*y1)/(3*(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2)) + sqrt(3)*a/6)^2,  # eq5
   ]
end;

a = 10
iniv = [big"0.09", 0.27, 0.19, 0.17, 0.09] .* a
res = nls(H, ini=iniv);
println(res);
  (BigFloat[1.142125629369641587973968514127294180492236375359075288930728366920198239037438, 2.50000000000000007079306505041260465903195757502042722173815452658909699459355, 2.045875760182909850411741235498302045060062420025658180702662505806048892327303, 1.510890190652600068372820337416072773295486747795635718444765203143264858414334, 1.173562901893666266337957465962928304869915646101974796747607212272768807266053], true)

   r = 1.14213;  a = 10;  x = 2.5;  y = 2.04588;  x1 = 1.51089;  y1 = 1.17356
   小円の直径 = 2r = 2.28425

解析解

正三角形の面積を 3 個の相似な三角形と中央の小さな正三角形の面積に分割する。
中央の正三角形の面積は 3√3r^2 である(注)。
下側の三角形に注目する。鈍角は 120°であるので,3 辺の和は 2a + 2*r/√3 で,面積は (2a + 2*r/√3)r/2である。
よって,
√3a^2/4 = 3(2a + 2*r/√3)r/2 + 3√3r^2
を r について解けばよい。

注: 正三角形の内接円の半径が r のとき,正三角形の面積は 3√3r^2 である。

@syms a::positive, r::positive
eq = 3(2a + 2*r/√Sym(3))r/2 + 3√Sym(3)r^2 - √Sym(3)a^2/4 |> simplify
eq |> println

   -sqrt(3)*a^2/4 + 3*a*r + 4*sqrt(3)*r^2

solve(eq, r)[1] |> println

   a*(-sqrt(3) + sqrt(7))/8

直径は (√7 - √3)/4 * a である。
a = 10 のとき 2.2842512587392836 である。

(√7 - √3)/4 * 10

   2.2842512587392836

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 10
   (r, x, y, x1, y1) = res[1]
   @printf("r = %g;  a = %g;  x = %g;  y = %g;  x1 = %g;  y1 = %g\n", r, a, x, y, x1, y1)
   @printf("小円の直径 = 2r = %g\n", 2r)
   plot([a/2, 0, -a/2, a/2], [0, √3a/2, 0, 0], color=:black, linewidth=0.25)
   circle(0, a/2√3, r)
   for deg = 0:120:240
       A = [cosd(deg) -sind(deg); sind(deg) cosd(deg)]
       (x2, y2, p1, p2, p3, p4) = A * [x 0 x1; y-a/2√3 √3a/2-a/2√3 y1-a/2√3]
       segment(p1, p2 + a/2√3, p3, p4 + a/2√3, :blue)
       circle(x2, y2 + a/2√3, r)
   end
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, a/2√3, " a/2√3", :red)
       point(x, y, "(x,y)", :red)
       point(x1, y1, "(x1,y1) ", :black, :right)
       point(0, √3a/2, " √3a/2", :black)
       point(a/2, 0, " a/2", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その380)

2023年08月10日 | Julia

算額(その380

東京都八王子市片倉町 片倉・住吉神社 嘉永8年(1851)
http://www.wasan.jp/tokyo/sumiyosi.html
鮮明な写真があります
https://utrblog.exblog.jp/10465003/
解説があります
http://tk2-230-24951.vs.sakura.ne.jp/sangaku/sangaku_kai.pdf

問題

正三角形内に,大円 1 個と小円 2 個がある。正三角形の一辺の長さと,大円の直径が与えられているとき,小円の直径を求めよ。

正三角形の一辺の長さを a
小円の半径と中心座標を r1, (r1, x1)
大円の半径と中心座標を r2, (0, y2)
として,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r1::positive, x1::positive, r2::positive,y2::positive;

eq1 = (sqrt(Sym(3))a/2 - y2)/2 - r2;
eq2 = r1/(a/2 - x1) - 1/sqrt(Sym(3));
eq3 = x1^2 + (y2 - r1)^2 - (r1 + r2)^2;

res = solve([eq1, eq2, eq3], (r1, x1, y2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (sqrt(3)*a/3 - 2*sqrt(r2)*sqrt(sqrt(3)*a - 2*r2)/3 - r2/3, -a/2 + 2*sqrt(3)*sqrt(r2)*sqrt(sqrt(3)*a - 2*r2)/3 + sqrt(3)*r2/3, sqrt(3)*a/2 - 2*r2)
    (sqrt(3)*a/3 + 2*sqrt(r2)*sqrt(sqrt(3)*a - 2*r2)/3 - r2/3, -a/2 - 2*sqrt(3)*sqrt(r2)*sqrt(sqrt(3)*a - 2*r2)/3 + sqrt(3)*r2/3, sqrt(3)*a/2 - 2*r2)

最初の組のものが適解である。

res[1][1] |> factor |> println
res[1][2] |> factor |> println
res[1][3] |> factor |> println

   -(-sqrt(3)*a + 2*sqrt(r2)*sqrt(sqrt(3)*a - 2*r2) + r2)/3
   (-3*a + 4*sqrt(3)*sqrt(r2)*sqrt(sqrt(3)*a - 2*r2) + 2*sqrt(3)*r2)/6
   -(-sqrt(3)*a + 4*r2)/2

   r1 = 0.124024;  a = 1;  r2 = 0.25;  x1 = 0.285184;  y2 = 0.366025
   小円の直径 = 2r1 = 0.248049

using Plots

function draw(a, r2, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x1, y2) = (
       (√3a - 2√r2√(√3a - 2r2) - r2)/3,
       (4√3√r2√(√3a - 2r2) + 2√3r2 - 3a)/6,
       √3a/2 - 2r2
       )
   @printf("r1 = %g;  a = %g;  r2 = %g;  x1 = %g;  y2 = %g\n", r1, a, r2, x1, y2)
   @printf("小円の直径 = 2r1 = %g\n", 2r1)
   plot([-a/2, a/2, 0, -a/2], [0, 0, sqrt(3)a/2, 0], color=:black, linewidth=0.25)
   circle(0, y2, r2)
   circle(x1, r1, r1, :green)
   circle(-x1, r1, r1, :green)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, √3a/2, " √3a/2", :black, :left, :vcenter)
       point(a/2, 0, " a/2", :black, :left, :bottom)
       point(0, y2, " 大円:r2,(0,y2)", :red, :left, :vcenter)
       point(x1, r1, "小円:r1,(x1,r1)", :green, :center, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

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

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