裏 RjpWiki

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

算額(その268)

2023年06月07日 | Julia

算額(その268)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(183)
長野県下高井郡木島平村西小路 満昌院 嘉永5年(1852)

左下の頂点が原点の正方形と,2つの頂点がその正方形の辺上にあり,最初の正方形の右上の頂点が45度回転した正方形がある。左下の二等辺直角三角形中に甲円,右上の2個の二等辺直角三角形中にそれぞれ乙円が内接している。甲円の直径を乙円の直径で表わせ。
甲円と乙円の入っている二等辺直角三角形は相似で,相似比が √2 である。
よって,「甲円の直径 = √2 乙円の直径」である。

   a = 10.000000;  b = 6.666667;  r1 = 1.952621;  r2 = 1.380712;  r1/r2 = 1.4142135623730951

include("julia-source.txt");

using SymPy
@syms b::positive, 甲円の直径::positive, 乙円の直径::positive;

eq1 = 甲円の直径 / 乙円の直径 - sqrt(Sym(2))b / b
res = solve(eq1, 甲円の直径)[1] |> println

   sqrt(2)*乙円の直径

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 10
   b = 2a/3
   r1 = (2b - sqrt(2b^2))/2
   r2 = (sqrt(2)b - b)/2
   @printf("a = %.6f;  b = %.6f;  r1 = %.6f;  r2 = %.6f;  r1/r2 = %.16f\n", a, b, r1, r2, r1/r2)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
   plot!([0, b, 2b, b, 0], [b, 0, b, 2b, b], color=:red, lw=0.5)
   circle(r1, r1, r1, :green)
   circle(a - b/2, a + r2, r2, :magenta)
   if more
       point(a, 0, " a", :black, :left, :bottom)
       point(b, 0, "  b", :black, :left, :bottom)
       point(r1, r1, " 甲:r1\n (r1,r1)", :green, :center, :top)
       point(b, a + r2, "乙:r2\n(b,a+r2)", :magenta, :center, :top)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その267)

2023年06月07日 | Julia

算額(その267)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(64)
長野県下高井郡木島平村往郷 水穂神社 寛政12年(1800)

甲,乙,丙の 3 個の正方形がある。甲,乙の面積の和を A,乙と丙の面積の和を B とする。このとき丙の一辺の長さを求めよ。ただし正方形の一辺の差は等しいとする。

甲,乙,丙の 3 個の正方形の一辺の長さを,甲,乙,丙とする。

include("julia-source.txt");

using SymPy
@syms 甲::positive, 乙::positive, 丙::positive, A::positive, B::positive;

eq1 = 甲^2 + 乙^2 - A
eq2 = 乙^2 + 丙^2 - B
eq3 = (甲 - 乙) - (乙 - 丙)
res = solve([eq1, eq2, eq3], (甲, 乙, 丙));

以下の 2 組の解が得られる。

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

    sqrt(2)*(-A + B - 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B - sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B))
    -(B/2 + sqrt(-A^2 + 6*A*B - B^2)/4)*sqrt(-2*A + 14*B - 2*sqrt(-A^2 + 6*A*B - B^2))/(A - 5*B)
    sqrt(-A/8 + 7*B/8 - sqrt(-A^2 + 6*A*B - B^2)/8)

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

    sqrt(2)*(-A + B + 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B + sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B))
    sqrt(2)*(-A + B + 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B + sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B))
    sqrt(-A/8 + 7*B/8 + sqrt(-A^2 + 6*A*B - B^2)/8)


それぞれを使って,実際に甲,乙,丙を与えて A, B を求め,数式により正しく甲,乙,丙が返ってくるか確かめる。

function func1(甲, 乙, 丙)  # 1 組目の解
    A = 甲^2 + 乙^2
    B = 乙^2 + 丙^2
    return (sqrt(2)*(-A + B - 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B - sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B)),
-(B/2 + sqrt(-A^2 + 6*A*B - B^2)/4)*sqrt(-2*A + 14*B - 2*sqrt(-A^2 + 6*A*B - B^2))/(A - 5*B),
sqrt(-A/8 + 7*B/8 - sqrt(-A^2 + 6*A*B - B^2)/8))
end;

function func2(甲, 乙, 丙)  # 2 組目の解
    A = 甲^2 + 乙^2
    B = 乙^2 + 丙^2
    return (sqrt(2)*(-A + B + 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B + sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B)),
sqrt(2)*(-A + B + 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B + sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B)),
sqrt(-A/8 + 7*B/8 + sqrt(-A^2 + 6*A*B - B^2)/8))
end;

func1(8, 6, 4) |> println

func1(9.34, 8.084, 6.828) |> println

    (8.000000000000002, 6.0, 4.0)
    (9.340000000000003, 8.084000000000001, 6.828000000000001)

func2(8, 6, 4) |> println
func2(9.34, 8.084, 6.828) |> println

    (-9.899494936611665, -9.899494936611665, 7.0710678118654755)
    (-12.32062855539441, -12.32062855539441, 10.544376321053798)


1 組目の解 func1 のほうが適切解である。

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

算額(その266)

2023年06月07日 | Julia

算額(その266)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(171)
長野県下高井郡木島平村上木島 天然寺 弘化2年(1845)

図のように,正三角形の内側に 3 個の正方形がある。それぞれの正方形の 2 つの頂点は正三角形の辺の上にあり,他の 2 個の頂点は他の正方形の頂点と一致している。正三角形の一辺が 8 寸のとき,正方形の一辺の長さを求めよ。

正三角形の一辺の長さを a,正方形の一辺の長さを x とすると,⊿ABC は一辺の長さが x の正三角形,∠ OAD は 30°,OB = a/2 ゆえ 2(x + √3/2x) = a。x について解くと x = (2 - √3)a。a に 8 を代入すると x = 2.14359353944898。

include("julia-source.txt");

using SymPy
@syms x::positive, a::positive;

eq1 = 2(x + sqrt(Sym(3))/2 * x) - a
res = solve(eq1, x)[1]
res |> println

   -sqrt(3)*a + 2*a

res(a => 8).evalf()

   2.14359353944898

正方形の一辺の長さは 2.143593 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 8
   x = (2 - √3)a
   @printf("正三角形の一辺の長さ = %.6f;  正方形の一辺の長さ = %.6f\n", a, x)
   plot([a/2, 0, -a/2, a/2], [0, a√3/2, 0, 0], color=:blue, lw=0.5)
   xy = x .* [ # 図形の頂点座標,一筆書き
           0                        sind(30)                 # 1
           cosd(60)                 sind(30)+sind(60)        # 2
           cosd(30) + cosd(60)      sind(60)                 # 3
           cosd(30)                 0                        # 4
           0                        sind(30)                 # 5
           -cosd(30)                0                        # 6
           -cosd(30) - cosd(60)     sind(60)                 # 7
           -cosd(60)                sind(30) + sind(60)      # 8
           -cosd(60)                sind(30) + sind(60) + 1  # 9
           cosd(60)                 sind(30) + sind(60) + 1  # 10
           cosd(60)                 sind(30) + sind(60)      # 11
           -cosd(60)                sind(30) + sind(60)      # 12
           0                        sind(30)                 # 13
       ]
   plot!(xy[:, 1], xy[:, 2], color=:red, lw=0.5)
   if more
       point(0, 0, " O", :red, :left, :bottom)
       point(x*cosd(30), 0, " A", :red, :left, :bottom)
       point(x*cosd(30) + x, 0, " B", :red, :left, :bottom)
       point(x*(cosd(30)+cosd(60)), x*sind(60), " C", :red, :left, :bottom) 
       point(0, x*sind(30), "  D", :red, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

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

算額(その265)

2023年06月07日 | Julia

算額(その265)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3(91)
長野県下高井郡木島平村穂高 一川谷大元神社 文化8年(1811)

直角三角形(鉤股弦)内に円と正方形がある。三角形の面積は 294 平方寸 ,鉤の3乗と股の3上の和に円の直径(径)と正方形の一辺の長さ(方面)の差を掛けると 62426 となる(数字の掛け算で次元は無視)。このとき,円の直径を求めよ。

include("julia-source.txt");

using SymPy
@syms 鉤::positive, 股::positive, 弦::positive,
   径::positive, 方面::positive;

eq1 = (鉤 + 股 - 弦) - 径
eq2 = (鉤 * 股 / 2) - 294;
eq3 = (鉤^3 + 股^3)*(径 - 方面) - 62426
eq4 = (鉤^2 + 股^2) - 弦^2
eq5 = (鉤 - 方面)/方面 - 方面/(股 - 方面)
res = solve([eq1, eq2, eq3, eq4, eq5], (鉤, 股, 弦, 径, 方面))

   2-element Vector{NTuple{5, Sym}}:
    (21, 28, 35, 14, 12)
    (28, 21, 35, 14, 12)

鉤 < 股 なので (21, 28, 35, 14, 12) が適切である。

鉤 = 21;  股 = 28;  弦 = 35;  径 = 14;  方面 = 12

円の直径は 14 である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, 弦, 径, 方面) = (21, 28, 35, 14, 12)
   r = 径/2
   @printf("鉤 = %.0f;  股 = %.0f;  弦 = %.0f;  径 = %.0f;  方面 = %.0f\n", 鉤, 股, 弦, 径, 方面)
   plot([0, 股, 0, 0], [0, 0, 鉤, 0], color=:black, lw=0.5)
   plot!([0, 方面, 方面, 0, 0], [0, 0, 方面, 方面, 0], color=:black, lw=0.5)
   circle(径/2, 径/2, 径/2)
   if more
       point(r, r, "半径=$(Int(r))")
       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アクセスランキング にほんブログ村