裏 RjpWiki

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

算額(その816)

2024年03月28日 | Julia

算額(その816)

文化三年丙寅正月 藤田貞資門人 石州津和野 二介泉尹
藤田貞資(1807):続神壁算法

http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

直角三角形の田がある。直角を挟む二辺の短い方(鈎),長い方(股)が 6 間,8 間のとき,鈎,股に平行な直線で区切ってできる長方形(長辺(長),短辺(平)で,頂点の一つは斜辺上にある)の田の面積が最大になるとき,長,平はいかほどか。

鈎,股を A, B,鈎,股上にある点を a(平), b(長) とおくと,長方形の面積 S は S = a*b である。
以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms A, B, a, b, S

eq1 = (A - a)/b - A/B
eq2 = a/(B - b) - A/B
eq3 = a*b - S
res = solve([eq1, eq2, eq3], (a, b, S))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (A*(B - b)/B, b, A*b*(B - b)/B)

S = res[1][3]
S |> println

   A*b*(B - b)/B

問のように,A = 6, B = 8 のとき,b の取り方で,面積 S は b の関数で,S(b) は以下のように変化する。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(S(A => 6, B => 8), xlims=(0, 8), xlabel="b", ylabel="S(b)")

面積が最大になるの上図の曲線の接線の傾きが 0 になるときである。数値的に解くには,S(b) を b で微分し,S'(b) = 0 になるときの b の値を求める。

diff(S, b) |> println

   -A*b/B + A*(B - b)/B

solve(diff(S, b), b)[1] |> println

   B/2

b = B/2,a = A/2 が解である。
すなわち,長方形の短辺が 3 間,長辺が 6 間である。

function draw(more)
    pyplot(size=(300, 300), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (A, B) = (6, 8)
   (a, b) = (3, 4)
   plot([0, B, 0, 0], [0, 0, A, 0])
   rect(0, 0, b, a, :red)
   if more == true
       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)
       plot!(showaxis=false)
       point(0, A, " A(鈎)", :blue, :left, :bottom, delta=delta/2)
       point(B, 0, " B(股)", :blue, :left, :bottom, delta=delta/2)
       point(0, a, " a(平)", :red, :left, :bottom, delta=delta/2)
       point(b, 0, " b(長)", :red, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その815)

2024年03月27日 | Julia

算額(その815)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

大円内に 5 個の等円を描く。外円と等円の直径をそれぞれ 16 寸,10 寸としたとき,5 個の円の共通部分(黒積)を求めよ。

求めたい面積 S は
S1 = B を中心とする半径 r1 の円の ∠CBE/360 倍
S2 = 三角形 BCE
S3 = 三角形 OCE
とすると,
S = S1 - S2 + S3

まず,2個の等円の交点座標 C:(x0, y0) を求める。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms r0, r1, R, x0, y0, OX1, OY1, OX2, OY2
R = r0 - r1
(OX1, OY1) = (R*cosd(Sym(18)), R*sind(Sym(18)))
(OX2, OY2) = (0, R)
eq1 = (x0 - OX1)^2 + (y0 - OY1)^2 - r1^2
eq2 = (x0 - OX2)^2 + (y0 - OY2)^2 - r1^2
res = solve([eq1, eq2], (x0, y0));

x0 = res[2][1] |> factor
y0 = res[2][2] |> factor;

角度 ∠COE を求める

θ = asind(-x0/r1);

S1, S2, S3 を求める。式は長いが,そのままにしておく。

S1 = π*r1^2*θ/360
S2 = (-x0)*(OY2 - y0)/2
S3 = x0*y0/2
S1 |> println
S2 |> println
S3 |> println

   -r1^2*asin((-5*sqrt(2) + sqrt(10))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(8*r1*(-5 + sqrt(5))*sqrt(sqrt(5) + 5)))/2
   -(-5*sqrt(2) + sqrt(10))*(r0 - r1 + (sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(4*(-5 + sqrt(5))))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(16*(-5 + sqrt(5))*sqrt(sqrt(5) + 5))
   -(-5*sqrt(2) + sqrt(10))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))^2/(64*(-5 + sqrt(5))^2*sqrt(sqrt(5) + 5))

面積を求める。

S = 10(S1 - S2 + S3)
S |> println

   -5*r1^2*asin((-5*sqrt(2) + sqrt(10))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(8*r1*(-5 + sqrt(5))*sqrt(sqrt(5) + 5))) + 5*(-5*sqrt(2) + sqrt(10))*(r0 - r1 + (sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(4*(-5 + sqrt(5))))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(8*(-5 + sqrt(5))*sqrt(sqrt(5) + 5)) - 5*(-5*sqrt(2) + sqrt(10))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))^2/(32*(-5 + sqrt(5))^2*sqrt(sqrt(5) + 5))

r0 = 16//2, r1 = 10//2 を代入して数値を求める。
面積は 13.6341847764931 平方寸(歩)である。

S(r0 => 16//2, r1 => 10//2).evalf() |> println

   13.6341847764931

S(r0 => 1872//2, r1 => 1466//2).evalf() |> println

   915183.000045518

function transform(x, y; deg=0, dx=0, dy=0)
  return [cosd(deg) -sind(deg); sind(deg) cosd(deg)] * [x, y]
end;

function rotatefigure(x2, y2)
   for i in 1:5
       plot!(x2, y2, seriestype=:shape, color=:gray70, fillcolor=:gray70, lw=0)
       i == 5 && return
       (x2, y2) = transform(x2, y2, deg = 72)
   end
end;

function fillblack(r1, θ, OX2, OY2)
   x2 = []
   y2 = []
   for i in 270 - θ:0.5:270
       append!(x2, r1*cosd(i) + OX2)
       append!(y2, r1*sind(i) + OY2)
   end
   x2 = vcat(x2, -reverse(x2), [0, x2[1]])
   y2 = vcat(y2,  reverse(y2), [0, y2[1]])
   rotatefigure(x2, y2)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x = Vector{Float64}(undef, 6)
   y = Vector{Float64}(undef, 6)
   (r0, r1) = (16, 10) .// 2
   (r0, r1) = (1872, 1466) .// 2
   R = r0 - r1
   for i = 1:6
       θ = 18 + (i - 1)*72
       x[i] = R*cosd(θ)
       y[i] = R*sind(θ)
   end
   (OX1, OY1) = (R*cosd(18), R*sind(18))
   (OX2, OY2) = (0, R)
   x0 = (-5*sqrt(2) + sqrt(10))*(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(8*(-5 + sqrt(5))*sqrt(sqrt(5) + 5))
   y0 = -(sqrt(5)*r0 + 5*r0 - 5*r1 - sqrt(5)*r1 - sqrt(10)*sqrt(-5*r0^2 + sqrt(5)*r0^2 - 2*sqrt(5)*r0*r1 + 10*r0*r1 + sqrt(5)*r1^2 + 3*r1^2))/(4*(-5 + sqrt(5)))
   θ = asind(-x0/r1)
   println("θ = $θ")
   S = 10(π*r1^2*θ/360 - (-x0)*(OY2 - y0)/2 + x0*y0/2)
   println("S = $S")
   plot()
   fillblack(r1, θ, OX2, OY2)
   segment(OX2, OY2, x0, y0, :black, lw=1)
   segment(OX2, OY2, 0, OY2 - r1, :black, lw=1.5)
   segment(0, y0, x0, y0, :black, lw=1)
   segment(0, 0, x0, y0, :black, lw=1)
   circle(0, 0, r0)
   circle(0, 0, r0 - r1, :magenta)
   color = [:blue, :green, :gray80, :gray80, :gray80]
   for i = 1:5
       # segment(x[i], y[i], x[i + 1], y[i + 1], :blue)
       circle(x[i], y[i], r1, color[i])
   end
   if more == true
       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(OX1, OY1, " A:r1,(OX1,OY1)", :blue, :left, :bottom, delta=delta/2)
       point(OX2, OY2, " B:r1,(OX2,OY2)", :green, :left, :bottom, delta=delta/2)
       point(x0, y0, "C:(x0,y0) ", :black, :right, :vcenter)
       point(0, OY2 - r1, " D:(0,OY2 - r1) ", :black, :left, delta=-delta/2)
       point(0, y0, " E:(0,y0) ", :black, :left, :bottom, delta=delta/2)
       point(0, 0, " O", :black, :left, :vcenter)
       point(r0 - r1, 0, " r0-r1", :magenta, :left, :bottom, delta=delta/2)
       point(0, r0, " r0", :red, :left, :bottom)
       # plot!(xlims=(-3, 5), ylims=(-3, 3))
   end
end;

   θ = 26.63149441418678
   S = 915183.0000455183

r0,r1 を与えて黒積を返す関数を定義する。

function func(r0, r1)
   p = √5 + 5
   q = √5 - 5
   s = p*(r0 - r1) - √10*sqrt(q*r0*(r0 - 2r1) + r1^2*(√5 + 3))
   t = (√10 - 5√2)*s/(8q√p)
   5(r0*t - r1*t - r1^2*asin(t/r1))
end;
func(8, 5)

   13.634184776493093

この関数を使って,黒積が整数に極めて近い値になるときの r0, r1 の組を探索することもできる。
そのような候補として,r0 = 1872/2, r1 = 1466/2 のとき,黒積が 915183.0000455179 になる

func(936, 733)

   915183.0000455179

 

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

算額・和算の挑戦状

2024年03月27日 | Julia

精要算法(下巻)にもある問題で,算法助術にも助術用例3として解説されているものである。

そこでは,等円が5個の場合を取り上げている。

これを一般化して,外円内に n  個の等円を入れる。外円と等円に挟まれる部分の面積(黒積)が与えられたとき,等円の直径を求めよ。

n = 9, 黒積 = 66522.03 のとき,等円の直径はいかほどか...

 

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

算額(その814)

2024年03月27日 | Julia

算額(その814)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

外円内に 5 個の等円を入れる。外円と等円に挟まれる部分の面積(黒積)が 2330.89 歩のとき,等円の直径を求めよ。

黒積は,「扇形 OCD - 直角三角形 OAB - 扇形ACB」の面積の 10 倍である。
等円の半径 OB = R,AB = r とすれば,∠AOB = 36° より,R = r/sind(36)
OC = OB + BC = R + r = Rr とすれば,
扇形 OCD = πRr^2/10
直角三角形 OAB = r*R*cosd(36)/2
扇形 ACB = πr^2*(126/360)
よって,黒積 = 10*(π*Rr^2/10 - r*R*cosd(36)/2 - π*r^2*126/360)
黒積 = 2332.89 を解いて r を求めると,r = 21.5000025771553,等円の直径は 43.0000051543107 寸である

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms r::positive, R::positive, Rr::positive,
     黒積::positive, BLK::positive

R = solve(R*sind(Sym(36)) - r, R)[1]
Rr = R + r
R |> println
Rr |> println
(R.evalf(), Rr.evalf()) |> println

   2*sqrt(2)*r/sqrt(5 - sqrt(5))
   r + 2*sqrt(2)*r/sqrt(5 - sqrt(5))
   (1.70130161670408*r, 2.70130161670408*r)

黒積 = PI*Rr^2/10 - R*r*cosd(Sym(36))/2 - PI*r^2*126//360
eq = 10黒積 - BLK  # 2332.89;
eq |> println

   -BLK - 7*pi*r^2/2 - 10*sqrt(2)*r^2*(1/4 + sqrt(5)/4)/sqrt(5 - sqrt(5)) + pi*(r + 2*sqrt(2)*r/sqrt(5 - sqrt(5)))^2

方程式を解いて,等円の直径を求める。

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

   2*sqrt(2)*sqrt(BLK)*(5 - sqrt(5))^(3/4)/sqrt(-8*sqrt(10)*pi - 20*sqrt(10) - 9*pi*sqrt(5 - sqrt(5)) + 5*sqrt(5)*pi*sqrt(5 - sqrt(5)) + 40*sqrt(2)*pi)

黒積が 2332.89 のときの等円の直径を求める。

res(BLK => 2332.89) |> println

   96.6*sqrt(2)*(5 - sqrt(5))^(3/4)/sqrt(-8*sqrt(10)*pi - 20*sqrt(10) - 9*pi*sqrt(5 - sqrt(5)) + 5*sqrt(5)*pi*sqrt(5 - sqrt(5)) + 40*sqrt(2)*pi)

等円の直径は 43 寸である。

res(BLK => 2332.89).evalf() |> println

   43.0000051543107

function transform(x, y; deg=0, dx=0, dy=0)
  return [cosd(deg) -sind(deg); sind(deg) cosd(deg)] * [x, y]
end;

function rotatefigure(x2, y2)
   for i in 1:5
       plot!(x2, y2, color=:gray80, seriestype=:shape, fillcolor=:gray80)
       i == 5 && return
       (x2, y2) = transform(x2, y2, deg = 72)
   end
end;

function fillblack(r, Rr, x, y)
   x2 = []
   y2 = []
   for i in 180:0.5:306
       append!(x2, r*cosd(i) + x[1])
       append!(y2, r*sind(i) + y[1])
   end
   for i in 306:-0.5:270
       append!(x2, Rr*cosd(i))
       append!(y2, Rr*sind(i))
   end
   x2 = vcat(x2, -reverse(x2))
   y2 = vcat(y2, reverse(y2))
   rotatefigure(x2, y2)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x = Vector{Float64}(undef, 6)
   y = Vector{Float64}(undef, 6)
   r = 21.5000025771553
   R = 2√2*r/sqrt(5 - √5)
   Rr = R + r
   blk = 10(π*Rr^2/10 - R*r*cosd(36)/2 - π*r^2*126/360)
   @printf("等円の直径 = %g;  黒積 = %g;  R = %g;  Rr = %g\n", 2r, blk, R, Rr)
   for i = 1:6
       θ = 306 + (i - 1)*72
       x[i] = R*cosd(θ)
       y[i] = R*sind(θ)
   end
   plot()
   fillblack(r, Rr, x, y)
   for i = 1:5
       segment(x[i], y[i], x[i + 1], y[i + 1], :blue)
       circle(x[i], y[i], r, :green)
   end
   circle(0, 0, R)
   circle(0, 0, Rr, :magenta)
   segment(0, 0, 0, -Rr, :black, lw=1)
   segment(0, 0, Rr*cosd(306), Rr*sind(306), :black, lw=1)
   if more == true
       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(x[1], y[1])
       point(R, 0, " R", :red, :left, :bottom, delta=delta)
       point(Rr, 0, " Rr", :magenta, :left, :bottom, delta=delta)
       point(0, -R*cosd(36), "A ", :black, :right, delta=-delta/2)
       point(x[1], y[1], "  B", :black, :left, delta=-delta/2)
       point(Rr*cosd(306), Rr*sind(306), "C", :black, :center, delta=-delta/2)
       point(0, -Rr, "D ", :black, :right, :bottom, delta=delta/2)
       point(0, 0, "O", :black, :center, :bottom, delta=delta/2)
   end
end;

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

算額(その813)

2024年03月26日 | Julia

算額(その813)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

不等辺三角形内に全円,大円,中円,小円を入れる。大円,中円,小円の直径がそれぞれ 256 寸,225 寸,144 寸のとき,全円の直径はいかほどか。

底辺の長さ a,頂点の座標を (b, h)
全円の半径と中心座標を r0, (x0, r0)
大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (x3, y3)
とおき,以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b::positive, h::positive,
     r0::positive, x0::positive,
     r1::positive, x1::positive,
     r2::positive, x2::positive,
     r3::positive, x3::positive, y3::positive,
     d
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = (x3 - x2)^2 + (y3 - r2)^2 - (r2 + r3)^2
eq4 = (a + sqrt(b^2 + h^2) + sqrt((a - b)^2 + h^2))r0 - a*h
eq5 = r0/(a - x0) - r1/(a - x1)
eq6 = r0/x0 - r2/x2
eq7 = dist(0, 0, b, h, x2, r2) - r2^2
eq8 = dist(0, 0, b, h, x3, y3) - r3^2
eq9 = dist(a, 0, b, h, x3, y3) - r3^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)
   (a, b, h, x0, r0, x1, x2, x3, y3) = u
   return [
       (r1 - r2)^2 - (r1 + r2)^2 + (x1 - x2)^2,  # eq1
       (-r1 + y3)^2 - (r1 + r3)^2 + (x1 - x3)^2,  # eq2
       (-r2 + y3)^2 - (r2 + r3)^2 + (-x2 + x3)^2,  # eq3
       -a*h + r0*(a + sqrt(b^2 + h^2) + sqrt(h^2 + (a - b)^2)),  # eq4
       r0/(a - x0) - r1/(a - x1),  # eq5
       r0/x0 - r2/x2,  # eq6
       -r2^2 + (-b*(b*x2 + h*r2)/(b^2 + h^2) + x2)^2 + (-h*(b*x2 + h*r2)/(b^2 + h^2) + r2)^2,  # eq7
       -r3^2 + (-b*(b*x3 + h*y3)/(b^2 + h^2) + x3)^2 + (-h*(b*x3 + h*y3)/(b^2 + h^2) + y3)^2,  # eq8
       -r3^2 + (-h*(h*y3 + (-a + b)*(-a + x3))/(h^2 + (-a + b)^2) + y3)^2 + (-a + x3 - (-a + b)*(h*y3 + (-a + b)*(-a + x3))/(h^2 + (-a + b)^2))^2,  # eq9
   ]
end;

(r1, r2, r3) = (256, 225, 144) .// 2
iniv = BigFloat[1000, 355, 358, 385, 162, 513, 269, 368, 269]
res = nls(H, ini=iniv)
res |> println
2res[1][5] |> println

   ([1014.0, 354.88757396449705, 357.8698224852071, 384.0, 160.0, 510.0, 270.0, 367.9881656804734, 268.8284023668639], true)
   320.0

大円,中円,小円の直径がそれぞれ 256 寸,225 寸,144 寸のとき,全円の直径は 320 寸である。

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

a = 1014;  b = 354.888;  h = 357.87;  x0 = 384;  r0 = 160;  x1 = 510;  x2 = 270;  x3 = 367.988;  y3 = 268.828

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (256, 225, 144) .// 2
   (a, b, h, x0, r0, x1, x2, x3, y3) = res[1]
   @printf("全円の直径 = %g\n", 2r0)
   @printf("a = %g;  b = %g;  h = %g;  x0 = %g;  r0 = %g;  x1 = %g;  x2 = %g;  x3 = %g;  y3 = %g\n", a, b, h, x0, r0, x1, x2, x3, y3)
   plot([0, a, b, 0], [0, 0, h, 0], color=:black, lw=0.5)
   circle(x0, r0, r0, :orange)
   circle(x1, r1, r1, :blue)
   circle(x2, r2, r2, :green)
   circle(x3, y3, r3)
   if more == true
       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(x0, r0, " 全円:r0,(x0,r0)", :black, :left, :vcenter)
       point(x1, r1, "大円:r1\n(x1,r1)", :blue, :center, delta=-delta)
       point(x2, r2, "中円:r2\n(x2,r2)", :green, :center, delta=-delta)
       point(x3, y3, " 小円:r3,(x3,y3)", :red, :left, :vcenter)
       point(a, 0, " a", :black, :left, :bottom, delta=delta)
       point(b, h, "(b,h)", :black, :left, :bottom, delta=delta)
   end
end;

 

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

算額(その812)

2024年03月26日 | Julia

算額(その812)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

長方形内に大円,中円,小円が入っている。長方形の長辺,短辺がそれぞれ 86634 寸,77008 寸のとき,大円の直径はいかほどか。

長方形の長辺,短辺を a, b
大円の半径と中心座標を r1, (a - r1, r1)
中円の半径と中心座標を r2, (r2, b - r2)
小円の半径と中心座標を r3, (r3, r3)
とおいて以下の連立方程式を解く。
なお,SymPy の能力的に,a, b を変数のままにして解くことはできない。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, r2::positive, r3::positive,
     a::positive, b::positive
(a, b) = (86634, 77008)
eq1 = (a - r1 - r2)^2 + (b - r2 - r1)^2 - (r1 + r2)^2
eq2 = (a - r1 - r3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (r2 - r3)^2 + (b - r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (r1, r2, r3))
res[2]

   (-57756*sqrt(338 + 240*sqrt(2)) + 24065 + 81821*sqrt(169 + 120*sqrt(2)), -81821*sqrt(169 + 120*sqrt(2)) + 24065 + 57756*sqrt(338 + 240*sqrt(2)), -4813*sqrt(169 + 120*sqrt(2)) + 24065 + 57756*sqrt(2))

6 通りの解が得られるが,2 番目のものが適解である。

長方形の長辺,短辺がそれぞれ 86634 寸,77008 寸のとき,大円の直径は 53345.0001130072 寸である。

2res[2][1].evalf() |> println

   53345.0001130072

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

r1 = 26672.5;  r2 = 21457.5;  r3 = 17166.1

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (86634, 77008)
   t = sqrt(169 + 120√2)
   u = 81821 - 57756√2
   (r1, r2, r3) = 24065 .+ (t*u, -t*u, 57756√2 - 4813t)
   @printf("大円の直径 = %g;  中円の直径 = %g;  小円の直径 = %g\n", 2r1, 2r2, 2r3)
   @printf("r1 = %g;  r2 = %g;  r3 = %g\n", r1, r2, r3)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:black, lw=0.5)
   circle(a - r1, r1, r1)
   circle(r2, b - r2, r2, :blue)
   circle(r3, r3, r3, :green)
   if more == true
       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 - r1, r1, "大円:r1,(a-r1,r1)", :red, :center, delta=-delta)
       point(r1, b - r2, "中円:r2,(r2,b-r2)", :blue, :center, delta=-delta)
       point(r3, r3, "小円:r3,(r3,r3)", :green, :center, delta=-delta)
       point(a, 0, " a", :black, :left, :bottom, delta=delta)
       point(0, b, " b", :black, :left, :bottom, delta=delta)
   end
end;

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

算額(その811)

2024年03月25日 | Julia

算額(その811)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

長方形に楕円と円が入っている。長方形の長辺と短辺が 16 寸,9 寸のとき,円の直径はいかほどか。

算法助術の公式90 そのものである。
楕円の長径,短径を a,b,円の直径を d とすれば,
d^2 +ab -2d√ab - 2d(a + b) = 0
が成り立つ。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b::positive, d::positive
formula90 = d^2 + a*b - 2d*sqrt(a*b) - 2d*(a + b);

d について解けば,a, b を与えて d を求める公式になる。

d = solve(formula90, d)[1]
d |> println

   sqrt(a)*sqrt(b) + a + b - sqrt(2*a^(3/2)*sqrt(b) + 2*sqrt(a)*b^(3/2) + a^2 + 2*a*b + b^2)

楕円の長径,短径を与えれば円の直径が求まる。

d(a => 16, b => 9) |> println

   2

関数にしてみる。

func(a, b) = sqrt(a*b) + a + b - sqrt(2a^(3/2)*√b + 2√a*b^(3/2) + a^2 + 2a*b + b^2);

長径,短径が 16, 9 のとき,円の直径は 2 である。

func(16, 9)

   2.0

この公式,関数を用いるとき,長径,短径をあたえれば直径が得られ,長半径,短半径を与えれば半径が得られる。

func(a, b) = sqrt(a*b) + a + b - sqrt(2a^(3/2)*√b + 2√a*b^(3/2) + a^2 + 2a*b + b^2);

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (16, 9) .// 2
   r = func(a, b)
   @printf("円の直径 = %g\n", 2r)
   # 図形描画システムと統一を図るため長半径,短半径,半径にする
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:blue, lw=0.5)
   ellipse(0, 0, a, b, color=:red)
   circle4(a - r, b - r, r, :green)
   if more == true
       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 - r, b - r, "(a-r,b-r) ", :black, :right, :bottom, delta=delta)
       point(a, 0, " a", :red, :left, :bottom, delta=delta)
       point(0, b, " b", :red, :left, :bottom, delta=delta)        
   end
end;

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

算額(その810)

2024年03月25日 | Julia

算額(その810)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

直角三角形内に,大円,中円,小円を入れる。中円,小円の直径がそれぞれ 4635 寸,2060 寸のとき,大円の直径はいかほどか。

直角三角形の中に3個の円を入れる類題は多い。算額(その740),算額(その453),算額(その188

大円の半径と中心座標を r1, (r1, r2)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (r3, y3)
直角を挟む短い方の辺と長いほうの辺の長さを a, b
とおき,以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, r2::positive, x2::positive,
     r3::positive, y3::positive, a::positive,
     b::positive
eq1 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = r1/(b - r1) - r2/(b - x2)
eq4 = r1/(a - r1) - r3/(a - y3)
eq5 = a + b - sqrt(a^2 + b^2) - 2r1;
#res = solve([eq1, eq2, eq3, eq4, eq5], (r1, x2, y3, a, b))

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)
   (r1, x2, y3, a, b) = u
   return [
       (-r1 + x2)^2 + (r1 - r2)^2 - (r1 + r2)^2,  # eq1
       (-r1 + y3)^2 + (r1 - r3)^2 - (r1 + r3)^2,  # eq2
       r1/(b - r1) - r2/(b - x2),  # eq3
       r1/(a - r1) - r3/(a - y3),  # eq4
       a + b - 2*r1 - sqrt(a^2 + b^2),  # eq5
   ]
end;

(r2, r3) = (4635, 2060) .// 2
iniv = BigFloat[3500, 9200, 7300, 8900, 20300]
res = nls(H, ini=iniv)

   ([3515.5003002014273, 9224.150559489277, 7321.26713972666, 8898.390201990349, 20267.383999142992], true)

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

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

r1 = 3515.5;  x2 = 9224.15;  y3 = 7321.27;  a = 8898.39;  b = 20267.4

function draw(more)
    pyplot(size=(700, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (4635, 2060) .// 2
   (r1, x2, y3, a, b) = res[1]
   @printf("大円の直径 = %g\n", 2r1)
   @printf("r1 = %g;  x2 = %g;  y3 = %g;  a = %g;  b = %g\n", r1, x2, y3, a, b)
   plot([0, b, 0, 0], [0, 0, a, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, r2, r2, :green)
   circle(r3, y3, r3, :magenta)
   if more == true
       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\n(r1,r1)", :red, :center, :bottom, delta=3delta)
       point(x2, r2, " 中円:r2\n(x2,r2)", :green, :center, delta=-delta)
       point(r3, y3, "小円:r3,(r3,y3) ", :blue, :left, delta=-3delta)
       point(0, a, " a", :black, :left, :bottom)
       point(b, 0, " b", :black, :left, :bottom)
   end
end;

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

算額(その809)

2024年03月25日 | Julia

算額(その809)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

不等辺三角形内に全円,大円,中円,小円を入れる。大円,中円,小円の直径がそれぞれ 9 寸,4 寸,1 寸のとき,全円の直径はいかほどか。

三角形の底辺の長さを a,頂点の座標を (x0, y0) 
全円の半径と中心座標を r1, (x1, r1)
大円の半径と中心座標を r2, (x2, r2)
中円の半径と中心座標を r3, (x3, r3)
小円の半径と中心座標を r4, (x4, y4)
とおき,以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, x1::positive, r2::positive,
     x2::positive, r3::positive, x3::positive,
     r4::positive, x4::positive, y4::positive,
     a::positive, x0::positive, y0::positive,
     d
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x1 - x4)^2 + (r1 - y4)^2 - (r1 + r4)^2
# eq4 = r1/(a - x1) - r3/(a - x3)
eq5 = r2/x2 - r1/x1
eq6 = (a + sqrt(x0^2 + y0^2) + sqrt((a - x0)^2 + y0^2))*r1 - y0*a
eq7 = dist(0, 0, x0, y0, x2, r2) - r2^2
eq8 = dist(0, 0, x0, y0, x4, y4) - r4^2
eq9 = dist(a, 0, x0, y0, x3, r3) - r3^2
eq10 = dist(a, 0, x0, y0, x4, y4) - r4^2
eq7 = numerator(apart(eq7, d));
eq8 = numerator(apart(eq8, d));
eq9 = numerator(apart(eq9, d));
eq10 = numerator(apart(eq10, d));
# res = solve([eq1, eq2, eq3, eq5, eq6, eq7, eq8, eq9, eq10],
#     (r1, x1, x2, x3, x4, y4, a, x0, y0))

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)
   (r1, x1, x2, x3, x4, y4, a, x0, y0) = u
   return [
       (r1 - r2)^2 - (r1 + r2)^2 + (x1 - x2)^2,  # eq1
       (r1 - r3)^2 - (r1 + r3)^2 + (x1 - x3)^2,  # eq2
       -(r1 + r4)^2 + (r1 - y4)^2 + (x1 - x4)^2,  # eq3
       -r1/x1 + r2/x2,  # eq5
       -a*y0 + r1*(a + sqrt(x0^2 + y0^2) + sqrt(y0^2 + (a - x0)^2)),  # eq6
       -r2^2*y0^2 - 2*r2*x0*x2*y0 + x2^2*y0^2,  # eq7
       -r4^2*x0^2 - r4^2*y0^2 + x0^2*y4^2 - 2*x0*x4*y0*y4 + x4^2*y0^2,  # eq8
       -2*a^2*r3*y0 + a^2*y0^2 + 2*a*r3*x0*y0 + 2*a*r3*x3*y0 - 2*a*x3*y0^2 - r3^2*y0^2 - 2*r3*x0*x3*y0 + x3^2*y0^2,  # eq9
       -a^2*r4^2 + a^2*y0^2 - 2*a^2*y0*y4 + a^2*y4^2 + 2*a*r4^2*x0 + 2*a*x0*y0*y4 - 2*a*x0*y4^2 - 2*a*x4*y0^2 + 2*a*x4*y0*y4 - r4^2*x0^2 - r4^2*y0^2 + x0^2*y4^2 - 2*x0*x4*y0*y4 + x4^2*y0^2,  # eq10
   ]
end;

(r2, r3, r4) = (9, 4, 1) .// 2
iniv = BigFloat[5.5, 54.7, 44.8, 61.4, 57, 11.1, 65.1, 57.21, 11.6]
res = nls(H, ini=iniv)

   ([5.5, 54.7243090408641, 44.7744346697979, 61.357558621574896, 56.97961389830577, 11.06, 65.14798695340964, 57.20514438404994, 11.616], true)

全円の直径は 11 寸である。

算額にはよくあることであるが,実際の図と算額に描かれている図には大きな違いがある。

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

r1 = 5.5;  x1 = 54.7243;  x2 = 44.7744;  x3 = 61.3576;  x4 = 56.9796;  y4 = 11.06;  a = 65.148;  x0 = 57.2051;  y0 = 11.616

function draw(more)
    pyplot(size=(800, 170), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, r4) = (9, 4, 1) .// 2
   (r1, x1, x2, x3, x4, y4, a, x0, y0) = res[1]
   @printf("全円の直径 = %g\n", 2r1)
   @printf("r1 = %g;  x1 = %g;  x2 = %g;  x3 = %g;  x4 = %g;  y4 = %g;  a = %g;  x0 = %g;  y0 = %g\n", r1, x1, x2, x3, x4, y4, a, x0, y0)
   plot([0, a, x0, 0], [0, 0, y0, 0], color=:black, lw=0.5)
   circle(x1, r1, r1)
   circle(x2, r2, r2, :green)
   circle(x3, r3, r3, :magenta)
   circle(x4, y4, r4, :blue)
   if more == true
       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(x1, r1, "全円:r1\n(x1,r1)", :red, :center, :bottom, delta=3delta)
       point(x2, r2, "大円:r2\n(x2,r2)", :green, :center, delta=-3delta)
       point(x3, r3, "中円:r3,(x3,r3) ", :black, :right, :vcenter)
       point(x4, y4, "小円:r4,(x4,y4) ", :blue, :right, :bottom, delta=2delta)
       point(x0, y0, " (x0,y0)", :black, :left, :vcenter)
       point(a, 0, " a", :black, :left, :bottom)
   end
end;

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

算額(その808)

2024年03月24日 | Julia

算額(その808)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

台形(半梯)の中に甲円と乙円が入っている。大頭,小頭,縦(台形の高さ)がそれぞれ 8.5 寸,4.9 寸,4.8 寸,甲円の直径が 4 寸のとき,乙円の直径はいかほどか。

大頭,小頭,縦(台形の高さ)を b, c, a
甲円の半径と中心座標を r1, (r1, r1)
乙円の半径と中心座標を r2, (a - r2, y2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b::positive, c::positive,
     r1::positive, r2::positive, y2::positive, d
(a, b, c, r1) = (48//10, 85//10, 49//10, 4//2)
eq1 = (a - r2 - r1)^2 + (y2 - r1)^2 - (r1 + r2)^2 |> expand
eq2 = distance(0, b, a, c, a - r2, y2) - r2^2
eq2 = numerator(apart(eq2, d))
res = solve([eq1, eq2], (r2, y2))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (1, 22/5)

大頭,小頭,縦(台形の高さ)がそれぞれ 8.5 寸,4.9 寸,4.8 寸,甲円の直径が 4 寸のとき,乙円の直径は 2 寸である。

実際の図は,算額の図とは大違いである。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, r1) = (48//10, 85//10, 49//10, 4//2)
   (r2, y2) = (1, 22/5)
   @printf("乙円の直径 = %.9g\n", 2r2)
   @printf("r2 = %g;  y2 = %g\n", r2, y2)
   plot([0, a, a, 0, 0], [0, 0, c, b, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(a - r2, y2, r2, :blue)
   if more == true
       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)
       point(a - r2, y2, "乙円:r2\n(a-r2,y2)", :blue, :center, delta=-delta)
       point(0, b, " b", :black, :left, :bottom)
       point(a, c, "(a,c) ", :black, :right, delta=-delta/2)
       point(a, 0, "a ", :black, :right, :bottom, delta=delta/2)
   end
end;

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

算額(その807)

2024年03月24日 | Julia

算額(その807)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

二等辺三角形内に隔線を入れ,3 個の等円を入れる。底辺が 120 寸,斜辺が 87 寸のとき,等円の直径はいかほどか。

二等辺三角形の底辺の長さを 2a, 高さを h
3 本の隔線の交点座標を (0, b)
等円の半径と中心座標を r, (0, r), (r, y)
とおき,以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b, h::positive, r::positive, y::positive, d
a = 120//2
eq1 = a^2 + h^2 - 87^2
eq2 = dist(0, h, a, 0, r, y) - r^2
eq3 = (b - r)*a/sqrt(a^2 + b^2)- r
eq4 = dist(0, b, a, 0, r, y) - r^2
res = solve([eq1, eq2, eq3, eq4], (r, y, b, h))

   2-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (12, 33, 25, 63)
    (180/7, 513/7, 63, 63)

2 組の解が得られるが,最初のものが適解である。

底辺が 120 寸,斜辺が 87 寸のとき,等円の直径は 24 寸である。

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

   r = 12;  y = 33;  b = 25;  h = 63

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 120//2
   (r, y, b, h) = (12, 33, 25, 63)
   @printf("等円の直径 = %.9g\n", 2r)
   @printf("r = %g;  y = %g;  b = %g;  h = %g\n", r, y, b, h)
   plot([a, 0, -a, a], [0, h, 0, 0], color=:black, lw=0.5)
   circle(0, r, r)
   circle2(r, y, r)
   segment(-a, 0, 0, b, :blue)
   segment(a, 0, 0, b, :blue)
   segment(0, b, 0, h, :blue)
   if more == true
       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, y, "等円:r\n(r,y)", :red, :center, delta=-delta)
       point(0, r, "等円:r\n(0,r)", :red, :center, delta=-delta)
       point(a, 0, " a", :blue, :center, :bottom, delta=delta)
       point(0, b, " b", :blue, :center, delta=-2delta)
       point(0, h, " h", :black, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その806)

2024年03月24日 | Julia

算額(その806)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

直角三角形内に,大,中,小の 3 円を入れる。直角を挟む二辺の短い方(鈎),長い方(股)の長さが 2352 寸,3689 寸のとき,小円の直径を求めよ。

鈎,股の長さを a, b
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (x3, r3)
とおき,以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive,
     a::positive, b::positive
(a, b) = (2352, 3689)
eq1 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x3 - r1)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq4 = a + b - sqrt(a^2 + b^2) - 2r1
eq4 = (a + b - 2r1)^2 - (a^2 + b^2)
eq5 = r1/(b - r1) - r2/(b - x2)
eq5 = r1*(b - x2) - r2*(b - r1)
res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, x2, r3, x3));
res[1]

   (833, 7497/16, 4165/2, 153, 1547)

5 組の解が得られるが,最初のものが適解である。

鈎,股が 2352 寸,3689 寸のとき,小円の直径は 306 寸である。

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

 r1 = 833;  r2 = 468.562;  x2 = 2082.5;  r3 = 153;  x3 = 1547

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x2, r3, x3) = (833, 7497/16, 4165/2, 153, 1547)
   @printf("小円の直径 = %.9g\n", 2r3)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  a = %gx3\n", r1, r2, x2, r3, x3)
   plot([0, b, 0, 0], [0, 0, a, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   if more == true
       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, a, " a", :black, :left, :bottom) 
       point(b, 0, " b", :black, :left, :bottom) 
       point(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta)
       point(x2, r2, "中円:r2,(x2,r2)", :blue, :center, delta=-delta)
       point(x3, r3, " 小円:r3,(x3,r3)", :black, :left, :vcenter)
   end
end;

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

算額(その805)

2024年03月24日 | Julia

算額(その805)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

台形内に甲円,乙円,丙円を入れる。大頭(下底)が 782 寸,小頭(上底)が 291 寸のとき,甲円の直径はいかほどか。

横倒しになっているが,台形の大頭,小頭,高さを b, c, a とおく。
甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (r2, y2)
丙円の半径と中心座標を r3, (r3, r3)
とおいて,以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, r2::positive, y2::positive, r3::positive, x3::positive,
     a::positive, b::positive, c::positive, d
eq1 = (a - r1 - r2)^2 + (y2 - r1)^2 - (r1 + r2)^2 |> expand
eq2 = (a - r1 - r3)^2 + (r1 - r3)^2 - (r1 + r3)^2 |> expand
eq3 = (r2 - r3)^2 + (y2 - r3)^2 - (r2 + r3)^2 |> expand
eq4 = dist(0, b, a, c, a - r1, r1) - r1^2
eq4 = numerator(apart(eq4, d))
eq5 = dist(0, b, a, c, r2, y2) - r2^2
eq5 = numerator(apart(eq5, d))

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)
   (r1, r2, y2, r3, a) = u
   return [
       a^2 - 2*a*r1 - 2*a*r2 + r1^2 - 2*r1*y2 + y2^2,  # eq1
       a^2 - 2*a*r1 - 2*a*r3 + r1^2 - 2*r1*r3 + r3^2,  # eq2
       -4*r2*r3 + r3^2 - 2*r3*y2 + y2^2,  # eq3
       a^2*c^2 - 2*a^2*c*r1 + 2*a*b*c*r1 - 2*a*b*r1^2 - 2*a*c^2*r1 + 2*a*c*r1^2,  # eq4
       a^2*b^2 - 2*a^2*b*y2 - a^2*r2^2 + a^2*y2^2 - 2*a*b^2*r2 + 2*a*b*c*r2 + 2*a*b*r2*y2 - 2*a*c*r2*y2,  # eq5
   ]
end;

(b, c) = (782, 391)
iniv = BigFloat[250, 200, 500, 175, 780]
res = nls(H, ini=iniv)

   ([242.00048108871312, 183.47778926651384, 484.00096217742623, 151.04933058240556, 775.4318754192511], true)

甲円の直径は 242 寸有奇である。

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

r1 = 242;  r2 = 183.478;  y2 = 484.001;  r3 = 151.049;  a = 775.432

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, y2, r3, a) = res[1]
   @printf("甲円の直径 = %.9g\n", 2r1)
   @printf("r1 = %g;  r2 = %g;  y2 = %g;  r3 = %g;  a = %g\n", r1, r2, y2, r3, a)
   plot([0, a, a, 0, 0], [0, 0, c, b, 0], color=:black, lw=0.5)
   circle(a - r1, r1, r1)
   circle(r2, y2, r2, :blue)
   circle(r3, r3, r3, :green)
   if more == true
       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 - r1, r1, "甲円:r1,(a-r1,r1)", :red, :center, delta=-delta)
       point(r2, y2, "乙円:r2,(r2,y2)", :blue, :center, delta=-delta)
       point(r3, r3, "丙円:r3,(r3,r3)", :green, :center, delta=-delta)
       point(0, b, " b", :black, :left, :bottom)
       point(a, 0, " a", :black, :left, :bottom)
       point(a, c, " (a,c)", :black, :left, :bottom)
       plot!(xlims=(0, 850))
   end
end;

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

算額(その804)

2024年03月24日 | Julia

算額(その804)

藤田貞資:精要算法(下巻),天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

長方形内に甲,乙,丙,丁の 4 円を入れる。甲円の直径が 5140 寸のとき,乙円の直径はいかほどか。

甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (r2, 2r1 - r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (r4, r4)
長方形の長辺を a とおく。短辺は 2r1 である。
以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, r1::positive, r2::positive, r3::positive, x3::positive, r4::positive
eq1 = (a - r1 - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (a - r1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x3 - r2)^2 + (2r1 - r2 - r3)^2 - (r2 + r3)^2
eq4 = (r2 - r4)^2 + (2r1 - r2 - r4)^2 - (r2 + r4)^2
eq5 = (x3 - r4)^2 + (r3 - r4)^2 - (r3 + r4)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (a, r2, r3, x3, r4));
res[3]

   (r1*(3*sqrt(28*sqrt(2) + 185) + 34*sqrt(2) + 8*sqrt(56*sqrt(2) + 370) + 187)/136, r1*(-43*sqrt(28*sqrt(2) + 185) - 187 + 10*sqrt(56*sqrt(2) + 370) + 272*sqrt(2))/(68*(-sqrt(28*sqrt(2) + 185) + 2*sqrt(2) + 7)), r1*(-19*sqrt(28*sqrt(2) + 185) - 68*sqrt(2) + 6*sqrt(56*sqrt(2) + 370) + 221)/(17*(-sqrt(28*sqrt(2) + 185) + 2*sqrt(2) + 7)), r1*(-58*sqrt(56*sqrt(2) + 370) - 391 + 25*sqrt(28*sqrt(2) + 185) + 612*sqrt(2))/(68*(-sqrt(28*sqrt(2) + 185) + 2*sqrt(2) + 7)), r1*(-sqrt(7*sqrt(2)/16 + 185/64) + sqrt(2)/4 + 15/8))

4 組の解が得られるが,3 番目のものが適解である。

乙円の直径は,t = √(28√2 + 185) とおくと,甲円の直径の ((7 - 4√2)t + 119 - 34√2)/136 倍である。

甲円の直径が 5140 寸のとき,乙円の直径は 3441.0001373070395 寸である。

「答」では「三千四百四十寸◯◯◯有奇」とあるが,誤記であろう。

t = √(28√2 + 185)
5140*((7 - 4√2)t + 119 - 34√2)/136

   3441.0001373070395

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

a = 8496.06;  r2 = 1720.5;  r3 = 959.736;  x3 = 2785.03;  r4 = 912.939

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 5140//2
   t = sqrt(28*sqrt(2) + 185)
   u = 2√2 + 7 - t
   (a, r2, r3, x3, r4) = r1 .* (
       (3*t + 8*√2t + 187 + 34√2)/136,
       (-43*t + 10*√2t - 187 + 272√2)/68u,
       (-19*t + 6*√2t + 221 - 68√2)/17u,
       (25*t -58*√2t - 391 + 612√2)/68u,
       (-sqrt(7√2/16 + 185/64) + √2/4 + 15/8))
   @printf("乙円の直径 = %g\n", 2r2)
   @printf("a = %g;  r2 = %g;  r3 = %g;  x3 = %g;  r4 = %g\n", a, r2, r3, x3, r4)
   plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   circle(a - r1, r1, r1)
   circle(r2, 2r1 - r2, r2, :blue)
   circle(x3, r3, r3, :green)
   circle(r4, r4, r4, :magenta)
   if more == true
       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 - r1, r1, "甲円:r1,(a-r1,r1)", :red, :center, delta=-delta)
       point(r2, 2r1 - r2, "乙円:r2,(r2,2r1-r2)", :blue, :center, delta=-delta)
       point(x3, r3, "丙円:x3\n(x3,r3)", :green, :center, delta=-delta)
       point(r4, r4, "丁円:r4\n(r4,r4)", :magenta, :center, delta=-delta)
       point(a, 0, " a", :black, :left, :bottom, delta=delta)
   end
end;

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

算額・和算ページのまとめ(2)

2024年03月24日 | ブログラミング

神壁算法

  1. 算額(その771)
    https://blog.goo.ne.jp/r-de-r/e/d86f8eafad5635b089840548e2577bfb
  2. 算額(その770)
    https://blog.goo.ne.jp/r-de-r/e/28303bccaf422cc2177040e433f30ab2
  3. 算額(その769)
    https://blog.goo.ne.jp/r-de-r/e/1fde4dc2613ad46677b38c699f071794
  4. 算額(その697)
    https://blog.goo.ne.jp/r-de-r/e/ab7929468b8b6664735d652fee031327
  5. 算額(その696)
    https://blog.goo.ne.jp/r-de-r/e/6b46a88a2e9e196cf9e0c035d3794167
  6. 算額(その695)
    https://blog.goo.ne.jp/r-de-r/e/5ce9b6cd24b57fd7d01497b6aa9a1331
  7. 算額(その694)
    https://blog.goo.ne.jp/r-de-r/e/d8227dd71b21273a0effe1e59029b573
  8. 算額(その693)
    https://blog.goo.ne.jp/r-de-r/e/61cdf4ecbf7104c1f85b9f86e3079710
  9. 算額(その692)
    https://blog.goo.ne.jp/r-de-r/e/1b3e7a7043254fc1b5c573ba784967f8
  10. 算額(その691)
    https://blog.goo.ne.jp/r-de-r/e/17c7f558bedea88e621af02d3541043f
  11. 算額(その690)
    https://blog.goo.ne.jp/r-de-r/e/30af0ab7a44f52601537b571991c8fac
  12. 算額(その689)
    https://blog.goo.ne.jp/r-de-r/e/4ae4441cff6e545596ad2a364157d53d
  13. 算額(その687)
    https://blog.goo.ne.jp/r-de-r/e/ebbd04fd41bc039e531d8e67f0eea1ae
  14. 算額(その686)
    https://blog.goo.ne.jp/r-de-r/e/20fbf15bdacfff2a94c780f3cf0fa0f3
  15. 算額(その650)
    https://blog.goo.ne.jp/r-de-r/e/141375a967a505289641a98d59bd5f58
  16. 算額(その649)
    https://blog.goo.ne.jp/r-de-r/e/a89c101131a973c5fcd4bd2c2b31feaf
  17. 算額(その575)
    https://blog.goo.ne.jp/r-de-r/e/36a7228a5c3c1b7ae0189f08334ed5e3
  18. 算額(その447)
    https://blog.goo.ne.jp/r-de-r/e/e3b4d9e0d5fdc6c578c6925d586b183d

続神壁算法

  1. 算額(その792)
    https://blog.goo.ne.jp/r-de-r/e/5fc1e2f9b4eeee6953c3cce3295a5fc0
  2. 算額(その791)
    https://blog.goo.ne.jp/r-de-r/e/d74e1016f980b8f6556f77af065a2b52
  3. 算額(その790)
    https://blog.goo.ne.jp/r-de-r/e/dbb15ebd39534aeb3b00793f7d864dc3
  4. 算額(その789)
    https://blog.goo.ne.jp/r-de-r/e/c1965824c3ffb56527ca0f472779b132
  5. 算額(その788)
    https://blog.goo.ne.jp/r-de-r/e/19175b3fed4b6ba45257a043104d7a49
  6. 算額(その787)
    https://blog.goo.ne.jp/r-de-r/e/01ca15ba267aea8104461af6a02a4afe
  7. 算額(その446)
    https://blog.goo.ne.jp/r-de-r/e/0743744b8d4a5067267df6f06f239e8b
  8. 算額(その445)
    https://blog.goo.ne.jp/r-de-r/e/9f35b3c1b37302b50c61355f43956949
  9. 算額(その444)
    https://blog.goo.ne.jp/r-de-r/e/62b830b5a07cbc2f693a80ab9d29d1e0
  10. 算額(その443)
    https://blog.goo.ne.jp/r-de-r/e/08f51fefb42a8dc2a610944e746bc186
  11. 算額(その442)
    https://blog.goo.ne.jp/r-de-r/e/6cad986da6c2eaad06db77453c767746
  12. 算額(その441)
    https://blog.goo.ne.jp/r-de-r/e/7effc64587864aa350523addf0442972
  13. 算額(その440)
    https://blog.goo.ne.jp/r-de-r/e/59a28abd7846fcf2d0cb8e257673a330
  14. 算額(その439)
    https://blog.goo.ne.jp/r-de-r/e/8a2c65a3828a2ca879d7fb5b86444a12
  15. 算額(その438)
    https://blog.goo.ne.jp/r-de-r/e/5411ee85b446d73ca00e82a2e26be622
  16. 算額(その437)
    https://blog.goo.ne.jp/r-de-r/e/7d06276e0130f0c09bb89c49d1fc668d
  17. 算額(その436)
    https://blog.goo.ne.jp/r-de-r/e/25ceda8b47f762c54a441ebd059cd781
  18. 算額(その394)
    https://blog.goo.ne.jp/r-de-r/e/10fe593496313a52b495d0080109b0c1

神壁算法追加

  1. 算額(その515)
    https://blog.goo.ne.jp/r-de-r/e/894cc97a5779e7592d5f0a58cabf578e
  2. 算額(その514)
    https://blog.goo.ne.jp/r-de-r/e/548cd85b7970aa08c54aa0cfc63cc96d

 

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

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

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