裏 RjpWiki

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

算額(その306)

2023年06月30日 | Julia

算額(その306)

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 永保寺 天保15年(1844)

問題文2 

三辺のの和が 26.14 丈となる鉤股弦(直角三角形)がある。鉤股の差は分からないが,弦は鉤より 4.44 丈長い。この直角三角形に内接するように正三角形を入れるとき,正三角形の一辺,中句,三角面を求めよ。

中句と弦の交点を (x, y),正三角形の他の 2 点を (x, y2), (x3, 0) とおき,以下の連立方程式を解く。

まず,鉤, 股, 弦, 中句 を求める。

include("julia-source.txt");

using SymPy
@syms 鉤::positive, 股::positive, 弦::positive,
   中句::positive, 三角面::positive;

eq1 = 弦 -  鉤 - 4.44
eq2 = 鉤 + 股 + 弦 - 26.14
eq3 = 鉤^2 + 股^2 - 弦^2
eq4 = 中句 * 弦 - 鉤 * 股  # 鉤股弦の面積算出法が2通りある

solve([eq1, eq2, eq3, eq4], (鉤, 股, 弦, 中句))

   1-element Vector{NTuple{4, Sym}}:
    (6.46022727742318, 8.77954544515363, 10.9002272774232, 5.20336480374436)

正確な数値はともかく,出題文を勘案すれば小数点以下3桁目で四捨五入すべきであろう。

(鉤, 股, 弦) = (6.46, 8.78, 10.90)
鉤 + 股 + 弦 |> println
弦 - 鉤      |> println
sqrt(6.46^2 + 8.78^2) |> println

   26.14
   4.44
   10.900458705944443

(鉤, 股, 弦, 中句) = (6.46, 8.78, 10.90, 5.20) とするのが妥当。

ちなみに,算額の「答」は記述間違い(6.46 を 6.64 とするなど),計算自体の間違い(弦を 11.1 とするなど)があり信頼性がない。

なお,以上の解法では図を描くことはできない(中句の長さだけでは図を描けない)ので,以下のような別解を探る。
まず,鉤,股,弦,中句の他に中句と弦の交点座標 (x, y) を求める。

using SymPy
@syms 鉤::positive, 股::positive, 弦::positive,
   中句::positive, x::positive, y::positive;

eq1 = 弦 -  鉤 - 4.44
eq2 = 鉤 + 股 + 弦 - 26.14
eq3 = 鉤^2 + 股^2 - 弦^2
eq4 = 中句 * 弦 - 鉤 * 股  # 鉤股弦の面積算出法が2通りある
eq5 = y - (股 - x)*鉤/股
eq6 = x^2 + y^2 - 中句^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (鉤, 股, 弦, x, y, 中句))

   1-element Vector{NTuple{6, Sym}}:
    (6.46022727742318, 8.77954544515363, 10.9002272774232, 3.08387324263936, 4.19102983813985, 5.20336480374436)

弦の上にある (x, y) を頂点の一つとして,鉤,股の上に他の2つの頂点 (0, y2),(x3, 0) がある正三角形の x3, y2 を求めるが,上述の eq1〜eq6 に2条件を加えて一挙に8元連立方程式を解くのは solve() では無理なので,

(鉤, 股, 弦, x, y, 中句) = (6.46022727742318, 8.77954544515363, 10.9002272774232, 3.08387324263936, 4.19102983813985, 5.20336480374436) が既知であるとして,以下の二元連立方程式を解く。

@syms x3::positive, y2::positive
(鉤, 股, 弦, x, y, 中句) = (6.46022727742318, 8.77954544515363, 10.9002272774232, 3.08387324263936, 4.19102983813985, 5.20336480374436)
eq7 = (x^2 + (y - y2)^2) - (x3^2 + y2^2)
eq8 = ((x - x3)^2 + y^2) - (x3^2 + y2^2)
res = solve([eq7, eq8], (x3, y2))

   1-element Vector{Tuple{Sym, Sym}}:
    (4.17520337305602, 1.15039530221371)

算額の三角面の答えは 3.939986... とあり,「三重県に現存する算額の研究」福島完では 3.809101139 と不一致を見せている。しかし,後者も現代的解法の 60 ページに,BD = DG + GB = √3/2 DF + 1/2 DF としているところがあるが,GB は DF/2 ではない。掲載している図も実際の位置関係を表すものではない。GB = GF としているがそれは成り立たない。他の算額も同じであるが,図を描くための全てのパラメータを求めて,実際に図を描いてみれば,得られたパラメータが妥当なものかどうか明らかになる。想定していた図と全く異なる図になることもすぐに分かるはず。三角面が 3.809101139 の正三角形は描けない

   鉤 = 6.460227
   股 = 8.779545
   弦 = 10.900227
   中句 = 5.203365
   x = 3.083873
   y = 4.191030
   三角面 = 4.330789
   y2 = 1.150395
   x3 = 4.175203

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, 弦, x, y, 中句) = (6.46022727742318, 8.77954544515363,
       10.9002272774232, 3.08387324263936, 4.19102983813985,
       5.20336480374436)
   (x3, y2) = (4.17520337305602, 1.15039530221371)
   三角面 = sqrt(x3^2 + y2^2)
   txt = @sprintf("鉤 = %.6f\n", 鉤) * @sprintf("股 = %.6f\n", 股) *
         @sprintf("弦 = %.6f\n", 弦) * @sprintf("中句 = %.6f\n", 中句) *
         @sprintf("x = %.6f\n", x) * @sprintf("y = %.6f\n", y) *
         @sprintf("三角面 = %.6f\n", 三角面) * @sprintf("y2 = %.6f\n", y2) *
         @sprintf("x3 = %.6f\n", x3)
   println(txt)
   plot([0, 0, 股, 0], [0, 鉤, 0, 0], color=:black, lw=0.5)
   segment(0, 0, x, y, :red)
   plot!([x, 0, x3, x], [y, y2, 0, y], color=:green, lw=0.5)
   if more
       point(0, 鉤, " 鉤", :black, :left, :bottom)
       point(股, 0, " 股", :black, :left, :bottom)
       point(x, y, " (x,y)", :black, :left, :bottom)
       point(0, y2, "  y2", :green, :left, :bottom)
       point(x3, 0, " x3", :green, :left, :bottom)
       annotate!(5.5, 4.0, text(txt, :left, :green, 10))
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その305)

2023年06月29日 | Julia

算額(その305)

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 菅原神社 嘉永7年(1854)

問題文4 

外大円の中に長方形を起き,菱形を切り取った残りの領域に円を置く。円の直径が 3 尺,長方形の残りの長さ BE が 4.5 尺である。このとき,菱形の一辺(菱面=BD),長軸(同立=AD),短軸(同横= BC),外大円の円周,矢(r0 - y)を求めよ。

外大円の半径を r0 とおく。r0 = 同立/2
長方形の右上の頂点の座標を (x, y) とする。
以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x::positive, y::positive,
   菱面::positive, 同立::positive, 同横::positive,
   外大円円周::positive, 矢::positive;

r1 = 3/2
eq1 = 4.5 + 2y - sqrt(4.5^2 + (2y)^2) - 2r1
eq2 = 4.5^2 + (2y)^2 - 菱面^2  # ⊿BDE において
eq3 = (菱面 + 4.5)^2 + (2y)^2 - 同立^2  # ⊿ADE において
eq4 = (同立/2)^2 + (同横/2)^2 - 菱面^2  # ⊿OBE において
eq5 = 外大円円周 - 同立*PI
eq6 = 矢 - (同立/2 - y);

eq1 〜 eq6 までの6元連立方程式にすると解けなくなってしまうので,eq1 〜 eq5 の5元連立方程式を解く。

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (y, 菱面, 同立, 同横, 外大円円周, 矢))

res = solve([eq1, eq2, eq3, eq4, eq5], (y, 菱面, 同立, 同横, 外大円円周))

   1-element Vector{NTuple{5, Sym}}:
    (3.00000000000000, 7.50000000000000, 13.4164078649987, 6.70820393249937, 42.1488883862444)

res[1][3]/2 - res[1][1] |> println  # 矢

   3.70820393249937

なお,solve(eq1) で y だけを求めれば,ドミノ倒しで残りの変数は確定する。

res = solve(eq1)[1] |> println  # y

   3.00000000000000

sqrt(6^2 + 4.5^2)  # 菱面 = sqrt((2y)^2 + 4.5^2)

   7.5

sqrt((7.5 + 4.5)^2 + (2*3)^2)  # 同立 = sqrt((菱面 + 4.5) + (2y)^2)

   13.416407864998739

sqrt(7.5^2 - (13.416407864998739/2)^2)*2  # 同横 = sqrt(菱面^2-(同立/2)^2)*2

   6.708203932499367

13.4164078649987*pi  # 外大円円周 = 同立*π

   42.14888838624424

(13.416407864998739/2) - 3  # 矢 = 同立/2 - y

   3.7082039324993694

   y = 3.000000
   菱面 = 7.500000
   同立 = 13.416408
   同横 = = 6.708204
   外大円円周 = = 42.148888
   矢 = 3.708204

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (y, 菱面, 同立, 同横, 外大円円周) = (3.00000000000000, 7.50000000000000, 13.4164078649987, 6.70820393249937, 42.1488883862444)
   矢 = 3.70820393249937
   @printf("y = %.6f\n", y)
   @printf("菱面 = %.6f\n", 菱面)
   @printf("同立 = %.6f\n", 同立)
   @printf("同横 = = %.6f\n", 同横)
   @printf("外大円円周 = = %.6f\n", 外大円円周)
   @printf("矢 = %.6f\n", 矢)
   plot([6, 6, -6, -6, 6], [-y, y, y, -y, -y], color=:black, lw=0.5)
   plot!([6, 1.5, -6, -1.5, 6], [-y, y, y, -y, -y], color=:red, lw=0.5)
   circle(4.5, 1.5, 1.5, :green)
   segment(1.5, y, -1.5, -y, :blue)
   segment(-6, y, 6, -y, :magenta)
   if more
       point(4.5, 1.5, "(4.5,1.5)", :green, :center)
       point(1.5, y, "B:(1.5,y) ", :red, :left, :bottom)
       point(6, y, "E:(6,y)", :black, :left, :bottom)
       point(-6, y, " A", :green, :left, :bottom)
       point(-1.5, -y, " C", :green, :left, :bottom)
       point(6, -y, " D", :green, :left, :bottom)
       point(0, 0, " O", :green, :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でシェアする

算額(その304)

2023年06月29日 | Julia

算額(その304)

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 菅原神社 嘉永7年(1854)

問題文1 

正方形に内接する円を描き,円の中に図のように菱形と2種類の円を2つずつ入れる。菱形の一辺と短軸は等しい。小円の直径が 3.55 尺のとき,正方形の一変,正方形に内接する円の円周,菱形の一変,菱形内の大円の直径,矢(菱形の上の頂点から正方形の上辺までの距離)を求めよ。

正方形の中心を原点とする。正方形の一辺の長さを 2r0 とする(内接する円の半径が r0)。菱形内の小円の半径,中心座標を r1, (r1, y1),大円の半径,中心座標を r2, (r2, 0) とする。
以下の連立方程式をとき,y1, r2, r0 を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, y1::positive, r2::positive, r0::positive, y::positive;

r1 = 3.55/2
sqrt3 = sqrt(Sym(3))
y = r0/sqrt3
y = r2*sqrt3
eq1 = (r2 - r1)^2 + y1^2 - (r2 + r1)^2
eq2 = y - y1 - sqrt3*r1
eq3 = 2y - sqrt(y^2 + r0^2);

res = solve([eq1, eq2, eq3], (y1, r2, r0))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (2*sqrt(3)*r1, 3*r1, 9*r1)

y1 = 2*sqrt(3)*r1
r2 = 3r1
r0 = 9r1
である。他の変数の値は,これらから求められる。

   中円の径 = 2r2 = 10.650000
   菱面 = 2y = 18.446341
   方面 = 2r0 = 31.950000
   円廻り = π*2r0 = 100.373885
   矢 = r0 - y = 6.751829

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 3.55/2
   y1 = 2*sqrt(3)*r1
   r2 = 3r1
   r0 = 9r1
   y = r0/sqrt(3)
   @printf("中円の径 = 2r2 = %.6f\n", 2r2)
   @printf("菱面 = 2y = %.6f\n", 2y)
   @printf("方面 = 2r0 = %.6f\n", 2r0)
   @printf("円廻り = π*2r0 = %.6f\n", pi*2r0)
   @printf("矢 = r0 - y = %.6f\n", r0 - y)
   plot([r0, r0, -r0, -r0, r0], [-r0, r0, r0, -r0, -r0], color=:black, lw=0.5)
   circle(0, 0, r0, :green)
   circle(r1, y1, r1, :blue)
   circle(-r1, y1, r1, :blue)
   circle(r2, 0, r2, :red)
   circle(-r2, 0, r2, :red)
   plot!([r0, 0, -r0, 0, r0], [0, y, 0, -y, 0], color=:magenta, lw=0.5)
   x = sqrt(r0^2 - y^2)
   segment(-x, y, x, y)
   if more
       point(0, y, " y", :magenta, :left, :bottom)
       point(0, r0, " r0", :green, :left, :bottom)
       point(r0, 0, " r0", :green, :left, :bottom)
       point(r2, 0, " 中円\n r2", :red, :left, :bottom)
       point(r1, y1, " 小円(r1,y1)", :blue, :left, :bottom)
       point(0, (r0 + y)/2, " 矢=r0-y", :black, mark=false)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その303)

2023年06月29日 | Julia

算額(その303)

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 菅原神社 嘉永7年(1854)

問題文3 

外側の正三角形の一辺の長さを 2a,一番大きい正方形,二番目,三番目の正方形の一辺の長さをそれぞれ 2b, 2c, 2d とし,以下の連立方程式を解く。2d = 7.8179 が与えられているが,d も未知数として方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive, d ::positive

d = 7.8179/2
sqrt3 = sqrt(Sym(3))
eq1 = sqrt3*a - 2(b + c + d) - sqrt3*d
eq2 = sqrt3*a - 2(b + c) - sqrt3*c
eq3 = sqrt3*a - 2b - sqrt3*b

res = solve([eq1, eq2, eq3], (a, b, c))

   Dict{Any, Any} with 3 entries:
     b => 4*sqrt(3)*d/3 + 7*d/3
     a => 5*d + 26*sqrt(3)*d/9
     c => d + 2*sqrt(3)*d/3

res[a] |> simplify |> println
res[b] |> simplify |> println
res[c] |> simplify |> println

   d*(45 + 26*sqrt(3))/9
   d*(4*sqrt(3) + 7)/3
   d*(3 + 2*sqrt(3))/3

二方(2 番目に大きい正方形の一辺の長さ)= 2c

res[c](d => 7.8179/2).evalf()*2

   16.8452333389952

一方(1 番目大きい正方形の一辺の長さ)= 2b

res[b](d => 7.8179/2).evalf()*2

   36.296433344657

外三角(外側の正三角形の一辺の長さ)

res[a](d => 7.8179/2).evalf()*2

   78.207944468979

三角面(内側の正三角形の一辺の長さ)
内側の正三角形は正方形に内接する円(半径 b)に内接する

(res[b](d => 7.8179/2) * sqrt3/2 * 2).evalf()

   31.4336333432415

小円径(内側の正三角形に内接する円の直径)b/2 * 2 である

(res[b](d => 7.8179/2) / 2 * 2).evalf()

   18.1482166723285

   a = 39.103972;  b = 18.148217;  c = 8.422617; d = 3.908950
   二方 = 2c = 16.845233
   一方 = 2b = 36.296433
   外三角 = 2a = 78.207944
   三角面 = √3b = 31.433633
   小円径 = b = 18.148217

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   d = 7.8179/2
   a = d*(45 + 26*sqrt(3))/9
   b = d*(4*sqrt(3) + 7)/3
   c = d*(3 + 2*sqrt(3))/3
   @printf("a = %.6f;  b = %.6f;  c = %6f; d = %.6f\n", a, b, c, d)
   @printf("二方 = 2c = %.6f\n", 2c)
   @printf("一方 = 2b = %.6f\n", 2b)
   @printf("外三角 = 2a = %.6f\n", 2a)
   @printf("三角面 = √3b = %.6f\n", √3b)
   @printf("小円径 = b = %.6f\n", b)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:black, lw=0.5)
   rect(-d, 2(b + c), d, 2(b + c + d), :red)
   rect(-c, 2b, c, 2(b + c), :blue)
   rect(-b, 0, b, 2b, :green)
   circle(0, b, b, :brown)
   circle(0, b, b/2, :orange)
   plot!([b√3/2, 0, -b√3/2, b√3/2], [b/2, 2b, b/2, b/2], color=:magenta, lw=0.5)
   if more
       point(b, 2b, " (b,2b)")
       point(c, 2(b + c), " (c,2b+2c)", :blue)
       point(d, 2(b + c + d), " (c,2b+2c+2d)", :red)
       point(a, 0, " a", :black, :left, :bottom)
       point(0, √3a, " √3a", :brown, :left, :bottom)
       point(0, b, " b", :brown)
       point(0, b/2, " b/2", :orange)
       point(√3b/2, b/2, "  (√3b/2,b/2)", :magenta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その302)

2023年06月29日 | Julia

算額(その302)

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 菅原神社 嘉永7年(1854)

問題文5 面積が 163350 歩平方の鈎股弦がある。鈎股弦の長さの差は分からないが,弦が 82.5 尺であるとき,鈎,股,中勾,方面,小面を求めよ。
注: 1 尺 = 10 歩

正方形の一辺の長さを a として,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive, a ::positive

eq1 = 鈎 * 股 / 2 - 1633.50
eq2 = 鈎^2 + 股^2 - 82.5^2
eq3 = (鈎 - a)/a - a/(股 - a)

solve([eq1, eq2, eq3], (鈎, 股, a))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (49.5000000000000, 66.0000000000000, 28.2857142857143)
    (66.0000000000000, 49.5000000000000, 28.2857142857143)

鈎 < 股 ゆえ,最初の組が適解。

鈎 = 49.500000
股 = 66.000000
中勾(正方形の対角線) = 40.002041
方面(正方形の一辺の長さ) = 28.285714
小面(小さな正方形の一辺の長さ) = 20.001020

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, a) = (49.5000000000000, 66.0000000000000, 28.2857142857143)
   @printf("鈎 = %.6f\n", 鈎)
   @printf("股 = %.6f\n", 股)
   @printf("中勾(正方形の対角線) = %.6f\n", √2a)
   @printf("方面(正方形の一辺の長さ) = %.6f\n", a)
   @printf("小面(小さな正方形の一辺の長さ) = %.6f\n", a/√2)
   plot([0, 股 , 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   rect(0, 0, a, a, :red)
   plot!([a/2, a, a/2, 0, a/2], [0, a/2, a, a/2, 0], color=:blue, lw=0.5)
   if more
       point(0, 鈎, " 鈎", :black, :left, :bottom)
       point(股, 0, " 股", :black, :left, :bottom)
       point(a, 0, " a", :red, :left, :bottom)
       point(a, a/2, " (a,a/2)", :blue, :left, :bottom)
       point(a/2, a, " (a/2,a)", :blue, :left, :bottom)
       point(a, a, " (a,a)", :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でシェアする

算額(その301)

2023年06月28日 | Julia

算額(その301)

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 永保寺 天保15年(1844)

正三角形の中に図のように円,正方形が入っている。正三角形の面積が 2625 平方寸のとき,それぞれの大きさを求めよ。

正三角形の一辺の長さを 2a とする。
下の同心円の中心座標は (0, a/2),大きい円の半径 r1 は a/2 である。
上にある円の半径を r2 とし,以下の方程式を解き,a, r1, r2 を求める。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive
eq1 = a*sqrt(Sym(3))a - 2625
eq2 = a*tan(PI/6) - r1
eq3 = (sqrt(Sym(3))a - r2 - 2r1)sin(PI/6) - r2;

今回は,正三角形の面積が 2625 平方寸ということで,a が特定されている。

res = solve(eq1, a)
res[1] |> println

   5*3^(1/4)*sqrt(35)  # a

一般解を求める場合には以下の2元連立方程式を解く。

res = solve([eq2, eq3], (r1, r2))

   Dict{Any, Any} with 2 entries:
     r2 => sqrt(3)*a/9
     r1 => sqrt(3)*a/3

ついで,図中の b 及び小円の半径 r3 を求める。

   a = 38.929994;  r1 = 22.476243;  r2 = 7.492081
   b = 15.893104;  r3 = 11.238121
   三角面(正三角形の一辺の長さ)= 77.85998861
   中勾(三角形の高さ)= 67.42872808
   下円径(大きい円の直径)= 44.95248538
   上円径(上の円の直径)= 14.98416179
   大方面(大きい正方形の一辺の長さ) = 31.78620725
   小方面(小さい正方形の一辺の長さ) = 22.47624269
   小円径(小円の直径) = 22.47624269
   小円廻り(小円の円周) = 70.61119892

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 5*3^(1/4)*sqrt(35)
   r1 = sqrt(3)*a/3
   r2 = sqrt(3)*a/9
   b = r1*cos(PI/4)
   r3 = b*cos(PI/4)
   @printf("a = %.6f;  r1 = %.6f;  r2 = %.6f\nb = %.6f;  r3 = %.6f\n", a, r1, r2, b, r3)
   @printf("三角面(正三角形の一辺の長さ)= %.8f\n", 2a)  # 77.85998861
   @printf("中句(三角形の高さ)= %.8f\n", a * √3)  # 67.42872808
   @printf("下円径(大きい円の直径)= %.8f\n", sqrt(3)*a/3 * 2)  # 44.95248538
   @printf("上円径(上の円の直径)= %.8f\n", sqrt(3)*a/9 * 2)  # 14.98416179
   @printf("大方面(大きい正方形の一辺の長さ) = %.8f\n", 2r1*cos(pi/4))  # 31.78620725
   @printf("小方面(小さい正方形の一辺の長さ) = %.8f\n", 2(b/√2))  # 22.47624269
   @printf("小円径(小円の直径) = %.8f\n", 2r3)  # 70.61119892087622
   @printf("小円廻り(小円の円周) = %.8f\n", 2r3*pi)  # 70.61119892087622
   plot([a, 0 , -a, a], [0, √3a, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1)
   circle(0, r2 + 2r1, r2, :blue)
   rect(-b, r1 - b, b, r1 + b, :magenta)
   plot!([b, 0, -b, 0, b], [r1, r1 + b, r1, r1 - b, r1], color=:green, lw=0.5)
   circle(0, r1, r3, :orange)
   if more
       point(0, √3a, " √3a", :black, :left, :bottom)
       point(a, 0, " a", :black, :left, :bottom)
       point(0, r1, " r1", :red, :left, :bottom)
       point(0, r1+r3, "\n r1+r3", :orange, :left, :top)
       point(0, r1 + b, " r1+b", :green, :left, :bottom)
       point(0, r1 - b, " r1-b", :green, :left, :top)
       point(b, r1, "\n(b,r1)", :green, :left, :top)
       point(0, 2r1 + r2, " 2r1+r2", :blue)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その300)

2023年06月26日 | Julia

算額(その300)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(80)
長野県長野市若槻東條 雁田八幡宮 文化3年(1806) 文面一部欠損
県内の算額(274)
長野県長野市北尾張部 連開庵 奉納年不明

二等辺三角形に全円が内接し,斜線で区切られた三角形にそれぞれ等円が内接している。全円の直径が 3 寸,中鉤(高さ)が 4 寸のとき,等円の直径を求めよ。

5元連立方程式は SymPy の solve() では解けないが,

include("julia-source.txt");

using SymPy

@syms h::positive, a::positive, w::positive, x0::positive,
     r0::positive, r1::positive, x1::positive, x2::positiv;

h = 4
r0 = 3//2
x = a/2
l = sqrt(x^2 + h^2)
m = sqrt((a - x)^2 + h^2)
eq1 = (2w + l)r0 - w*h
eq2 = (a + 2l)r1 + (2w - a + l)r1 - w*h
eq3 = (w - x)^2 + h^2 - w^2
eq4 = distance(x, h, w, 0, x2, r1) - r1^2
eq5 = distance(x, h, w, 0, x0, r0) - r0^2;
@syms d
eq4 = apart(eq4, d) |> numerator
eq5 = apart(eq5, d) |> numerator;

# solve([eq1, eq2, eq3, eq4, eq5], (a, w, x0, x2, r1))

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")

   -w + 3*sqrt(a^2/4 + 16)/2,  # eq1
   r1*(a + 2*sqrt(a^2/4 + 16)) + r1*(-a + 2*w + sqrt(a^2/4 + 16)) - 4*w,  # eq2
   -w^2 + (-a/2 + w)^2 + 16,  # eq3
   16*a*r1*w - 16*a*r1*x2 - 64*r1^2 - 32*r1*w^2 + 32*r1*w*x2 + 64*w^2 - 128*w*x2 + 64*x2^2,  # eq4
   24*a*w - 24*a*x0 + 16*w^2 - 80*w*x0 + 64*x0^2 - 144,  # eq5

[eq1, eq3] は他とは独立しているので,これをまず解く。

solve([eq1, eq3], (a, w))

   1-element Vector{Tuple{Sym, Sym}}:
    (2*sqrt(2), 9*sqrt(2)/2)

a, w を既知として方程式 [eq2, eq4, eq5] を組み直す。

@syms h::positive, x0::positive,
     r0::positive, r1::positive, x1::positive, x2::positiv;

h = 4
r0 = 3//2
a = 2sqrt(Sym(2))    # 既知となった
w = 9sqrt(Sym(2))/2  # 既知となった
x = a/2
l = sqrt(x^2 + h^2)
m = sqrt((a - x)^2 + h^2)
eq1 = (2w + l)r0 - w*h
eq2 = (a + 2l)r1 + (2w - a + l)r1 - w*h
eq3 = (w - x)^2 + h^2 - w^2
eq4 = distance(x, h, w, 0, x2, r1) - r1^2
eq5 = distance(x, h, w, 0, x0, r0) - r0^2;
@syms d
eq4 = apart(eq4, d) |> numerator
eq5 = apart(eq5, d) |> numerator;

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")

   0,  # eq1
   18*sqrt(2)*r1 - 18*sqrt(2),  # eq2
   0,  # eq3
   -32*r1^2/81 + 56*sqrt(2)*r1*x2/81 - 56*r1/9 + 32*x2^2/81 - 32*sqrt(2)*x2/9 + 16,  # eq4
   32*x0^2/81 - 68*sqrt(2)*x0/27 + 52/9,  # eq5

解は 4 組得られる。

res = solve([eq2, eq4, eq5], (x0, x2, r1))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    (3*sqrt(2)/2, 5*sqrt(2)/2, 1)
    (3*sqrt(2)/2, 19*sqrt(2)/4, 1)
    (39*sqrt(2)/8, 5*sqrt(2)/2, 1)
    (39*sqrt(2)/8, 19*sqrt(2)/4, 1)

names = ("x0", "x2", "r1")
for i in 1:length(res)
   println("result $i")
   for (j, name) in enumerate(names)
       @printf("  %s = %.10f\n", name, res[i][j])
   end
end

   result 1
     x0 = 2.1213203436
     x2 = 3.5355339059
     r1 = 1.0000000000
   result 2
     x0 = 2.1213203436
     x2 = 6.7175144213
     r1 = 1.0000000000
   result 3
     x0 = 6.8942911166
     x2 = 3.5355339059
     r1 = 1.0000000000
   result 4
     x0 = 6.8942911166
     x2 = 6.7175144213
     r1 = 1.0000000000

x0 , x2 < 6 < w ゆえ,1 組目のものが適解で,等円の半径は 1 ゆえ,直径は 2 寸
a = 2*sqrt(2)
w = 9*sqrt(2)/2
x0 = 3*sqrt(2)/2
x2 = 5*sqrt(2)/2
r1 = 1

数値解は nlsolve() で求めることもできる。

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)
   (a, w, x0, x2, r1) = u
   return [
       -w + 3*sqrt(a^2/4 + 16)/2,  # eq1
       r1*(a + 2*sqrt(a^2/4 + 16)) + r1*(-a + 2*w + sqrt(a^2/4 + 16)) - 4*w,  # eq2
       -w^2 + (-a/2 + w)^2 + 16,  # eq3
       16*a*r1*w - 16*a*r1*x2 - 64*r1^2 - 32*r1*w^2 + 32*r1*w*x2 + 64*w^2 - 128*w*x2 + 64*x2^2,  # eq4
       24*a*w - 24*a*x0 + 16*w^2 - 80*w*x0 + 64*x0^2 - 144,  # eq5
   ]
end;

r3 = 1
h = 4
r0 = 3//2
iniv = [big"2.8", 6.3, 2.15, 3.5, 0.98]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[2.828427124746190097603377448419396157139390481511502720472215083707481854103854, 6.363961030678927719607599258943641353563445321336646285922337820157100414824585, 2.121320343559642573202533086283528300895624077402693816025820728389215128629599, 3.535533905932737622004221810468008036286461682938140594670174743019746380147886, 1.000000000000000000000000000000000000000569878492122589984801745785913706379155], true)

names = ("a", "w", "x0", "x2", "r1")
for (i, name) in enumerate(names)
   @printf("%2s = %.6f\n", name, res[1][i])
end

    a = 2.828427
    w = 6.363961
   x0 = 2.121320
   x2 = 3.535534
   r1 = 1.000000

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (h, r0) = (4, 3//2)
   (a, w, x0, x2, r1) = (2*sqrt(2), 9*sqrt(2)/2, 3*sqrt(2)/2, 5*sqrt(2)/2, 1)
   x = a/2
   plot([0, w, x, 0], [0, 0, h, 0], color=:black, lw=0.5)
   circle(x0, r0, r0)
   circle(x, r1, r1, :blue)
   circle(x2, r1, r1, :blue)
   segment(x, h, a, 0)
   if more
       point(x0, r0, " 全円:r0,(x0,r0)\n", :red, :left, :bottom)
       point(x, r1, " 等円:r1\n (x,r1)", :blue)
       point(x2, r1, " 等円:r1\n (x2,r1)", :blue)
       point(x, h, " (x,h)", :black, :left, :bottom)
       point(a, 0, " a", :black, :left, :bottom)
       point(w, 0, " w", :black, :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でシェアする

算額(その299)

2023年06月25日 | Julia

算額(その299)

愛媛県 伊佐爾波神社 明治6年
高阪金次郎峻則
http://www.wasan.jp/ehime/isaniwa22.html
https://www.city.matsuyama.ehime.jp/kanko/kankoguide/rekishibunka/bunkazai/ken/isaniwa_sangaku.html
https://isaniwa.official.jp/2012/12/28/%E7%AE%97%E9%A1%8D%E3%81%AB%E6%8C%91%E6%88%A6/

扇面(中心角 120 度)に東円 1 個,西円,南円,北円が 2 個ずつ入っている。南円の直径を知って北円の直径を求めよ。

以下の方程式を解けばよいが,一挙に解くのが難しい。

方程式は [eq4, eq7] と残りの方程式は(r0 を除いて)互いに独立である。

まず [eq4, eq7] を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, r3::positive, r4::positive,
     x2::positive, x3::positive, x4::positive, y4::positive;

r1 = r0/4
eq1 = x2^2 + (r0 - r1 - r0/2 - r2)^2 - (r1 + r2)^2
eq2 = x4^2 + (y4 - r0 + r1)^2 - (r1 + r4)^2
eq3 = (x2 - x4)^2 + (y4 - r0/2 - r2)^2 - (r2 + r4)^2
eq4 = x3^2 + (r0/2 - r3)^2 - (r3 + r0//2)^2
eq5 = x2^2 + (r0/2 + r2)^2 - (r0 - r2)^2
eq6 = x4^2 + y4^2 - (r0 - r4)^2
eq7 = r3/(sqrt(Sym(3))*r0/2 - x3) - tan(PI/12);

まず,eq4, eq7 は 未知数 r3, x3 を r0 で表して解くことができる。

res = solve([eq4, eq7], (r3, x3))

   1-element Vector{Tuple{Sym, Sym}}:
    (3*r0*(7 - 4*sqrt(3))/2, r0*(-3 + 2*sqrt(3)))

次に,残りの方程式のセット [eq1, eq2, eq3, eq5, eq6] から (r2, x2, r4, x4, y4) を求める(r0 を記号として含む)。

東円の半径は扇の半径の 1/4 である。

@syms r0::positive, r1::positive, r2::positive, r3::positive, r4::positive,
     x2::positive, x3::positive, x4::positive, y4::positive;

r1 = r0//4
eq1 = x2^2 + (r0 - r1 - r0//2 - r2)^2 - (r1 + r2)^2
eq2 = x4^2 + (y4 - r0 + r1)^2 - (r1 + r4)^2
eq3 = (x2 - x4)^2 + (y4 - r0//2 - r2)^2 - (r2 + r4)^2
eq5 = x2^2 + (r0//2 + r2)^2 - (r0 - r2)^2
eq6 = x4^2 + y4^2 - (r0 - r4)^2;

res = solve([eq1, eq2, eq3, eq5, eq6], (r2, x2, r4, x4, y4))

   1-element Vector{NTuple{5, Sym}}:
    (3*r0/16, r0*(-3 + 14*sqrt(3))*sqrt(28*sqrt(3) + 199)/772, 3*r0*(25 - 12*sqrt(3))/193, r0*sqrt(336*sqrt(3)/37249 + 2388/37249), r0*(68/193 + 60*sqrt(3)/193))

res[1][1] |> simplify |> println  # r2 西円の半径

   3*r0/16

res[1][3] |> simplify |> println  # r4 北円の半径

   3*r0*(25 - 12*sqrt(3))/193

以上の結果から,問題への答えは,「北円の半径は南円の半径の (32√3 + 62)/193 ≒ 0.6084229318 倍」である。

(3*r0*(25 - 12*sqrt(Sym(3)))/193) / (3*r0*(7 - 4*sqrt(Sym(3)))/2) |> simplify |> println

   32*sqrt(3)/193 + 62/193

扇形の半径が 20 のとき,各円の直径を求める。

   南円の直径 = 4.307806
   北円の直径 = 2.620968
   東円の直径 = 10.000000
   西円の直径 = 7.500000
   扇形の半径 = 20.000000
   北円の直径は南円の直径の 0.6084229318 倍

一挙に数値解を求めることもできる。

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)
   (r0, r2, x2, x3, r4, x4, y4) = u
   return [
x2^2 + (r0/4 - r2)^2 - (r0/4 + r2)^2,  # eq1
x4^2 + (-3*r0/4 + y4)^2 - (r0/4 + r4)^2,  # eq2
-(r2 + r4)^2 + (x2 - x4)^2 + (-r0/2 - r2 + y4)^2,  # eq3
x3^2 + (r0/2 - r3)^2 - (r0/2 + r3)^2,  # eq4
x2^2 + (r0/2 + r2)^2 - (r0 - r2)^2,  # eq5
x4^2 + y4^2 - (r0 - r4)^2,  # eq6
r3/(sqrt(3)*r0/2 - x3) - 2 + sqrt(3),  # eq7
   ]
end;

r3 = 1
iniv = [big"157.0", 29.1, 68.2, 73.2, 10.3, 45.2, 140]
iniv = [big"9.0", 2, 4, 4, 0.5, 3, 8]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[9.285468820183671312972329900930808394772674258716870318199501414526902889954248, 1.74102540378443837117861943355628490397136248947194920156030216985869867358414, 4.020725942163689539353471755593597407534879935141037673297496272009488982864332, 4.309401076758502716985195569314300770820850779232595561630515863568051039690857, 0.6084229318248914707848103579780154432784296686583510598363298988001496181698688, 2.621938437530752623591216631288396179019582647806301923444247679455871388075398, 8.271430600475518861666281929830511546667832242439599342793407445063645710020953], true)

@printf("r3 = %.6f\n", r3)
names = ("r0", "r2", "x2", "x3", "r4", "x4", "y4")
for (i, name) in enumerate(names)
   @printf("%2s = %.6f\n", name, res[1][i])
end

   r3 = 1.000000
   r0 = 9.285469
   r2 = 1.741025
   x2 = 4.020726
   x3 = 4.309401
   r4 = 0.608423
   x4 = 2.621938
   y4 = 8.271431

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 20
   r1 = r0/4
   (r3, x3) = (3*r0*(7 - 4*sqrt(3))/2, r0*(-3 + 2*sqrt(3)))
   (r2, x2, r4, x4, y4) = (3*r0/16, r0*(-3 + 14*sqrt(3))*sqrt(28*sqrt(3) + 199)/772, 3*r0*(25 - 12*sqrt(3))/193, r0*sqrt(336*sqrt(3)/37249 + 2388/37249), r0*(68/193 + 60*sqrt(3)/193))
   @printf("南円の直径 = %.6f\n", 2r3)
   @printf("北円の直径 = %.6f\n", 2r4)
   @printf("東円の直径 = %.6f\n", 2r1)
   @printf("西円の直径 = %.6f\n", 2r2)
   @printf("扇形の半径 = %.6f\n", r0)
   @printf("北円の直径は南円の直径の %.10f 倍\n", (32√3 + 62)/193)
   plot()
   circle(0, 0, r0, :black, beginangle=30, endangle=150)
   circle(0, 0, r0/2, :black, beginangle=30, endangle=150)
   circle(x3, r0/2 - r3, r3, :blue)
   circle(-x3, r0/2 - r3, r3, :blue)
   circle(0, r0 - r1, r1)
   circle(x2, r0/2 + r2, r2, :green)
   circle(-x2, r0/2 + r2, r2, :green)
   circle(x4, y4, r4, :orange)
   circle(-x4, y4, r4, :orange)
   segment(-√3r0/2, r0/2, √3r0/2, r0/2)
   segment(-√3r0/2, r0/2, 0, 0)
   segment(√3r0/2, r0/2, 0, 0)
   if more
       point(0, r0 - r1, " 東円:r1\n r0-r1", :red)
       point(x2, r0/2 + r2, " 西円:r2\n(x2,r0/2+r2)", :green, :center)
       point(x3, r0/2 - r3, "   南円:r3\n(x3,r0/2-r3)", :blue, :left, :bottom)
       point(x4, y4, "  北円:r4,(x4,y4)\n", :orange, :left, :bottom)
       point(0, r0/2, " r0/2", :black)
       point(0, r0, " r0", :black, :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でシェアする

算額(その298)

2023年06月24日 | Julia

算額(その298)

愛媛県 伊佐爾波神社 昭和12年(1937)
算額22 中村正教

http://www.wasan.jp/ehime/isaniwa22.html
https://isaniwa.official.jp/%E5%BE%A1%E5%AE%9D%E7%89%A9/%E7%AE%97%E9%A1%8D/

大円の中に中円 4 個,小円 1 個が入っている。大円の径を知って,中円,小円の径を求めよ。

大円,中円,小円の半径を r0, r1, r2 とし,図に示すように記号を定め,以下の連立方程式を r1, r2 について解く(r0 は未知数として記号のまま残る)。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r0::positive

eq1 = 2r1^2 - (r1 + r2)^2
eq2 = 2r1^2 - (r0 - r1)^2

res = solve([eq1, eq2], (r1, r2))

   1-element Vector{Tuple{Sym, Sym}}:
    (-r0 + sqrt(2)*r0, r0*(3 - 2*sqrt(2)))

大円の半径を r0 とすると,中円の半径は (√2 - 1)r0,小円の半径は (3 - √8)r0 である。

r0 = 2 のとき,r1 = 0.82843;  r2 = 0.34315

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 2  # 大円の直径が 4 のとき
   (r1, r2) = ((√2 - 1)r0, (3 - √8)r0)
   @printf("r0 = %.5f; r1 = %.5f;  r2 = %.5f\n", r0, r1, r2)
   plot()
   circle(0, 0, r0, :black)
   circle4(r1, r1, r1, :blue)
   circle(0, 0, r2)
   if more
       point(r1, r1, " 中円:r1\n (r1,r1)", :blue)
       point(0, 0, " 小円:r2,(0,0)", :red, :center)
       point(r0, 0, "r0 ", :black, :right, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       zoomin && plot!(showaxis=true, xlims=(0, 47), ylims=(0, 35))
   else
       plot!(showaxis=false)
   end
end;

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

算額(その297)

2023年06月23日 | Julia

算額(その297)

中村信弥「算額への招待 追補2」
http://www.wasan.jp/syotai/syotai.html
長野県千曲市八幡 武水別神社(八幡神社) 天保8年(1837) 

5個の正方形と1個の円が直角三角形に内接している。円の直径が 7 寸,一番小さい正方形の一辺の長さが 6 寸のとき,鉤,股,弦の長さを求めよ。

(復元)算額の図が小さく(あるいは誤解?)分かりづらいが,円と一番小さい正方形の関係は,以下の拡大図のようになっている(算額では円が正方形に内接しているように描かれている)。


また,「算額への招待 追補2」は算額の術の解説をしており,上に示す 1 つの解のみを提示しているが,解は二通りある。もう一つの解の拡大図は以下のとおりである。

記号の説明のために全体図も二通り示しておく。
この 2 つの解は,図を裏返してみれば相似であることがわかる。

何れにせよ,SymPy では二通りの解を与える。


また,図に現れる直角三角形は全て相似なので,x3, y3 以降は SymPy で多元連立方程式を解くことでも得られるが,逐次的に繰り返して求めることができる。直角三角形の頂点 (a, b) と弦の長さ c も計算によって求める。

図に示すように記号を定め,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms x1::positive, x2::positive, y2::positive, r1::positive;

r1 = 7 // 2
x2 = x1 + 6
eq1 = (x2 - x1) / x1 - y2 / x2
eq2 = x2 + y2 - sqrt(x2^2 + y2^2) - 2r1
res = solve([eq1, eq2], (x1, y2))

   2-element Vector{Tuple{Sym, Sym}}:
    (9/2, 14)
    (8, 21/2)

● x1, y2 が (9/2, 14) の場合の解

x1 = 9/2;  x2 = 21/2;  y2 = 14
鉤 = 1023.63621399177
股 = 767.727160493827
弦 = 1279.54526748971

● x1, y2 が (8, 21/2) の場合の解

x1 = 8;  x2 = 14;  y2 = 21/2
鉤 = 182.185253906250
股 = 242.913671875000
弦 = 303.642089843750

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (s, r) = (6, 7/2)
   (x1, y2) = res[2]
   println((x1, y2))
   x2 = x1 + s
   tanθ = y2/x2
   sinθ = y2/sqrt(x2^2 + y2^2)
   cosθ = x2/sqrt(x2^2 + y2^2)
   xs = zeros(7)
   ys = zeros(7)
   xs[1:2] = [x1, x2]
   ys[1:2] = [s, y2]
   for i = 3:7
       xs[i] = xs[i - 1] + ys[i - 1]
       ys[i] = xs[i]*tanθ
   end
   println("ys[6]= $(ys[6])")
   c = xs[7] + ys[6] * tanθ
   a = xs[6] + ys[6]*cosθ^2
   b = ys[6] + ys[6]*cosθ*sinθ
   println("x1 = $x1;  x2 = $(x1 + 6);  y2 = $y2")
   println("鉤 = $(sqrt(b^2 + (c - a)^2))")
   println("股 = $(sqrt(b^2 + a^2))")
   println("弦 = $c")
   if zoomin
       plot()
   else
       plot([0, c, a, 0], [0, 0, b, 0], color=:black, lw=0.5)
   end
   circle(xs[2] - r, r, r)
   for i = 1:6
       println((xs[i], ys[i], xs[i+1], 0))
       rect(xs[i], 0, xs[i+1], ys[i], :blue)
   end
   abline(0, 0, tanθ, 0, 60)
   if more
       if !zoomin
           point(c, 0, " c", :black, :left, :bottom)
           point(a, b, "(a,b) ", :black, :right, :bottom)
       end
       point(x1, 0, "x1 ", :blue, :right, :bottom)
       point(x1, s, "(x1,s) ", :blue, :right, :bottom) 
       point(x2-r, r, "(x2-r,r)")
       point(x2, 0, " x2", :blue, :left, :bottom)
       point(x2, y2, "(x2,y2) ", :blue, :right, :bottom) 
       point(xs[3], 0, " x3", :blue, :left, :bottom)
       point(xs[3], ys[3], "(x3,y3) ", :blue, :right, :bottom) 
       if zoomin
           txt = @sprintf(" x1 = %.1f\n x2 = %.1f\n y2 = %.1f", x1, x1 + 6, y2)
           annotate!(xs[3], ys[3]/2, text(txt, 10, :left))
       else
           point(xs[4], 0, " x4", :blue, :left, :bottom)
           point(xs[4], ys[4], "(x4,y4) ", :blue, :right, :bottom) 
           point(xs[5], 0, " x5", :blue, :left, :bottom)
           point(xs[5], ys[5], "(x5,y5) ", :blue, :right, :bottom) 
           point(xs[6], 0, " x6", :blue, :left, :bottom)
           point(xs[6], ys[6], "(x6,y6) ", :blue, :right, :bottom) 
           point(xs[7], 0, " x7", :blue, :left, :bottom)
       end
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       zoomin && plot!(showaxis=true, xlims=(0, 47), ylims=(0, 35))
   else
       plot!(showaxis=false)
   end
end;

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

算額(その296)

2023年06月22日 | Julia

算額(その296)

中村信弥「算額への招待 追補2」
http://www.wasan.jp/syotai/syotai.html
長野県千曲市八幡 武水別神社(八幡神社) 天保8年(1837) 

天円,地円,人円が互いに外接し,三角形に内接している。天円と地円の直径がそれぞれ 8 寸,4.5 寸のとき,三角形の三辺の長さを求めよ。

図に示すように記号を定め,以下の連立方程式を nlsolve() により解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive,
     x1::positive, x2::positive, x3::positive, x32::positive,
     a::positive, b::positive, c::positive;

(r1, r2) = (8//2, 9//4)
eq1 = r3/x3 - r2/x2
eq2 = r3/x3 - r1/x1
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq4 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq5 = (x32 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq6 = (a + sqrt((a - b)^2 + c^2) + sqrt(b^2 + c^2)) * r1 - a * c
eq7 = r3 / (a - x32) - r1 / (a - x1)
eq8 = distance(b, c, a, 0, x1, r1) - r1^2

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])

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)
   (r3, x1, x2, x3, x32, a, b, c) = u
   return [
r3/x3 - 9/(4*x2),  # eq1
r3/x3 - 4/x1,  # eq2
(9/4 - r3)^2 - (r3 + 9/4)^2 + (x2 - x3)^2,  # eq3
(x1 - x2)^2 - 36,  # eq4
(4 - r3)^2 - (r3 + 4)^2 + (-x1 + x32)^2,  # eq5
-a*c + 4*a + 4*sqrt(b^2 + c^2) + 4*sqrt(c^2 + (a - b)^2),  # eq6
r3/(a - x32) - 4/(a - x1),  # eq7
(x1 - (a^2*x1 - 2*a*b*x1 + a*c^2 - 4*a*c + b^2*x1 + 4*b*c)/(a^2 - 2*a*b + b^2 + c^2))^2 + (-c*(a^2 - a*b - a*x1 + b*x1 + 4*c)/(a^2 - 2*a*b + b^2 + c^2) + 4)^2 - 16,  # eq8
      ]
end;

iniv = [big"1.3", 14, 8, 4, 18, 20, 15, 10]
res = nls(H, ini=iniv);

using Printf
names =  ("r3", "x1", "x2", "x3", "x32", "a", "b", "c")
for (i, name) in enumerate(names)
   @printf("%s = %.6f\n", name, res[1][i])
end
println(res[2])

   r3 = 1.265625
   x1 = 13.714286
   x2 = 7.714286
   x3 = 4.339286
   x32 = 18.214286
   a = 20.297143
   b = 15.250421
   c = 9.723228
   true

大斜 = a = 20.29714285714286
中斜 = sqrt(b^2 + c^2) = 18.08636238036625
小斜 = sqrt((a - b)^2 + c^2)) = 10.954933808937678

算額の答えは,大斜 = 41.427,中斜 = 29.533,小斜 = 15.534 であるが,これはありえない。図を描くとわかる。

中村信弥「算額への招待 追補2」では,大斜の式として,天地人の直径をそれぞれ a, x, b として
l = a/(a - b)*(sqrt(a*sqrt(a*b)) + sqrt(b*sqrt(a*b)) + sqrt(a*b))
を導いている。この式は正しいが,
「a = 8, b = 4.5 とすると l = 16*(6 + 7√3)/7 = 41.427097...」としているが,4.5 は地円の直径であり,人円の直径はこの時点では未知である。
実は,途中で x = sqrt(a*b) とあるので,既知の a, x から b を求めればよい。すなわち,b = x^2/a である。b = 4.5 を代入したために間違った解を求めてしまったのである。
上の式で,a = 8, b = 4.5^2/8 = 2.53125 を代入すれば,大斜 = 20.297142857142855 になる。元の算額の術も同じ間違いをおかしている。

a = 8
b = 4.5^2/8  # = 2.53125
l = a/(a - b)*(sqrt(a*sqrt(a*b)) + sqrt(b*sqrt(a*b)) + sqrt(a*b))

   20.297142857142855

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, x1, x2, x3, x32, a, b, c) = res[1]
   println("r1 = $(Float64(r1));  r2 = $(Float64(r2))")
   println("大斜 = $(Float64(a));  中斜 = $(Float64(sqrt(b^2 + c^2)));  小斜 = $(Float64(sqrt((a - b)^2 + c^2)))")
   plot([0, a, b, 0], [0, 0, c, 0], color=:black, lw=0.5)
   circle(x1, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   circle(x32, r3, r3, :green)
   if more
       point(x1, r1, "天:r1,(x1,r1)", :red, :center)
       point(x2, r2, "地:r2,(x2,r2)", :blue, :center)
       point(x3, r3, "人:r3,(x3,r3)", :green, :center)
       point(x32, r3, "人:r3,(x32,r3)", :green, :center)
       point(b, c, " (b,c)", :green, :left, :bottom)
       point(a, 0, " a", :green, :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でシェアする

算額(その295)

2023年06月22日 | Julia

算額(その295)

中村信弥「算額への招待 追補2」
http://www.wasan.jp/syotai/syotai.html
長野県千曲市八幡 武水別神社(八幡神社) 天保8年(1837) 

直角三角形(鉤股弦)の中に,天円,地円,人円それぞれ 1 個,小円 2 個が入っている。
天地人円の直径の和と大円の直径の積,および直角三角形の周の二乗の和が 46216 歩のとき,鉤および股を求めよ。

図に示すように記号を定め,以下の連立方程式を nlsolve() により解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, r4::positive, r5::positive,
     x2::positive, x3::positive, x4::positive, x5::positive,
     鉤::positive, 股::positive, 弦::positive;

eq1 = 2(r2 + r3 + r4) * 2r1 + (鉤 + 股 + 弦)^2 - 46216
eq2 = 鉤 + 股 - 弦 - 2r1
eq3 = 鉤^2 + 股^2 - 弦^2
eq4 = 2(r1 - r5)^2 - (r1 + r5)^2
eq5 = r1 / (股 - r1) - r2 / (股 -x2)
eq6 = r1 / (股 - r1) - r3 / (股 -x3)
eq7 = r1 / (股 - r1) - r4 / (股 -x4)
eq8 = r1 / (股 - r1) - r5 / (股 -x5)
eq9 = (r1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq10 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq11 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
eq12 = (x4 - x5)^2 + (r4 - r5)^2 - (r4 + r5)^2;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])

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, r3, r4, r5, x2, x3, x4, x5) = u
   return [
       2*r1*(2*r2 + 2*r3 + 2*r4) + (弦 + 股 + 鉤)^2 - 46216,  # eq1
       -2*r1 - 弦 + 股 + 鉤,  # eq2
       -弦^2 + 股^2 + 鉤^2,  # eq3
       2*(r1 - r5)^2 - (r1 + r5)^2,  # eq4
       r1/(-r1 + 股) - r2/(-x2 + 股),  # eq5
       r1/(-r1 + 股) - r3/(-x3 + 股),  # eq6
       r1/(-r1 + 股) - r4/(-x4 + 股),  # eq7
       r1/(-r1 + 股) - r5/(-x5 + 股),  # eq8
       (r1 - r2)^2 - (r1 + r2)^2 + (r1 - x2)^2,  # eq9
       (r2 - r3)^2 - (r2 + r3)^2 + (x2 - x3)^2,  # eq10
       (r3 - r4)^2 - (r3 + r4)^2 + (x3 - x4)^2,  # eq11
       (r4 - r5)^2 - (r4 + r5)^2 + (x4 - x5)^2,  # eq12
      ]
end;

iniv = [big"39.0", 83, 91, 15, 10, 6, 4, 3, 39, 54, 63, 70]
res = nls(H, ini=iniv);

using Printf
names =  ("鉤", "股", "弦", "r1", "r2", "r3", "r4", "r5", "x2", "x3", "x4", "x5")
for (i, name) in enumerate(names)
   @printf("%s = %.6f\n", name, res[1][i])
end

   鉤 = 38.566769
   股 = 82.527544
   弦 = 91.094408
   r1 = 14.999953
   r2 = 9.653883
   r3 = 6.213184
   r4 = 3.998769
   r5 = 2.573585
   x2 = 39.067174
   x3 = 54.556700
   x4 = 64.525669
   x5 = 70.941641

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, 弦, r1, r2, r3, r4, r5, x2, x3, x4, x5) = res[1] 
   plot([0, 股, 0, 0], [0, 0, 鉤, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   circle(x4, r4, r4, :magenta)
   circle(x5, r5, r5, :brown)
   circle(r5, r5, r5, :brown)
   if more
       point(r1, r1, "大:r1\n(r1,r1)", :red)
       point(x2, r2, "天:r2\n(x2,r2)", :blue)
       point(x3, r3, "地:r3\n(x3,r3)", :green)
       point(x4, r4, "  人:r4,(x4,r4)\n", :magenta, :left, :bottom)
       point(x5, r5, "    小:r5,(x5,r5)", :brown, :left, :bottom)
       point(r5, r5, "小:r5,(r5,r5)", :brown, :left, :top)
       point(0, 鉤, " 鉤", :green, :left, :bottom)
       point(股, 0, " 股", :green, :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でシェアする

算額(その294)

2023年06月22日 | Julia

算額(その294)

中村信弥「算額への招待 追補2」
http://www.wasan.jp/syotai/syotai.html
長野県千曲市八幡 武水別神社(八幡神社) 天保8年(1837) 

大円内に互いに外接し大円に内接する 3 個の小円がある。3 個の小円の面積と,小円に囲まれる面積(水色)の和に小円の直径をかけた積に小円の直径の3倍の 3 乗を加えると 794050 立方寸になる。このとき,大円の直径を求めよ。

大円の半径を r0,小円の半径,中心座標を r1, (r1, y), (0, r1) として,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, y::negative;
eq1 = r1^2 + (r0 - r1 - y)^2 - 4r1^2
eq2 = r1^2 + y^2 - (r0 - r1)^2
eq3 = (3PI*r1^2 + r1^2*sqrt(Sym(3)) - PI*r1^2/2)*2r1 + (6r1)^3 - 794050
res = solve([eq1, eq2, eq3], (r0, r1, y))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (sqrt(3)*794050^(1/3)*((20*sqrt(3)*pi + 25*pi^2 + 864*sqrt(3) + 2160*pi + 46668)^(1/3)*(450*pi^2 + 125*sqrt(3)*pi^3 + 38880*pi + 16200*sqrt(3)*pi^2 + 839880 + 700020*sqrt(3)*pi + 10085472*sqrt(3)) + 2*(20*sqrt(3)*pi + 25*pi^2 + 864*sqrt(3) + 2160*pi + 46668)^(11/6))/(3*(20*sqrt(3)*pi + 25*pi^2 + 864*sqrt(3) + 2160*pi + 46668)^2), 794050^(1/3)*(2*sqrt(3) + 5*pi + 216)/(20*sqrt(3)*pi + 25*pi^2 + 864*sqrt(3) + 2160*pi + 46668)^(2/3), -794050^(1/3)/(540*sqrt(3)*pi + 675*pi^2 + 23328*sqrt(3) + 58320*pi + 1260036)^(1/6))

術では円周率を 3.162 として求めているが,pi の場合についても計算しておく。

f(pi0) = 2*(sqrt(3)*794050^(1/3)*((20*sqrt(3)*pi0 + 25*pi0^2 + 864*sqrt(3) + 2160*pi0 + 46668)^(1/3)*(450*pi0^2 + 125*sqrt(3)*pi0^3 + 38880*pi0 + 16200*sqrt(3)*pi0^2 + 839880 + 700020*sqrt(3)*pi0 + 10085472*sqrt(3)) + 2*(20*sqrt(3)*pi0 + 25*pi0^2 + 864*sqrt(3) + 2160*pi0 + 46668)^(11/6))/(3*(20*sqrt(3)*pi0 + 25*pi0^2 + 864*sqrt(3) + 2160*pi0 + 46668)^2))
f(3.162) |> println
f(pi) |> println

   64.64101362909618
   64.65036111699513

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, y) = res[1]
   plot([r1, 0, -r1, r1], [y, r0 - r1, y, y], linecolor=:red, linewidth=0.5, seriestype=:shape, fillcolor=:cyan)
   circle(0, 0, r0, :black)
   circlef(0, r0 - r1, r1)
   circlef(r1, y, r1)
   circlef(-r1, y, r1)
   if more
       point(r1, y, "(r1,y)", :yellow)
       point(0, r0 - r1, " r0-r1", :yellow)
       point(r0, 0, " r0")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その293)

2023年06月21日 | Julia

算額(その293)

中村信弥「算額への招待 追補1」
http://www.wasan.jp/syotai/syotai.html
長野県飯山市滝澤家所蔵
http://www.wasan.jp/nagano/takizawa.html

長方形内に 3 個の円が内接している。円の直径が 5 寸のとき,長方形の短辺の長さを求めよ。

問では3個の円について述べており,右上,左上の2個の円についての記述はない。意図は良くわからないがここではこの円についても半径を求めておく。

長方形の短辺の長さを y とする。
半径,中心座標を以下のように決める。
円  r1, (r1, r1), (0, y - r1)
小円 r2, (r2, y - r2)

以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, y::positive;

eq1 = r1^2 + (y - 2r1)^2 - 4r1^2
eq2 = (2r1 - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
res = solve(eq1, y)
res |> println

   Sym[r1*(2 - sqrt(3)), r1*(sqrt(3) + 2)]

解が 2 つ求まるが,r1*(sqrt(3) + 2) が適解である。
すなわち短辺の長さは 5/2 * (sqrt(3) + 2) = 9.330127018922193 寸である。

5/2 * (sqrt(3) + 2)

   9.330127018922193

res2 = solve([eq1, eq2], (y, r2))

   4-element Vector{Tuple{Sym, Sym}}:
    (r1*(2 - sqrt(3)), 2*r1*(2 - sqrt(3)))
    (r1*(2 - sqrt(3)), 2*r1*(sqrt(3) + 2))
    (r1*(sqrt(3) + 2), 2*r1*(2 - sqrt(3)))
    (r1*(sqrt(3) + 2), 2*r1*(sqrt(3) + 2))

三番目の組が適解である。

   y = 9.33013;  r2 = 1.33975

小円の半径は 2*r1*(2 - sqrt(3)) = 1.33975 である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 5 / 2
   (y, r2) = (r1*(sqrt(3) + 2), 2*r1*(2 - sqrt(3)))
   @printf("y = %.5f;  r2 = %.5f\n", y, r2)
   plot([0, 4r1, 4r1, 0, 0], [0, 0, y, y, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(3r1, r1, r1)
   circle(2r1, y - r1, r1)
   circle(r2, y - r2, r2, :blue)
   circle(4r1 - r2, y - r2, r2, :blue)
   if more
       point(r1, r1, "(r1,r1)", :red)
       point(3r1, r1, "(3r1,r1)", :red)
       point(2r1, y - r1, "(0,y-r1)", :red)
       point(r2, y - r2, "(r2,y-r2)", :blue, :center)
       point(4r1-r2, y - r2, "(4r1-r2,y-r2)", :blue, :center)
       point(0, y, " y")
       point(4r1, 0, "4r1 ", :green, :right, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

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

算額(その292)

2023年06月21日 | Julia

算額(その292)

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(282)
長野県中野市片塩 大徳寺 天保9年(1838)

大円と小円は互いに外接し,また 2 本の平行線の上下に接している。斜線は大円の共通接線である。
大円と小円の直径がそれぞれ 100 寸,61 寸のとき,黒く塗られた部分の面積(歩)を求めよ。

原点と平行線の距離を y とする。
小円の半径と中心座標 r1, (0, r1 + y)
大円の半径と中心座標 r2, (x2, r2 + y)
共通接線の傾きは ±((r2 + y) / x2),共通接線と平行線の交点の x 座標は ±x である。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::positive, y::positive;

(r1, r2) = (61//2, 100//2)
eq1 = x2^2 + (r2 - r1)^2 - (r2 + r1)^2
eq2 = distance(0, r2 + y, x2, 0, x2, r2 + y) - r2^2
res = solve([eq1, eq2], (x2, y))

   1-element Vector{Tuple{Sym, Sym}}:
    (10*sqrt(61), -50 + 25*sqrt(61)/3)

   x2 = 78.10250;  y = 15.08541
   黒い部分の面積 = 4166.666666666663

黒の部分の面積は,⊿PAO - ⊿PBC = 2((r2 + y)\*x2 - r2\*x) = 4166.666666666663 = 4166+2/3 である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (61//2, 100//2)
   (x2, y) = (10*sqrt(61), -50 + 25*sqrt(61)/3)
   @printf("x2 = %.5f;  y = %.5f\n", x2, y)
   plot()
   circle4(x2, r2 + y, r2)
   circle(0, r1 + y, r1, :blue)
   circle(0, -r1 - y, r1, :blue)
   slope = (r2 + y) / x2
   abline(x2, r2 + y, slope, -140, 80, :black)
   abline(x2, r2 + y, -slope, -80, 140, :black)
   abline(-x2, -r2 - y, slope, -80, 140, :black)
   abline(-x2, -r2 - y, -slope, -140, 80, :black)
   segment(-140, y, 140, y, :black)
   segment(-140, -y, 140, -y, :black)
   if more
       point(x2, 0, "  x2")
       point(0, r2 + y, " r2+y")
       point(x2, r2 + y, " 大円:r2\n (x2,r2+y)")
       point(0, r1 + y, "小円:r1\n (0,r1+y)")
       x = r2/slope
       plot!([x2, x, -x, -x2, -x, x, x2], [0, y, y, 0, -y, -y, 0],
           linecolor=:black, linewidth=0.5, seriestype=:shape, fillcolor=:gray)
       point(0, y, " y", :yellow, :left, :top)
       point(r2/slope, y, "x ", :yellow, :right, :top)
       print(2((r2 + y)*x2 - r2*x))
       point(-x2, 0, "  A", :yellow)
       point(0, 0, "O ", :yellow, :right)
       point(-x, -y, "  B", :green, :left, :top)
       point(0, -y, " C", :green, :right, :top)
       point(0, -r2-y, "    P", :green, :left, :center)
       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アクセスランキング にほんブログ村