裏 RjpWiki

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

算額(その445)

2023年09月26日 | Julia

算額(その445)

十二 群馬県高崎市石原町 清水寺 享和3年(1803)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.

享和3年癸亥4月(1803) 上州中里邑 黒崎九兵祀則
藤田貞資(1807):続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

キーワード:円4個,外円,菱形
#Julia, #SymPy, #算額, #和算

外円内に甲円 1 個,乙円 2 個,菱形 1 個が入っている。外円と甲円の直径がそれぞれ 7 寸,3 寸のとき,乙円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (r0 - r1, 0)
乙円の半径と中心座標を r2, (x2, y2)
とおき,連立方程式の解を求める。

include("julia-source.txt")

using SymPy

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

(r0, r1) = (7, 3) .// 2
eq1 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + y2^2 - (r0 - r2)^2
eq3 = distance(0, r0 - 2r1, sqrt(r0^2 - r1^2), -r1, x2, y2) - r2^2 
res = solve([eq1, eq2, eq3], (r2, x2, y2));

解は 1 組得られるが,係数が分数の長い式になるので,evalf() した値を示しておく。

乙円の直径は 2 * 1.2 = 2寸4分である。

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

   (1.20000000000000, 2.24499443206436, 0.500000000000000)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1) = (7, 3) .// 2
   (r2, x2, y2) = (res[1][1].evalf(), res[1][2].evalf(), res[1][3].evalf())
   @printf("r2 = %g;  x2 = %g;  y2 = %g\n", r2, x2, y2)
   @printf("乙円の直径 = %g\n", 2r2)
   plot()
   circle(0, 0, r0, :red)
   circle(0, r0 - r1, r1, :blue)
   circle2(x2, y2, r2, :brown)
   b = sqrt(r0^2 - r1^2)
   plot!([0, b, 0, -b, 0], [-r0, -r1, r0 - 2r1, -r1, -r0], color=:orange, lw=0.5)
   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, r0 - r1, " 甲円:r1, (0,r0-r1)", :black, :center, delta=-delta)

       point(x2, y2, "乙円:r2,(x2,y2)", :brown, :center, delta=-delta)
       point(0, r0 - 2r1, " r0-2r1")
       point(0, -r1, " -r1")
       point(0, r0, " r0", :red, :left, :bottom)
       point(sqrt(r0^2 - r1^2), -r1, "(√(r0^2-r1^2),-r1) ", :orange, :right, :vcenter)
   end
end;

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

算額(その444)

2023年09月26日 | Julia

算額(その444)

十二 群馬県高崎市石原町 清水寺 享和3年(1803)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.

享和3年癸亥4月(1803) 上州棟高邑 坂本丹治亮春
藤田貞資(1807):続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

キーワード:円5個,長方形

長方形内に大円 1 個,中円,小円がそれぞれ 2 個入っており,互いに内接・外接しあっている。長方形の長辺が 40寸8分のとき,それぞれの円の直径を求めよ。

原点を図のようにおき,長方形の長辺を a,短辺を 2b とする。
大円の半径と中心座標を r1, (a - r1, 0)
中円の半径と中心座標を r2, (r2, b - r2)
小円の半径と中心座標を r3, (r3, r3)
とおき,連立方程式の解を求める。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive,
     r1::positive, r2::positive, r3::positive;

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

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

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

a = 40寸8分 のとき各円の直径は 34寸1分あまり, 10寸2分, 7寸あまりである。

2(a*(9 - 4*sqrt(2))/8)(a => 40.8).evalf() |> println
2(a/8)(a => 40.8).evalf() |> println
2(a*(3/2 - sqrt(2)))(a => 40.8).evalf() |> println

   34.1000866551777
   10.2000000000000
   7.00017331035544

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 408//10
   (r1, r2, r3) = (a*(9 - 4*sqrt(2))/8, a/8, a*(3/2 - sqrt(2)))
   b = r1
   @printf("r1 = %g;  r2 = %g;  r3 = %g\n", r1, r2, r3)
   plot([0, a, a, 0, 0], [-b, -b, b, b, -b], color=:black, lw=0.5)
   circle(a - r1, 0, r1, :red)
   circle(r2, b - r2, r2, :blue)
   circle(r2, r2 - b, r2, :blue)
   circle(r3, r3, r3, :brown)
   circle(r3, -r3, r3, :brown)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(a - r1, 0, " 大円:r1,(a-r1,0)", :red, delta=-delta)
       point(r2, b - r2, " 中円:r2,(r2,b-r2)", :black)
       point(r3, r3, " 小円:r3,(r3,r3)", :black)
       point(a, 0, "a ", :black, :right, :bottom, delta=delta)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

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

算額(その443)

2023年09月26日 | Julia

算額(その443)

天明五年乙巳秋九月(1785) 藤田貞資門人 筑後州久留米藩 清水與市道香

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

求めるものと条件として与えられる数値が異なるだけで,算額(その31)と基本的に同じである。

図のように,直線の上に大円,中円,小円,甲円,乙円がそれぞれ外接しあっている。甲円,乙円の径を98寸1厘,121寸 としたとき,中円の径を求めよ。

それぞれの円の中心座標と半径を (0, r1, r1), (x2, r2, r2), (x3, r3, r3), (x4, r4, r4), (x5, r5, r5) とする。

include("julia-source.txt")

using SymPy

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

(r4, r5) = (981, 1210) .// 20;
eq1 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2;
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2;
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2;
eq4 = x5^2 + (r1 - r5)^2 - (r1 + r5)^2;
eq5 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2;
eq6 = (x2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2;
eq7 = (x5 - x3)^2 + (r3 - r5)^2 - (r3 + r5)^2;

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r2, r3, x2, x3, x4, x5))

   2-element Vector{NTuple{7, Sym}}:
    (3917133*sqrt(1090)/3682898 + 304705467/7365796, 121/2, 3917133*sqrt(1090)/52441 + 260073891/104882, 3993*sqrt(1090)/2714 + 118701/1357, 188259786/310753 + 11751399*sqrt(1090)/621506, 118701/1357 + 32373*sqrt(1090)/6785, 3993*sqrt(1090)/2714 + 118701/1357)
    (35254197*sqrt(1090)/3682898 + 2742349203/7365796, 70508394*sqrt(1090)/14891881 + 6218626689/29783762, 2340665019/104882 - 35254197*sqrt(1090)/52441, 2340665019/5236663 + 176270985*sqrt(1090)/10473326, -401684184/310753 + 35254197*sqrt(1090)/621506, 356103/1357 + 97119*sqrt(1090)/6785, 11979*sqrt(1090)/2714 + 356103/1357)

2番めのものが適解である。

name = ["r1", "r2", "r3", "x2", "x3", "x4", "x5"]
j = 2
for i in 1:7
   println("$(name[i]) = $(res[j][i].evalf())")
end

   r1 = 688.343020749222
   r2 = 365.108908025960
   r3 = 122.232157448001
   x2 = 1002.63686078867
   x3 = 580.129821644955
   x4 = 734.990886123080
   x5 = 408.140920542540

中円の直径は 730.21781605192 である。算額の答えでは 729 となっており,かなり違いがある。

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

   730.21781605192

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r4, r5) = (981, 1210) ./ 20
   r1, r2, r3, x2, x3, x4, x5 = (35254197*sqrt(1090)/3682898 + 2742349203/7365796, 70508394*sqrt(1090)/14891881 + 6218626689/29783762, 2340665019/104882 - 35254197*sqrt(1090)/52441, 2340665019/5236663 + 176270985*sqrt(1090)/10473326, -401684184/310753 + 35254197*sqrt(1090)/621506, 356103/1357 + 97119*sqrt(1090)/6785, 11979*sqrt(1090)/2714 + 356103/1357)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  r4 = %g;  r5 = %g\n", r1, r2, r3, r4, r5)
   @printf("中円の直径 = %g\n", 2r2)
   plot()
   circle(0, r1, r1, :red)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :brown)
   circle(x4, r4, r4, :green)
   circle(x5, r5, r5, :red)
   hline!([0], color=:black, linewidth=0.25)
   if more
       point(0, r1, "大円 r1")
       point(x2, r2,"中円 r2")
       point(x3, r3,"小")
       point(x4, r4,"乙")
       point(x5, r5,"甲")
   end
end;

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

Plotlyを使用した対話型データ可視化の作成

2023年09月26日 | Python
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

NumPyとPandasを組み合わせて高速なデータ操作を実現する方法

2023年09月26日 | Python
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

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

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