裏 RjpWiki

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

算額(その1121)

2024年07月04日 | Julia

算額(その1121)

二十八 一関市萩荘 赤萩観音寺 天保2年(1831)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html
キーワード:円1個,四分円(半円),斜線3本

正方形の中に斜線3本と,円弧,小円を容れる。円弧の直径が 1 寸のとき,小円の直径はいかほどか。

注:「問」では「大半円」といっているが,明らかに半円ではなく,四分円である。

正方形の一辺の長さを a
大円の半径と中心座標を R, (a/2, a + a/√2); a = √2R
小円の半径と中心座標を r, (x, r)
とおき,以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms a::positive, R::positive, r::positive, x::positive
a = √Sym(2)R
eq1 = dist2(0, a, a, 0, x, r, r)
eq2 = dist2(0, a, a/2, 0, x, r, r)
res = solve([eq1, eq2], (r, x))[2]

   (-9*sqrt(2)*R*(-1 + sqrt(13 - 4*sqrt(10))/3)^2/8 + sqrt(2)*R*(-1 + sqrt(13 - 4*sqrt(10))/3)/2 - 9*sqrt(2)*R*(-1 + sqrt(13 - 4*sqrt(10))/3)^3/16 + sqrt(2)*R, -3*sqrt(2)*R*(-1 + sqrt(13 - 4*sqrt(10))/3)/4)

ans_r = res[1] |> sympy.sqrtdenest |> simplify
ans_r |> println

   R*(-2*sqrt(5) - 2 + sqrt(10) + 3*sqrt(2))/4

小円の半径 r は,大円の半径 R の (√10 + 3√2 - 2√5 - 2)/4 = 0.2331955980720215 倍である。
大円の直径が 1 寸のとき,小円の直径は 0.2331955980720215 である。

「術」(山村も)では 1 - √(1 - √0.5) = 0.4588038998538031 倍としているが,根拠も不明である。そもそも,図を描いてみても,大円のほぼ半分の大きさというのは理解できない。

function draw(R, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = √2R
   (r, x) = (-9*sqrt(2)*R*(-1 + sqrt(13 - 4*sqrt(10))/3)^2/8 + sqrt(2)*R*(-1 + sqrt(13 - 4*sqrt(10))/3)/2 - 9*sqrt(2)*R*(-1 + sqrt(13 - 4*sqrt(10))/3)^3/16 + sqrt(2)*R, -3*sqrt(2)*R*(-1 + sqrt(13 - 4*sqrt(10))/3)/4)
   @printf("大円の直径が %g のとき,小円の直径は %g である。\n", 2R, 2r)
   #@printf("r2 = %g;  a = %g;  b = %g;  c = %g;  r1 = %g;  x1 = %g;  y2 = %g\n", r2, a, b, c, r1, x1, y2)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:magenta, lw=0.5)
   circle(a/2, a + R/√2, R)
   circle(x, r, r, :green)
   segment(0, a, a, 0, :blue)
   segment(0, 0, a, a, :blue)
   segment(0, a, a/2, 0, :blue)
   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(a, 0, " a", :magenta, :left, :bottom, delta=delta/2)
       point(x, r, "小円:r\n(x,r)", :green, :center, delta=-delta/2)
       point(a/2, a + R/√2, "大円:R,(a/2,a+R/√2)", :red, :center, delta=-delta/2)
   end
end;

 

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

算額(その1119)

2024年07月04日 | Julia

算額(その1119)

十 岩手県胆沢町若柳市野々 個人宅 安政2年(1855)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html
キーワード:円3個,正方形,斜線2本

正方形の中に斜線 2 本を引き,甲円 2 個,乙円 1 個を容れる。乙円の直径が 1 寸のとき,短い方の斜線の長さはいかほどか。

正方形の一辺の長さを a,斜線と正方形の斜線との交点を (0, b), (c, a)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (a - r2, y2)
とおき,以下の連立方程式を解く。
短い方の斜線の長さは sqrt((a - c)^2  + a^2) である。

include("julia-source.txt")

using SymPy
@syms a::positive, b::positive, c::positive,
     r1::positive, x1::positive,
     r2::positive, y2::positive
eq1 = dist2(0, b, a, a, x1, r1, r1)
eq2 = dist2(0, b, a, a, r1, a - r1, r1)
eq3 = dist2(0, b, a, a, a - r2, y2, r2)
eq4 = dist2(c, a, a, 0, x1, r1, r1)
eq5 = dist2(c, a, a, 0, r1, a - r1, r1)
eq6 = dist2(c, a, a, 0, a - r2, y2, r2);

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (a, b, c, r1, x1, y2) = u
   return [
       a^2*b^2 - 2*a^2*b*r1 + 2*a^2*b*x1 - a^2*r1^2 - 2*a^2*r1*x1 + a^2*x1^2 - 2*a*b^2*x1 + 2*a*b*r1^2 + 2*a*b*r1*x1 - 2*a*b*x1^2 - b^2*r1^2 + b^2*x1^2,  # eq1
       a*(a^3 - 2*a^2*b - 4*a^2*r1 + a*b^2 + 6*a*b*r1 + 2*a*r1^2 - 2*b^2*r1 - 2*b*r1^2),  # eq2
       a*(a^3 - 2*a^2*r2 - 2*a^2*y2 + 2*a*b*r2 - a*r2^2 + 2*a*r2*y2 + a*y2^2 - 2*b*r2*y2),  # eq3
       a*(a^3 - 2*a^2*r1 - 2*a^2*x1 + 2*a*c*r1 - a*r1^2 + 2*a*r1*x1 + a*x1^2 - 2*c*r1*x1),  # eq4
       a*(a*c^2 - 2*a*r1^2 - 2*c^2*r1 + 2*c*r1^2),  # eq5
       -a^2*r2^2 - 2*a^2*r2*y2 + a^2*y2^2 + 2*a*c*r2^2 + 2*a*c*r2*y2 - 2*a*c*y2^2 - c^2*r2^2 + c^2*y2^2,  # eq6
   ]
end;

r2 = 1/2
iniv = BigFloat[2.9, 0.5, 1.2, 0.8, 1.6, 1.9]
res = nls(H, ini=iniv)

   ([2.927373520970951, 0.4800319114286444, 1.2284447668135445, 0.7795448985438439, 1.5736412260872523, 1.8576489200442854], true)

乙円の直径が 1 寸のとき,小斜は 3.3846528098144275 寸である。

「答」では 2寸3分6厘となっているが,そんな長さはどこにもない。

sqrt(2.927373520970951^2 + (2.927373520970951 - 1.2284447668135445)^2)

   3.3846528098144275

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

   r2 = 0.5;  a = 2.92737;  b = 0.480032;  c = 1.22844;  r1 = 0.779545;  x1 = 1.57364;  y2 = 1.85765

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (a, b, c, r1, x1, y2) = res[1]
   小斜 = sqrt(a^2 + (a - c)^2)
   @printf("乙円の直径が %g のとき,小斜の長さは %g である。\n", 2r2, 小斜)
   @printf("r2 = %g;  a = %g;  b = %g;  c = %g;  r1 = %g;  x1 = %g;  y2 = %g\n", r2, a, b, c, r1, x1, y2)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
   circle(x1, r1, r1)
   circle(r1, a - r1, r1)
   circle(a - r2, y2, r2, :blue)
   segment(0, b, a, a, :magenta)
   segment(c, a, a, 0, :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(0, b, "b ", :magenta, :right, :vcenter)
       point(0, a, "a ", :magenta, :right, :vcenter)
       point(a, 0, " a", :magenta, :center, delta=-delta/2)
       point(c, a, "(c,a)", :magenta, :center, :bottom, delta=delta/2)
       point(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta/2)
       point(r1, a - r1, "甲円:r1,(r1,a-r1)", :red, :center, delta=-delta/2)
       point(a - r2, y2, "乙円:r2,(a-r2,y2)", :blue, :center, delta=-delta/2)
       plot!(xlims=(-6delta, a+3delta), ylims=(-4delta, a+3delta))
   end
end;

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

算額(その1118)

2024年07月04日 | Julia

算額(その1118)

八 岩手県水沢市南津田字化粧坂 中尊寺末寺薬師堂 明治18年(1885)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html
キーワード:円,直角三角形,正方形,矢,鈎,股,弦

算額では基本的な三題である。

1. 直角三角形内の正方形

鈎 8.4 寸,股 14 寸のとき,内接する正方形の一辺の長さ(方面)を求めよ。

include("julia-source.txt")

using SymPy
@syms 鈎, 股, 方面
eq1 = (鈎 - 方面)/方面 - 鈎/股
res1 = solve(eq1, 方面)[1]
res1 |> println
res1(鈎 => 8.4, 股 => 11.2).evalf() |> println

   股*鈎/(股 + 鈎)
   4.80000000000000

正方形の一辺の長さは,股*鈎/(股 + 鈎) で得られる。

2. 直角三角形内の円

鈎 8.4 寸,弦 14 寸のとき,内接する円の径(直径)を求めよ。

@syms 鈎, 股, 弦, 径
股 = sqrt(弦^2 - 鈎^2)
eq5 = 鈎 + 股 - 弦 - 径
res2 = solve(eq5, 径)[1]
res2 |> println
股 |> println
股(鈎 => 8.4, 弦 => 14).evalf() |> println
res2(鈎 => 8.4, 弦 => 14).evalf() |> println

   -弦 + 鈎 + sqrt(弦^2 - 鈎^2)
   sqrt(弦^2 - 鈎^2)
   11.2000000000000
   5.60000000000000

円の直径は 鈎 + 股 - 弦 で得られる。

3. 弦と矢

直径が 10 寸,弦が 6 寸のとき,矢を求めよ。

@syms 径, 矢, 弦

eq6 = (径/2)^2 - ((径/2) - 矢)^2 - (弦/2)^2
res3 = solve(eq6, 矢)[1]
res3 |> println
res3(径 => 10, 弦 => 6) |> println

   径/2 - sqrt(-弦^2 + 径^2)/2
   1

矢は 径/2 - sqrt(-弦^2 + 径^2)/2 で得られる。

4. 図

function ex1()
   (鈎, 弦) = (8.4, 11.2)
   (方面, 股) = (鈎*(弦^2 - 鈎^2 - 鈎*sqrt(弦^2 - 鈎^2))/(弦^2 - 2*鈎^2), sqrt((弦 - 鈎)*(弦 + 鈎)))
   p1 = plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5, showaxis=false)
   delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
   plot!([0, 方面, 方面, 0,  0], [0, 0, 方面, 方面, 0], color=:red, lw=0.5)
   dimension_line(-8delta, 0, -8delta, 鈎, "鈎", deltax=-6delta)
   dimension_line(0, -8delta, 股, -8delta, "股", delta=-6delta)
   plot!(xlims=(-20delta, 股 + 3delta))
   return p1
end;

function ex2()
   (鈎, 弦) = (8.4, 14)
   股 = sqrt(弦^2 - 鈎^2)
   径 = -弦 + 鈎 + sqrt(弦^2 - 鈎^2)
   p2 = plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5, showaxis=false)
   delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
   r = 径/2
   circle(r, r, r)
   dimension_line(-10delta, 0, -10delta, 鈎, "鈎", deltax=-6delta)
   dimension_line(0, -10delta, 股, -10delta, "股", delta=-6delta)
   return p2
end;

function ex3()
   (径, 弦) = (10, 6)
   r = 径/2
   矢 = 径/2 - sqrt(-弦^2 + 径^2)/2
   p3 = plot(showaxis=false)
   delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
   circle(0, 0, r)
   x = sqrt(r^2 - (r - 矢)^2)
   segment(-x, r - 矢, x, r - 矢, :green)
   segment(-x, r - 矢, x, 矢 - r)
   segment(0, r, 0, r - 矢, :blue)
   point(0, r - 矢/2, "矢", :blue, :left, :vcenter, mark=false)
   point(0, r - 矢, "弦", :green, :left, delta=-5delta, deltax=100delta, mark=false)
   point(0, 0, "径", :black, :left, :vcenter, deltax=20delta, mark=false)
   return p3
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   p1 = ex1()
   p2 = ex2()
   p3 = ex3()
   plot(p1, p2, p3, layout=(1,3))
   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)
   end
end;

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

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

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