裏 RjpWiki

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

算額(その113)

2023年01月27日 | Julia

算額(その113)

(16) 京都府京都市右京区山ノ内宮脇町 山王神社 明治23年(1890)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

愛媛県松山市吉藤町 三島神社 明治13年(1880)7月
http://www.wasan.jp/ehime/misima.html

兵庫県上郡町 上郡天満神社 文化5年(1808)
https://kamigori-town-museum.jp/cultural-asetts/559/

鉤が 108間6合,股が144間8合の直角三角形の内部に円と正三角形が入っている。正三角形の一辺の長さはいかほどか?

鉤股を 1086, 1448 とする。
円の半径を r,正三角形の頂点座標を,x1, (x2, y2), x3 として以下のように記号を定め方程式を解く。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す!!
end;

@syms r, x1, x2, y2;
eq1 = distance(0, 1086, 1448, 0, r, r) - r^2;
eq2 = (x2 - x1)^2 + y2^2  - 4(x2 - x1)^2 |> expand
eq3 = distance(x1, 0, x2, y2, r, r) - r^2
eq4 = 1448y2 - 1086*(1448 - x2)

res = solve([eq1, eq2, eq3, eq4], (r, x1, x2, y2));
res

   8-element Vector{NTuple{4, Sym}}:
    (362, 362 - 362*sqrt(3), 5792/13 - 1448*sqrt(3)/13, 1086*sqrt(3)/13 + 9774/13)
    (362, 362 - 362*sqrt(3)/3, -1448*sqrt(3)/3, 362*sqrt(3) + 1086)
    (362, 362 + 362*sqrt(3), 1448*sqrt(3)/13 + 5792/13, 9774/13 - 1086*sqrt(3)/13)
    (362, 362*sqrt(3)/3 + 362, 1448*sqrt(3)/3, 1086 - 362*sqrt(3))
    (2172, 2172 - 2172*sqrt(3), 4344 - 2896*sqrt(3), -2172 + 2172*sqrt(3))
    (2172, 2172 - 724*sqrt(3), 21720/13 - 8688*sqrt(3)/13, -2172/13 + 6516*sqrt(3)/13)
    (2172, 2172 + 2172*sqrt(3), 4344 + 2896*sqrt(3), -2172*sqrt(3) - 2172)
    (2172, 724*sqrt(3) + 2172, 8688*sqrt(3)/13 + 21720/13, -6516*sqrt(3)/13 - 2172/13)

for i = 1:8
   println("r  = $(res[i][1].evalf())")
   println("x1 = $(res[i][2].evalf())")
   println("x2 = $(res[i][3].evalf())")
   println("y2 = $(res[i][4].evalf())")
   println()
end

   r  = 362.000000000000
   x1 = -265.002392339934
   x2 = 252.614648510790
   y2 = 896.539013616908

   r  = 362.000000000000
   x1 = 152.999202553355
   x2 = -836.003189786578
   y2 = 1713.00239233993

   r  = 362.000000000000
   x1 = 989.002392339934
   x2 = 638.462274566133
   y2 = 607.153294075400

   r  = 362.000000000000
   x1 = 571.000797446645
   x2 = 836.003189786578
   y2 = 458.997607660066

   r  = 2172.00000000000
   x1 = -1590.01435403960
   x2 = -672.019138719469
   y2 = 1590.01435403960

   r  = 2172.00000000000
   x1 = 917.995215320133
   x2 = 513.226352603200
   y2 = 701.080235547600

   r  = 2172.00000000000
   x1 = 5934.01435403960
   x2 = 9360.01913871947
   y2 = -5934.01435403960

   r  = 2172.00000000000
   x1 = 3426.00478467987
   x2 = 2828.31210893526
   y2 = -1035.23408170145

4 番目の解が妥当である。

正三角形の一辺の長さは x1 + 2(x2-x1) = 530.004784679867,もとの単位では 53間強である。

144間8合 × (√3 - 1)/2 = 53.00047846798671 間

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
  plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x1, x2, y2) = res[4]
   x3 = x1 + 2(x2 - x1)
   plot([0, 1448, 0, 0], [0, 0, 1086, 0], color=:black, lw=0.5)
   circle(r, r, r)
   segment(x1, 0, x2, y2, :green)
   segment(x2, y2, x3, 0, :green)
   if more == true
       println("x = $(sqrt((x2-x1)^2 + y2^2)) = $((sqrt((x2-x1)^2 + y2^2)).evalf())")
       println("x3 - x1 = $(x3 - x1) = $((x3 - x1).evalf())")
       point(r, r, "(r,r)", :green, :center)
       point(x1, 0, " x1", :green, :left, :bottom)
       point(x3, 0, " x3", :green, :left, :bottom)
       point(x2, y2, " (x2,y2)", :green, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   x = sqrt((-362 + 362*sqrt(3))^2 + (1086 - 362*sqrt(3))^2) = 530.004784679867
   x3 - x1 = -724 + 724*sqrt(3) = 530.004784679867

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

算額(その112)

2023年01月27日 | Julia

算額(その112)

愛媛県松山市 伊佐爾波神社 嘉永3年(1850)5月
http://www.wasan.jp/ehime/isaniwa14.html

円内に 2 つの甲円,3 つの乙円が入っている。それぞれの円の半径を求めよ。

算額(その111)に似ているのだが,写真が不鮮明であるが,下の 2 つの乙円と甲円がそれぞれ直線に接している。

外円の半径を 1 とする。甲円,乙円の中心座標と半径を (r1, y1, r1),(r2, y2, r2) とおき,以下の方程式を立てて解く。
算額(その111)と違うのは eq4 である。しょうもない条件であるが,これでは以下のように無意味な解しか得られない。

using SymPy
@syms r1, r2, y1, y2;
eq3 = r1^2 + (1 - r2 - y1)^2 - (r1 + r2)^2 |> expand
eq4 = y1 - (y2 + r2) - r1
eq5 = r1^2 + y1^2 - (1 - r1)^2 |> expand
eq6 = r2^2 + y2^2 - (1 - r2)^2 |> expand;

solve([eq3, eq4, eq5, eq6], (r1, r2, y1, y2))

   1-element Vector{NTuple{4, Sym}}:
    (0, 0, 1, 1)

そこで, nlsolve() で数値解を求めることにする。

eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println

   -2*r1*r2 + 2*r2*y1 - 2*r2 + y1^2 - 2*y1 + 1
   -r1 - r2 + y1 - y2
   2*r1 + y1^2 - 1
   2*r2 + y2^2 - 1

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

function H(u)
  (r1, r2, y1, y2) = u
  return [
       -2*r1*r2 + 2*r2*y1 - 2*r2 + y1^2 - 2*y1 + 1,
       -r1 - r2 + y1 - y2,
       2*r1 + y1^2 - 1,
       2*r2 + y2^2 - 1
   ]
end;

iniv = [0.5, 0.3, 0.1, -0.6]
res = nls(H, ini=iniv)
println(res)

ちゃんと解が求まった。

   ([0.49309632821478017, 0.28307747532588534, 0.11750465339908704, -0.6586691501415785], true)

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, y1, y2) = [0.49309632821478017, 0.28307747532588534, 0.11750465339908704, -0.6586691501415785]
   plot()
   circle(0, 0, 1)
   circle(r2, y2, r2, :blue)
   circle(-r2, y2, r2, :blue)
   circle(0, 1 - r2, r2, :blue)
   circle(r1, y1, r1, :green)
   circle(-r1, y1, r1, :green)
   hline!([y2+r2], color=:black, lw=0.5)
   if more == true
       @printf("r1 = %.5f; r2 = %.5f;  y1 = %.5f;  y2 = %.5f", r1, r2, y1, y2)
       point(r1, y1, "(r1,y1)", :green, :center)
       point(0, 1 - r2, "1-r2 ", :blue, :right)
       point(r2, y2, "(r2,y2)", :blue, :center)
       point(0, y2+r2, "y2+r2", :blue, :center)
       point(r1, y2+r2, "(r1, y2+r2)")
       hline!([0, y2+r2], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.49310; r2 = 0.28308;  y1 = 0.11750;  y2 = -0.65867

 

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

算額(その111)

2023年01月27日 | Julia

算額(その111)

愛媛県松山市 伊佐爾波神社 嘉永3年(1850)正月
http://www.wasan.jp/ehime/isaniwa7.html

円内に 2 つの甲円,3 つの乙円が入っている。それぞれの円の半径を求めよ。

外円の半径を 1 とする。甲円,乙円の中心座標と半径を (r1, y1, r1),(r2, y2, r2) とおき,以下の方程式を立てて解く。

using SymPy
@syms r1, r2, y1, y2;
eq3 = r1^2 + (1 - r2 - y1)^2 - (r1 + r2)^2
eq4 = (r1 - r2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq5 = r1^2 + y1^2 - (1 - r1)^2
eq6 = r2^2 + y2^2 - (1 - r2)^2;

res = solve([eq3, eq4, eq5, eq6], (r1, r2, y1, y2))

   5-element Vector{NTuple{4, Sym}}:
    (0, 0, 1, 1)
    ((-142*sqrt(2) - 197 + 65*sqrt(-13 + 16*sqrt(2)) + 46*sqrt(-26 + 32*sqrt(2)))/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -16*sqrt(2)/49 + 2*sqrt(-13 + 16*sqrt(2))/49 + 13/49 + 4*sqrt(-26 + 32*sqrt(2))/49, (-10*sqrt(-26 + 32*sqrt(2)) - 13*sqrt(-13 + 16*sqrt(2)) + 30*sqrt(2) + 43)/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -4*sqrt(2)/7 - 2/7 + sqrt(-13 + 16*sqrt(2))/7)
    (-(197 + 142*sqrt(2) + 65*sqrt(-13 + 16*sqrt(2)) + 46*sqrt(-26 + 32*sqrt(2)))/(sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -16*sqrt(2)/49 - 4*sqrt(-26 + 32*sqrt(2))/49 - 2*sqrt(-13 + 16*sqrt(2))/49 + 13/49, (13*sqrt(-13 + 16*sqrt(2)) + 30*sqrt(2) + 43 + 10*sqrt(-26 + 32*sqrt(2)))/(sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9), -4*sqrt(2)/7 - sqrt(-13 + 16*sqrt(2))/7 - 2/7)
    ((-65*sqrt(13 + 16*sqrt(2)) + 46*sqrt(26 + 32*sqrt(2)) - 142*sqrt(2)*I + 197*I)/(sqrt(13 + 16*sqrt(2)) - 9*I + 4*sqrt(2)*I), 13/49 + 16*sqrt(2)/49 - 2*sqrt(-16*sqrt(2) - 13)/49 + 4*sqrt(-32*sqrt(2) - 26)/49, (-10*sqrt(26 + 32*sqrt(2)) + 13*sqrt(13 + 16*sqrt(2)) - 43*I + 30*sqrt(2)*I)/(sqrt(13 + 16*sqrt(2)) - 9*I + 4*sqrt(2)*I), -2/7 + 4*sqrt(2)/7 - I*sqrt(13 + 16*sqrt(2))/7)
    ((-65*sqrt(13 + 16*sqrt(2)) + 46*sqrt(26 + 32*sqrt(2)) - 197*I + 142*sqrt(2)*I)/(sqrt(13 + 16*sqrt(2)) - 4*sqrt(2)*I + 9*I), 13/49 + 16*sqrt(2)/49 - 4*sqrt(-32*sqrt(2) - 26)/49 + 2*sqrt(-16*sqrt(2) - 13)/49, (-10*sqrt(26 + 32*sqrt(2)) + 13*sqrt(13 + 16*sqrt(2)) - 30*sqrt(2)*I + 43*I)/(sqrt(13 + 16*sqrt(2)) - 4*sqrt(2)*I + 9*I), -2/7 + 4*sqrt(2)/7 + I*sqrt(13 + 16*sqrt(2))/7)

5 組の解が得られるが,2番目のものだけが適切である。なお,@syms において各変数を ::positive と宣言すると解が求まらないので注意。

for i = 1:5
   println("r1 = $(res[i][1].evalf())")
   println("r2 = $(res[i][2].evalf())")
   println("y1 = $(res[i][3].evalf())")
   println("y2 = $(res[i][4].evalf())")
   println()
end

   r1 = 0
   r2 = 0
   y1 = 1.00000000000000
   y2 = 1.00000000000000

   r1 = 0.494520180083009
   r2 = 0.288374102478171
   y1 = 0.104688298457764
   y2 = -0.650578046850383

   r1 = -45.1219371780525
   r2 = -0.681329898313661
   y1 = 9.55216595103462
   y2 = -1.53709459586173

   r1 = 0.31370849898476 - 0.463999436044365*I
   r2 = 0.727090142815704 + 0.445454898991966*I
   y1 = -0.82842712474619 - 0.560096865715887*I
   y2 = 0.522407749927483 - 0.852695809075959*I

   r1 = 0.31370849898476 + 0.463999436044365*I
   r2 = 0.727090142815704 - 0.445454898991966*I
   y1 = -0.82842712474619 + 0.560096865715887*I
   y2 = 0.522407749927483 + 0.852695809075959*I

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")

   (r1, r2, y1, y2) = (
       (-142*sqrt(2) - 197 + 65*sqrt(-13 + 16*sqrt(2)) + 46*sqrt(-26 + 32*sqrt(2)))/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9),
       -16*sqrt(2)/49 + 2*sqrt(-13 + 16*sqrt(2))/49 + 13/49 + 4*sqrt(-26 + 32*sqrt(2))/49,
       (-10*sqrt(-26 + 32*sqrt(2)) - 13*sqrt(-13 + 16*sqrt(2)) + 30*sqrt(2) + 43)/(-sqrt(-13 + 16*sqrt(2)) + 4*sqrt(2) + 9),
       -4*sqrt(2)/7 - 2/7 + sqrt(-13 + 16*sqrt(2))/7)

   plot()
   circle(0, 0, 1)
   circle(r2, y2, r2, :blue)
   circle(-r2, y2, r2, :blue)
   circle(0, 1 - r2, r2, :blue)
   circle(r1, y1, r1, :green)
   circle(-r1, y1, r1, :green)
   if more == true
       @printf("r1 = %.5f; r2 = %.5f;  y1 = %.5f;  y2 = %.5f", r1, r2, y1, y2)
       point(r1, y1, "(r1,y1)", :green, :center)
       point(0, 1 - r2, "1-r2 ", :blue, :right)
       point(r2, y2, "(r2,y2)", :blue, :center)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.49452; r2 = 0.28837;  y1 = 0.10469;  y2 = -0.65058

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

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

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