二年ぶりにたくさん蕾がついた。嬉しさと，寂しさ。ないまぜ。
https://www.543life.com/72seasons/hangesyouzu.html
半夏生 ハンゲショウズ と読むのは初めて知った（？）
田植えが終わったお疲れ様パーティーで，地方によっては団子やタコを食べる慣わしがあるとか。先日も，スーパーへ行ったらタコが推しだったとか聞いた。
植物の世界では，文字通り「ハンゲショウ」というのがある。「ランマン」で出たかな?（見てないので未確認）
不具合に関わらず，マイナンバーカード関連の政府発の記事で，必ず，イラストに女性のイラストが付いたマイナンバーカードのイメージ図がつくが，あれは，逆差別か？
そもそも，イメージ図なんぞ不要と思うし。記事は大抵の場合，マイナスイメージのことが多いので，女性のイラストを使うのはどうだかなあと思う次第。
たまには，髭面の汚いおじさんの図でも使うと，親しみが湧くな(んなわけない)。
算額（その312）
「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県菰野町 広幡神社 嘉永5年(1852)
問題
半円の中に円弧があり，半円と円弧に挟まれた領域にいくつかの円が挟まれている。半円の直径，丙円，丁円の直径が与えられたとき，戊円の直径はいくつか。
円弧の半径と中心座標を r0, (0, y0)
半円の半径と中心座標を r, (0, 0)
丙円の半径と中心座標を r1, (x1, y1)
丁円の半径と中心座標を r2, (x2, y2)
戊円の半径と中心座標を r3, (x3, y3)
として，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms r0::positive, y0::positive, r::positive,
r1::positive, x1::positive, y1::negative,
r2::positive, x2::positive, y2::negative,
r3::positive, x3::positive, y3::negative;
(r, r1, r2) = (50, 10, 5)
eq1 = x1 ^2 + (y0 - y1)^2 - (r0 + r1)^2
eq2 = x2 ^2 + (y0 - y2)^2 - (r0 + r2)^2
eq3 = x3 ^2 + (y0 - y3)^2 - (r0 + r3)^2
eq4 = x1^2 + y1^2 - (r - r1)^2
eq5 = x2^2 + y2^2 - (r - r2)^2
eq6 = x3^2 + y3^2 - (r - r3)^2
eq7 = (x2 - x1)^2 + (y2 - y1)^2 - (r2 + r1)^2
eq8 = (x3 - x2)^2 + (y3 - y2)^2 - (r3 + r2)^2
eq9 = r^2 + y0^2 - r0^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])
連立方程式が複雑すぎて solve() では解けないので，nlsolve() で数値解を求める。
println(eq1, ", # eq1")
println(eq2, ", # eq2")
println(eq3, ", # eq3")
println(eq4, ", # eq4")
println(eq5, ", # eq5")
println(eq6, ", # eq6")
println(eq7, ", # eq7")
println(eq8, ", # eq8")
println(eq9, ", # eq9")
x1^2 - (r0 + 10)^2 + (y0 - y1)^2, # eq1
x2^2 - (r0 + 5)^2 + (y0 - y2)^2, # eq2
x3^2 - (r0 + r3)^2 + (y0 - y3)^2, # eq3
x1^2 + y1^2 - 1600, # eq4
x2^2 + y2^2 - 2025, # eq5
x3^2 + y3^2 - (50 - r3)^2, # eq6
(-x1 + x2)^2 + (-y1 + y2)^2 - 225, # eq7
-(r3 + 5)^2 + (-x2 + x3)^2 + (-y2 + y3)^2, # eq8
-r0^2 + y0^2 + 2500, # 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)
(r0, y0, x1, y1, x2, y2, r3, x3, y3) = u
return [
x1^2 - (r0 + 10)^2 + (y0 - y1)^2, # eq1
x2^2 - (r0 + 5)^2 + (y0 - y2)^2, # eq2
x3^2 - (r0 + r3)^2 + (y0 - y3)^2, # eq3
x1^2 + y1^2 - 1600, # eq4
x2^2 + y2^2 - 2025, # eq5
x3^2 + y3^2 - (50 - r3)^2, # eq6
(-x1 + x2)^2 + (-y1 + y2)^2 - 225, # eq7
-(r3 + 5)^2 + (-x2 + x3)^2 + (-y2 + y3)^2, # eq8
-r0^2 + y0^2 + 2500, # eq9
]
end;
iniv = [big"66.0", 45.0, 31.2, -24.5, 42, -13, 3, 47, -7]
res = nls(H, ini=iniv);
println(res);
(BigFloat[76.12612612612612612612612611757250855206487617275733389169965733209449470577765, 57.40372007954591399580886043401747772776331606928505070760935723414311831179716, 33.42516087186933536638372029599284981641171160639776150802648249557181481914727, -21.97176872010205673632683968491622029305014289094789997653090282550024009759404, 43.6384044716071878394454126067690155978138981377281139789361700897187230488509, -10.9858843600510283681634198397671466786061193950162591863412958011058978766673, 2.226463104325699745547073791346257958834702049836892015992527235142629533435189, 47.52241383500506650373257942562149702505122375947591050606909041806671507003448, -4.8919332392084731919302760073024864943304938091114520230459338837681434939186], true)
円弧を構成する円の直径 = 152.252252, 中心座標 = (0.000000, 57.403720)
半円の直径 = 100.000000, 中心座標 = (0.000000, 0.000000)
円弧を構成する円の直径 = 152.252252, 中心座標 = (0.000000, 57.403720)
丙円の直径 = 20.000000, 中心座標 = (33.425161, -21.971769)
丁円の直径 = 10.000000, 中心座標 = (43.638404, -10.985884)
戊円の直径 = 4.452926, 中心座標 = (47.522414, -4.891933)
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r, r1, r2) = (50, 10, 5)
(r0, y0, x1, y1, x2, y2, r3, x3, y3) = res[1]
@printf("円弧を構成する円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r0, 0, y0)
@printf("半円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r, 0, 0)
@printf("円弧を構成する円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r0, 0, y0)
@printf("丙円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r1, x1, y1)
@printf("丁円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r2, x2, y2)
@printf("戊円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r3, x3, y3)
plot()
circle(0, y0, r0, :black, beginangle=180, endangle=360)
circle(0, 0, r, beginangle=180, endangle=360)
circle(x1, y1, r1, :blue)
circle(x2, y2, r2, :orange)
circle(x3, y3, r3, :magenta)
if more
point(0, y0, " y0", :black)
point(r, 0, "r ", :red, :right, :bottom)
point(x1, y1, " 丙円:r1,(x1,y1)", :blue, :left, :bottom)
point(x2, y2, " 丁円:r2,(x2,y2)", :orange, :left, :bottom)
point(x3, y3, " 戊円:r3,(x3,y3)", :magenta, :left, :bottom)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
算額（その311）
「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県菰野町 伎留太神社 寛政9年(1797)
問題文1
扇形の中に直径の等しい 3 個の円が入っている。条件として，弦（扇形を形成する左右の辺の端を結んだ線分），矢（円弧の中点から弦におろした垂線），左右斜（扇形の半径から扇形と中心が一致する内部の円の半径を引いた線分）が与えられているとき円の直径を求めよ。
扇の中心から扇の先端までの距離（扇の外円の半径）を R，等円の半径を r，右下の等円の中心座標を(r, R - 矢) として，以下の方程式を解く。
include("julia-source.txt");
using SymPy
@syms R::positive, r::positive, y1::positive,
左右斜::negative, 矢::positive, 弦::positive;
(左右斜, 矢, 弦) = (50, 40, 130)
eq1 = r^2 + (R - r - y1)^2 - 4r^2
eq2 = r^2 + y1^2 - (r + (R - 左右斜))^2
eq3 = (弦/2)^2 - (R^2 - (R - 矢)^2);
eq1 を y1 について解く。
solve(eq1, y1) |> println
Sym[R - r + sqrt(3)*r, R - sqrt(3)*r - r]
eq2 の y1 に代入する。
eq22 = eq2(y1 => R - r - sqrt(Sym(3))r) |> expand
eq22 |> println
-4*R*r - 2*sqrt(3)*R*r + 2*R*左右斜 + 2*sqrt(3)*r^2 + 4*r^2 + 2*r*左右斜 - 左右斜^2
eq22 を R について解く。
solve(eq22, R)[1] |> println
(sqrt(3)*r^2 + 2*r^2 + r*左右斜 - 左右斜^2/2)/(sqrt(3)*r + 2*r - 左右斜)
eq3 に代入する。
eq32 = eq3(R => (sqrt(Sym(3))*r^2 + 2*r^2 + r*左右斜 - 左右斜^2/2)/(sqrt(Sym(Sym(3)))*r + 2*r - 左右斜))
eq32 |> println
弦^2/4 + (-矢 + (sqrt(3)*r^2 + 2*r^2 + r*左右斜 - 左右斜^2/2)/(sqrt(3)*r + 2*r - 左右斜))^2 - (sqrt(3)*r^2 + 2*r^2 + r*左右斜 - 左右斜^2/2)^2/(sqrt(3)*r + 2*r - 左右斜)^2
左右斜, 矢, 弦 を代入すると r のみの式になる。
eq33 = eq32(左右斜 => 50, 矢 => 40, 弦 => 130)
(-40 + (sqrt(3)*r^2 + 2*r^2 + 50*r - 1250)/(sqrt(3)*r + 2*r - 50))^2 + 4225 - (sqrt(3)*r^2 + 2*r^2 + 50*r - 1250)^2/(sqrt(3)*r + 2*r - 50)^2
二通りの解（r）が得られる。
res = solve(eq33);
length(res)
2
いずれも虚数として表示されるが，虚部は誤差範囲で 0 である。
したがって，r は 45.2629281760725，14.1521122023713 であるが，後者が適解である。
solve(eq33)[1].evalf(), solve(eq33)[2].evalf()
(45.2629281760725 - 0.e-22*I, 14.1521122023713 + 0.e-21*I)
このあとさらに y1, R も求めなければならない。
---
最初から連立方程式に，左右斜, 矢, 弦 を代入して解を求めれば，R, r, y1 が一度に求まる。
using SymPy
@syms R::positive, r::positive, y1::positive,
左右斜::negative, 矢::positive, 弦::positive;
(左右斜, 矢, 弦) = (50, 40, 130)
eq1 = r^2 + (R - r - y1)^2 - 4r^2
eq2 = r^2 + y1^2 - (r + (R - 左右斜))^2
eq3 = (弦/2)^2 - (R^2 - (R - 矢)^2);
res = solve([eq1, eq2, eq3], (R, r, y1))
1-element Vector{Tuple{Sym, Sym, Sym}}:
(72.8125000000000, 14.1521122023713, 34.1482104287061)
左右斜= 50; 矢 = 40; 弦 = 130
R = 72.812500; r = 14.152112; y1 = 34.148210
等円の直径は 2r = 28.304224
等円の直径は 2×14.1521122023713 = 28.3042244047426 である。
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r, y1) = (72.8125000000000, 14.1521122023713, 34.1482104287061)
@printf("左右斜= %d; 矢 = %d; 弦 = %d\n", 左右斜, 矢, 弦)
@printf("R = %.6f; r = %.6f; y1 = %.6f\n", R, r, y1)
@printf("等円の直径は 2r = %.6f\n", 2r)
y = R - 矢
x = sqrt(R^2 - y^2)
θ = round(Int, atand(y/x))
plot()
circle(0, 0, R, beginangle=θ, endangle=180-θ)
circle(0, 0, R - 左右斜, :blue, beginangle=θ, endangle=180-θ)
circle(0, R - r, r, :green)
circle(r, y1, r, :green)
circle(-r, y1, r, :green)
plot!([-x, 0, x, -x], [y, 0, y, y], color=:black, lw=0.5)
if more
point(0, R - r, " R-r", :green)
point(0, y, " R-矢 ", :black, :right, :bottom)
point(r, y1, " (r,y1)", :green, :left, :bottom)
point(sqrt(R^2 - y^2), y, "(√(R^2-(R-矢)^2),R-矢) ", :black, :right)
point(0, R, " R", :red, :left, :bottom)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
PVアクセスランキング にほんブログ村
