裏 RjpWiki

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

算額(その924)

2024年05月06日 | Julia

算額(その924)

香川県善通寺市与北町 皇美屋神社 明治11年(1878)
本田益夫:算額随想-香川県内の算額について-,私家版,昭和45年(1970).

直角三角形内に菱形を入れる。鈎,股,弦の長さがそれぞれ 3 寸,4 寸,5 寸のとき,菱形の一辺の長さはいかほどか。

鈎,股,弦の長さをそのまま変数名として使う。
菱形の左側の一辺が鈎,股と交差する座標を (0, b), (a, 0)
菱形の一辺の長さを diamond
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     a::positive, b::positive, diamond::positive
eq1 = a^2 + b^2 - diamond^2
eq2 = (鈎 - b)/diamond - 鈎/股
eq3 = b/a - 鈎/股
(a, b, diamond) = solve([eq1, eq2, eq3], (a, b, diamond))[1];

菱形の一辺の長さは 股*(股^2 - 股*sqrt(股^2 + 鈎^2) + 鈎^2)/鈎^2 である。

diamond |> println

   股*(股^2 - 股*sqrt(股^2 + 鈎^2) + 鈎^2)/鈎^2

鈎,股,弦の長さがそれぞれ 3 寸,4 寸,5 寸のとき,菱形の一辺の長さは 20/9 = 2.22222222222222 寸である。

ちなみに,a = 16/9 = 1.77777777777778 寸, b = 4/3 = 1.33333333333333 寸である。

a(鈎 => 3, 股 => 4, 弦 => 5).evalf() |> println
b(鈎 => 3, 股 => 4, 弦 => 5).evalf() |> println

   1.77777777777778
   1.33333333333333

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, 弦) = (3, 4, 5)
   (a, b, diamond) = (股^2*(-股 + sqrt(股^2 + 鈎^2))/鈎^2, 股*(-股 + sqrt(股^2 + 鈎^2))/鈎, 股*(股^2 - 股*sqrt(股^2 + 鈎^2) + 鈎^2)/鈎^2)
   @printf("鈎 %g;  股 = %g;  弦 = %g のとき, 菱形の一辺の長さ = %g\n", 鈎, 股, 弦, diamond)
   @printf("a = %g;  b = %g;  diamond = %g\n", a, b, diamond)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   plot!([0, a, 股, diamond, 0], [b, 0, 0, b, b], color=:red, lw=1)
   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, " 股", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, "  a", :blue, :left, :bottom, delta=delta/2)
       point(0, 鈎, " 鈎", :blue, :left, :bottom, delta=delta/2)
       point(0, b, "  b", :blue, :left, :vcenter)
   end
end;

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

算額(その923)

2024年05月06日 | Julia

算額(その923)

香川県善通寺市与北町 皇美屋神社 明治11年(1878)
本田益夫:算額随想-香川県内の算額について-,私家版,昭和45年(1970).

鈎股弦内に正方形を入れる。鈎,股,弦が 12,16,25 寸のとき,正方形の一辺の長さを求めよ。



変数名をそのまま鈎,股,弦とする。
正方形の一辺の長さを a
鈎,股の辺上にある正方形の頂点座標を (0, c), (b ,0) とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     a::positive, b::positive, c::positive
eq1 = c/b - 鈎/股
eq1 = b^2 + c^2 - a^2
eq2 = a/(鈎 - c) - 股/弦
eq3 = a/(股 - b) - 鈎/弦
(a, b, c) = solve([eq1, eq2, eq3], (a, b, c))[1]

   (股*鈎*(弦*股^2 + 弦*鈎^2 - 股*鈎*sqrt(股^2 + 鈎^2))/(弦^2*股^2 + 弦^2*鈎^2 - 股^2*鈎^2), 股^2*鈎*(弦*sqrt(股^2 + 鈎^2) - 股*鈎)/(弦^2*股^2 + 弦^2*鈎^2 - 股^2*鈎^2), 股*鈎^2*(弦*sqrt(股^2 + 鈎^2) - 股*鈎)/(弦^2*股^2 + 弦^2*鈎^2 - 股^2*鈎^2))

鈎,股,弦が 12,16,25 寸のとき,正方形の一辺の長さは 寸である。

a(鈎 => 12, 股 => 16, 弦 => 20).evalf() |> println
b(鈎 => 12, 股 => 16, 弦 => 20).evalf() |> println
c(鈎 => 12, 股 => 16, 弦 => 20).evalf() |> println

   6.48648648648649
   5.18918918918919
   3.89189189189189

正方形の一辺の長さ a の一般式はかなり長い。

a |> println
b |> println
c |> println

   股*鈎*(弦*股^2 + 弦*鈎^2 - 股*鈎*sqrt(股^2 + 鈎^2))/(弦^2*股^2 + 弦^2*鈎^2 - 股^2*鈎^2)
   股^2*鈎*(弦*sqrt(股^2 + 鈎^2) - 股*鈎)/(弦^2*股^2 + 弦^2*鈎^2 - 股^2*鈎^2)
   股*鈎^2*(弦*sqrt(股^2 + 鈎^2) - 股*鈎)/(弦^2*股^2 + 弦^2*鈎^2 - 股^2*鈎^2)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, 弦) = (12, 16, 20)
   # (鈎, 股, 弦) = (8, 15, 17)
   t = 弦^2*股^2 + 弦^2*鈎^2 - 股^2*鈎^2
   u = sqrt(股^2 + 鈎^2)
   v = 股*鈎*(弦*u - 股*鈎)/t
   (a, b, c) = (
       股*鈎*(弦*股^2 + 弦*鈎^2 - 股*鈎*u)/t,
       股*v,
       鈎*v)
   @printf("鈎 %g;  股 = %g;  弦 = %g のとき,正方形の一辺の長さ = %g\n", 鈎, 股, 弦, a)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   plot!([b, b + a*(鈎/弦), b - a*(股 - 鈎)/弦, 0, b], [0, a*(股/弦), a*(鈎 + 股)/弦, c, 0], color=:red, 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, " 股", :blue, :left, :bottom, delta=delta/2)
       point(b, 0, "  b", :blue, :left, :bottom, delta=delta/2)
       point(0, 鈎, " 鈎", :blue, :left, :bottom, delta=delta/2)
       point(0, c, "  c", :blue, :left, :vcenter)
   end
end;

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

算額(その922)

2024年05月06日 | Julia

算額(その922)

香川県善通寺市与北町 皇美屋神社 明治11年(1878)
本田益夫:算額随想-香川県内の算額について-,私家版,昭和45年(1970).

大円内に中円 2 個と小円 1 個を入れる。中円の直径が 8 寸,小円の直径が 4 寸のとき,大円の直径を求めよ。

大円の半径と中心座標を R, (0, 0)
中円の半径と中心座標を r1, (0, R - r1)
小円の半径と中心座標を r2, (0, 0)
とおいて,方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive
eq1 = R - 2r1 + r2;
R = solve(eq1, R)[1]
R |> println

   2*r1 - r2

大円の半径 R は 2*r1 - r2 である。
中円の直径が 8 寸,小円の直径が 4 寸のとき,大円の直径は 12 寸である。

2R(r1 => 8/2, r2 => 4/2) |> println

   12.0000000000000

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (8, 4) .// 2
   R = 2*r1 - r2
   @printf("中円の直径が %g,小円の直径が %g のとき,大円の直径は %g である\n", 2r1, 2r2, 2R)
   plot()
   circle(0, 0, R)
   circle2(R - r1, 0, r1, :blue)
   circle(0, 0, r2, :green)
   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(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(R - r1, 0, " 中円:r1,(R-r1,0)", :blue, :left, delta=-delta/2)
       point(0, 0, "小円:r2,(0,0)", :green, :center, :bottom, delta=delta/2)
   end
end;

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

算額(その921)

2024年05月06日 | Julia

算額(その921)

香川県善通寺市与北町 皇美屋神社 明治11年(1878)
本田益夫:算額随想-香川県内の算額について-,私家版,昭和45年(1970).

大円内に正方形と小円を入れる。正方形の一辺の長さが 6 寸,小円の直径が 3 寸のとき,大円の直径を求めよ。

本算額の図形は基本的に「 算額(その877)」と同じであるが,条件の与え方が違うので別途取り上げる。

大円の半径と中心座標を R, (0, 0)
小円の半径と中心座標を r, (0, R - r)
正方形の一辺の長さを 2a
とおいて,方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, a::positive
eq1 = a^2 + (R - 2r - 2a)^2 - R^2;
res = solve(eq1, R)[1]
res |> println

   (a^2/4 + (a + r)^2)/(a + r)

大円の半径 R は (a^2/4 + (a + r)^2)/(a + r) である。
正方形の一辺の長さ 2a が 6 寸,小円の直径 2r が 3 寸のとき,大円の直径は 10 寸である。

2res(a => 6/2, r => 3/2) |> println

   10.0000000000000

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r) = (6, 3) .// 2
   R = (a^2/4 + (a + r)^2)/(a + r)
   @printf("正方形の一辺の長さが %g,小円の直径が %g のとき,大円の直径は %g である\n", 2a, 2r, 2R)
   plot([a, a, -a, -a, a], [R - 2r - 2a, R - 2r, R - 2r, R - 2r - 2a, R - 2r - 2a], color=:green, lw=0.5)
   circle(0, 0, R)
   circle(0, R - r, r, :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(0, R, " R", :red, :left, :bottom, delta=delta/2)
       point(0, R - r, " 小円:r2,(0,R-r)", :red, :center, delta=-delta/2)
       point(0, R - 2r, " R-2r", :green, :left, delta=-delta/2)
       point(0, R - 2r - 2a, " R-2r-2a", :green, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その920)

2024年05月06日 | ブログラミング

算額(その920)

香川県善通寺市与北町 皇美屋神社 明治11年(1878)
本田益夫:算額随想-香川県内の算額について-,私家版,昭和45年(1970).

直線上に大円 2 個,小円 3 個,内円 1 個が積み上がっている。大円,小円の直径がそれぞれ 14 寸,7 寸のとき,内円の直径はいかほどか。

大円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (r2, y2), (0, y22)
内円の半径と中心座標を r3, (0, y3)
とおき,以下の連立方程式を解く。

大円,小円の直径が特定の値のときには,直接数値を代入して連立方程式を解けばよい。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, y2::positive, y22::positive,
     r3::positive, y3::positive
(r1, r2) = (14, 7) .// 2
eq1 = (r1 - r2)^2 + (y2 - r1)^2 - (r2 + r1)^2
eq2 = r1^2 + (y3 - r1)^2 - (r3 + r1)^2
eq3 = r2^2 + (y2 - y3)^2 - (r2 + r3)^2
eq4 = r2^2 + (y22 - y2)^2 - 4r2^2;
res = solve([eq1, eq2, eq3, eq4], (y2, y22, r3, y3))[2]  # 2 of 2

   (7 + 7*sqrt(2), 7*sqrt(3)/2 + 7 + 7*sqrt(2), 2, 4*sqrt(2) + 7)

大円,小円の直径がそれぞれ 14 寸,7 寸のとき,内円の直径は 4 寸である。

---

以下では,一般解を求める。

@syms r1::positive, r2::positive, y2::positive, y22::positive,
     r3::positive, y3::positive
eq1 = (r1 - r2)^2 + (y2 - r1)^2 - (r2 + r1)^2 |> expand
eq2 = r1^2 + (y3 - r1)^2 - (r3 + r1)^2 |> expand
eq3 = r2^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand
eq4 = r2^2 + (y22 - y2)^2 - 4r2^2 |> expand;

y2, y3 は簡単に求めることができる。

println(eq1, " から y2 を求める")
ans_y2 = solve(eq1, y2)[2]  # 2 of 2
ans_y2 |> println

   r1^2 - 4*r1*r2 - 2*r1*y2 + y2^2 から y2 を求める
   2*sqrt(r1)*sqrt(r2) + r1

println(eq2, " から y3 を求める")
ans_y3 = solve(eq2, y3)[2]  # 2 of 2
ans_y3 |> println

   r1^2 - 2*r1*r3 - 2*r1*y3 - r3^2 + y3^2 から y3 を求める
   r1 + sqrt(r3)*sqrt(2*r1 + r3)

r3 は,上で求めた y2, y3 を eq3 に代入し,方程式を解けばよい。もう少し簡約化することもできるが,大差ない。

eq13 = eq3(y2 => ans_y2, y3 => ans_y3) |> expand
println(eq13, " から r3 を求める")
ans_r3 = solve(eq13, r3)[1]
ans_r3 |> println

   -4*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(2*r1 + r3) + 4*r1*r2 + 2*r1*r3 - 2*r2*r3 から r3 を求める
   2*(-2*sqrt(2)*r1^(3/2)*r2^(3/2) + r1*r2*(r1 + r2))/(r1^2 - 6*r1*r2 + r2^2)

y22 は eq4 を解き,式中の y2 に,すでに求めた ans_y2 を代入すればよい。

println(eq4, " から y22 を求める")
ans_y22 = solve(eq4, y22)[2]
ans_y22(y2 => ans_y2) |> simplify |> println

   -3*r2^2 + y2^2 - 2*y2*y22 + y22^2 から y22 を求める
   2*sqrt(r1)*sqrt(r2) + r1 + sqrt(3)*r2

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (14, 7) .// 2
   r3 = 2*(-2√(2(r1^3*r2^3)) + r1*r2*(r1 + r2))/(r1^2 - 6*r1*r2 + r2^2)
   y3 = r1 + √(2r1*r3 + r3^2)
   y2 = r1 + 2√(r1*r2)
   y22 = 2√(r1*r2) + r1 + √3r2
   plot()
   circle2(r1, r1, r1)
   circle2(r2, y2, r2, :blue)
   circle(0, y22, r2, :blue)
   circle(0, y3, r3, :green)
   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(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta/2)
       point(r2, y2, "小円:r2,(r2,y2)", :blue, :center, delta=-delta/2)
       point(0, y22, "小円:r2,(0,y22)", :blue, :center, delta=-delta/2)
       point(0, y3, "内円:r3,(0,y3)", :black, :center, delta=-delta/2)
       segment(-2r1, 0, 2r1, 0, :magenta, lw=0.5)
   end
end;

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

算額(その919)

2024年05月06日 | Julia

算額(その919)

一一〇 久喜市菖蒲町小林 小林神社 大正5年(1916)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

全円内に楕円と小円をそれぞれ 4 個入れる。楕円の長径,短径がそれぞれ 3 寸,2 寸のとき,全円の直径はいかほどか。

図形を 45°時計方向に回転させ,楕円の中心が x, y 軸上になるようにして以下の連立方程式を解く。
全円の半径と中心座標を R, (0, 0)
楕円の短半径,長半径と中心座標を a, b, (R - a, 0)
小円の半径と中心座標を r, (x, y); y = x
隣り合う 2 個の楕円の接点座標を (x0, y0); y0 = x0
として,以下の連立方程式を解くと,全円の半径 R を求めることができる。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive,
     a::positive, b::positive, x0, y0, x1, y1
y0 =  x0
eq1 = (x0 - R + a)^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*(x0 - R + a)/(a^2*y0) - 1
(R, x0) = solve([eq1, eq2], (R, x0))[2]

   (a + sqrt(a^2 + b^2), b^2/sqrt(a^2 + b^2))

全円の直径は 2(a + sqrt(a^2 + b^2)) である。
楕円の長径,短径がそれぞれ 3 寸,2 寸のとき,a = 1.5, b = 1 なので,6.60555127546399 寸である。

a = 1.5
b = 1
2(a + sqrt(a^2 + b^2))

   6.60555127546399

算額の「問」の答えとしてはここまでである。

---

以下では,図を描くために,小円の半径と中心座標の数値解を求める。

using SymPy
@syms R::positive, r::positive, x::positive,
     a::positive, b::positive, x0, y0, x1, y1
@syms R, r, x, a, b, x0, y0, x1, y1
R = a + sqrt(a^2 + b^2)
y = x
eq3 = (x1 - x)^2 + (y1 - y)^2 - r^2
eq4 = (x1 - R + a)^2/a^2 + y1^2/b^2 - 1
eq5 = -b^2*(x1 - R + a)/(a^2*y1) - (x - x1)/(y - y1)
eq6 = 2x^2 - (R - r)^2;

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)
   (r, x, x1, y1) = u
   return [
       -r^2 + (-x + x1)^2 + (-x + y1)^2,  # eq3
       -1 + y1^2/b^2 + (x1 - sqrt(a^2 + b^2))^2/a^2,  # eq4
       -(x - x1)/(x - y1) - b^2*(x1 - sqrt(a^2 + b^2))/(a^2*y1),  # eq5
       2*x^2 - (a - r + sqrt(a^2 + b^2))^2,  # eq6
   ]
end;

(a, b) = (3, 2) .// 2
iniv = BigFloat[0.78, 1.78, 2, 1]
res = nls(H, ini=iniv)

   ([0.7824436199506534, 1.7821438606147606, 1.7711417205363233, 0.999777596442297], true)

function draw(more=false)
   pyplot(size=(600, 600), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 3/2
   b = 2/2
   R = a + sqrt(a^2 + b^2)
   r = 0.7824436199506534
   @printf("楕円の長径が %g 寸,短径が %g 寸のとき,外円の直径は %g,等円の直径は %g である。\n", 2a, 2b, 2R, 2r)
   @printf("a = %g;  b = %g;  R = %g;  r = %g\n", a, b, R, r)
   plot()
   circle(0, 0, R)
   circle42(0, R - r, r, :blue)
   for theta = 45:90:315
       ox = (R - a)cosd(theta)
       oy = (R - a)sind(theta)
       ellipse(ox, oy, a, b, φ=theta, color=:green)
   end
   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, R, " R", :green, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :green, :left, :bottom, delta=delta/2)
       point(R - r, 0, "等円:r,(R-r,0)", :green, :center, delta=-delta/2)
   end
end;

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

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

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