今日，別のプログラムを書いていて，初期値設定のときにエラー発生。なんで?？
以下はなんの問題もないですね。
julia> x = [1, -2, 3 + 3]
3-element Vector{Int64}:
1
-2
6
しかし，以下はエラーになりました。
julia> x = [1, -2, 3 +3]
ERROR: syntax: missing separator in array expression
Stacktrace:
[1] top-level scope
@ none:1
まあ，「それはエラーだよねえ」とも思いますが。どうも，しっくり来ません。か?？？
算額（その141）
福島県田村郡三春町御木沢 厳島神社 明治18年(1885)1月
http://www.wasan.jp/fukusima/miharuitukusima5.html
半円の中に大中小の正方形が入っている。それぞれの正方形の一辺の長さを求めよ。
半円の半径を1，大中小の正方形の一辺の長さを a, b, c とおく。方程式を立て，解く。
using SymPy
@syms a::positive, b::positive, c::positive;
eq2 = (1 - 2a)/a - c/(b - c); # 直角三角形が相似であることから
eq3 = (1 - b)/(a + b) - (2a - b)/b; # 直角三角形が相似であることから
eq5 = (a + b)^2 + b^2 - 1; # 右の中正方形の右上隅が円周上にあることから
res = solve([eq2, eq3, eq5], (a, b, c))
1-element Vector{Tuple{Sym, Sym, Sym}}:
(-2^(2/3)/5 + 1/5 + 2*2^(1/3)/5, -2/5 + 2^(1/3)/5 + 2*2^(2/3)/5, -4/5 + 3*2^(2/3)/10 + 2*2^(1/3)/5)
(-2^(2/3)/5 + 1/5 + 2*2^(1/3)/5, b-2/5 + 2^(1/3)/5 + 2*2^(2/3)/5, -4/5 + 3*2^(2/3)/10 + 2*2^(1/3)/5)
(0.3864882095643094, b + 0.486944630766254, 0.18018873554840903)
正方形の一辺は，大: 0.77298, 中: 0.48694, 小: 0.18019 である。
算額では，中の一辺が与えられたとき小の一辺はいくつになるかとあるので，中の一辺の 0.18018873554840903/0.4869446307662544 = 0.3700394750525633 倍が，小の一辺である。
算額の答えは「2.5 を解立法し（3乗根を取り），1を引いたものを中の一辺に掛ける」と読める。2.5^(1/3) - 1 = 0.35720880829745316 なので，あまりよい近似ではない。
using Plots
using Printf
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
if fill
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
else
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end
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 filledsquare(x1, y1, x2, y2, color=:black)
plot!([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], seriestype=:shape, fillcolor=color, color=color, lw=0.25)
end
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b, c) = res[1]
@printf("大: %.5f, 中: %.5f, 小: %.5f\n", 2a, b, c)
plot()
circle(0, 0, 1, :blue, beginangle=0, endangle=180, fill=true)
plot!([a + b, 0, -a - b, a + b], [b, 1, b, b], seriestype=:shape, fillcolor=:snow2, color=:snow2, lw=0.25)
filledsquare(-a, 0, a, 2a, :red)
filledsquare(a, 0, a + b, b, :navajowhite)
filledsquare(-a, 0, -a - b, b, :navajowhite)
filledsquare(a, b, a + c, b + c, :white)
filledsquare(-a, b, -a - c, b + c, :white)
segment(0, 1, a + b, b)
segment(0, 1, -a - b, b)
segment(-1, 0, 1, 0, :black)
if more == true
point(a, 0, "a")
point(0, 2a, "2a ", :black, :right)
point(a + c, b + c, "((a+c,b+c) ", :black, :right, :bottom)
point(a + b, 0, "a+b")
point(a + b, b, "(a+b,b) ", :black, :right)
point(0, a, "大", :black, :right)
point(a + b/2, b/2, "中", :black, :right)
point(a + c/2, b + c/2, "小", :black, :right)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5, ylims=(-0.075, 1))
else
plot!(showaxis=false)
end
end;
draw(false)
大: 0.77298, 中: 0.48694, 小: 0.18019
算額（その140）
福島県白河市明神 境明神 万延元年(1860)9月
http://www.wasan.jp/fukusima/sakai1.html
長方形の中に，半円 2 個，天円 4 個，地円 2 個，人円 2 個が入っている。人円の径が 1 寸のとき，天円の径はいくらか。
半円，天円，地円，人円の半径 をそれぞれ r0, r1, r2, r3 = 1 とおく。天円，地円の中心の x 座標を x1, x2 とおく。方程式を立て，解く。
using SymPy
@syms r0::positive, r1::positive, r2::positive, r3::positive, x1::positive, x2::positive;
r3 = 1
eq1 = x2^2 + r3^2 - (r2 + r3)^2 # 地円と人円が外接する
eq2 = x2^2 + (2r3 + r0)^2 - (r0 + r2)^2 # 半円と地円が外接する
eq3 = (x1 - x2)^2 + r1^2 - (r1 + r2)^2 # 天円と地円が外接する
eq4 = x1^2 + (2r3 + r0 - r1)^2 - (r0 + r1)^2 # 半円と天円が外接する
eq5 = 2r3 + r0 - 2r1; # 半径の関係
res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, r2, x1, x2))
1-element Vector{NTuple{5, Sym}}:
(6, 4, 14/5, 2*sqrt(21), 4*sqrt(21)/5)
天円の半径は 4，元の単位では天円の直径は 4 寸。
ちなみに半円の直径は 6 寸，地円の直径は 14/5 = 2寸8分。
using Plots
using Printf
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, lw=0.5, linestyle=:solid)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
plot!(ox .+ x, oy .+ y, color=color, linewidth=lw, linestyle=linestyle)
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")
(r0, r1, r2, x1, x2) = res[1]
@printf("半円の半径 = %.3f\n", r0)
@printf("天円の半径 = %.3f\n", r1)
@printf("地円の半径 = %.3f\n", r2)
plot([x1 + r1, x1 + r1, -x1 - r1, -x1 - r1, x1 + r1], [-2r1, 2r1, 2r1, -2r1, -2r1])
circle(0, 2r3 + r0, r0, beginangle=180, endangle=360)
circle(0, -2r3 - r0, r0, beginangle=0, endangle=180)
circle(x1, r1, r1, :blue)
circle(x1, -r1, r1, :blue)
circle(-x1, r1, r1, :blue)
circle(-x1, -r1, r1, :blue)
circle(x2, 0, r2, :green)
circle(-x2, 0, r2, :green)
circle(0, r3, r3, :orange)
circle(0, -r3, r3, :orange)
if more == true
point(0, 2r3 + r0, "2r3+r0 ", :red, :right)
point(0, r3, "r3", :orange, :right)
point(x2, 0, "x2", :green, :center)
point(x2+r1, 0, "x2+r2 ", :green, :right)
point(x1, 0, "x1", :blue, :center)
point(x1, r1, "(x1,r1)", :blue, :center)
point(x1 + r1, 0, "x1+r1", :blue, :right)
point(x1 + r1, 2r1, "(x1+r1,2r1) ", :blue, :right)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
draw(false)
半円の半径 = 6.000
天円の半径 = 4.000
地円の半径 = 2.800
算額（その139）
大阪府茨木市 井於神社 弘化3年(1846)
http://www.wasan.jp/osaka/iyo.html
鉤，股，弦がそれぞれ 3寸1分，3寸6分，4寸8分である三角形が2本の斜線で区切られ，径の同じ円が 2 個入っている。円の径を求めよ。
鉤股弦とはいっているが，直角三角形ではない。鉤と股をそのままの数値として受け入れれば，弦は 4.750789408087881 寸でなければならない。
鉤 = 31, 股 = 36 円の半径を r，それぞれの中心座標を (r, y), (x, r) とおき，直線までの距離に関する方程式を立て，解く。
solve() では解けないので，nlsolve() で数値解を求める。
using SymPy
function distance(x1, y1, x2, y2, x0, y0)
l = sympy.Line(sympy.Point(x1, y1), sympy.Point(x2, y2))
l.distance(sympy.Point(x0, y0))^2 # 二乗距離を返す！！
end;
@syms r::positive, x::positive, y::positive, x0::positive, y0::positive;
eq1 = distance(0, 31, x0, y0, r, y) - r^2 # 円1の中心から BD への距離
eq2 = distance(0, 31, x0, y0, x, r) - r^2 # 円2の中心から BD への距離
eq3 = distance(0, 31, 36, 0, x, r) - r^2 # 円2の中心から AB への距離
eq4 = distance(0, 0, x0, y0, r, y) - r^2 # 円1の中心から OC への距離
eq5 = distance(0, 0, x0, y0, x, r) - r^2; # 円2の中心から OC への距離
# res = solve([eq1, eq2, eq3, eq4, eq5], (r, x, y, x0, y0))
println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")
println(eq5, ",")
-r^2 + (r - x0*(r*x0 + y*y0 - 31*y - 31*y0 + 961)/(x0^2 + y0^2 - 62*y0 + 961))^2 + (y - (r*x0*y0 - 31*r*x0 + 31*x0^2 + y*y0^2 - 62*y*y0 + 961*y)/(x0^2 + y0^2 - 62*y0 + 961))^2,
-r^2 + (r - (r*y0^2 - 62*r*y0 + 961*r + x*x0*y0 - 31*x*x0 + 31*x0^2)/(x0^2 + y0^2 - 62*y0 + 961))^2 + (x - x0*(r*y0 - 31*r + x*x0 - 31*y0 + 961)/(x0^2 + y0^2 - 62*y0 + 961))^2,
-r^2 + (1116*r/2257 + 961*x/2257 - 34596/2257)^2 + (1296*r/2257 + 1116*x/2257 - 40176/2257)^2,
-r^2 + (r - x0*(r*x0 + y*y0)/(x0^2 + y0^2))^2 + (y - y0*(r*x0 + y*y0)/(x0^2 + y0^2))^2,
-r^2 + (r - y0*(r*y0 + x*x0)/(x0^2 + y0^2))^2 + (x - x0*(r*y0 + x*x0)/(x0^2 + y0^2))^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=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)
(r, x, y, x0, y0) = u
return [
-r^2 + (r - x0*(r*x0 + y*y0 - 31*y - 31*y0 + 961)/(x0^2 + y0^2 - 62*y0 + 961))^2 + (y - (r*x0*y0 - 31*r*x0 + 31*x0^2 + y*y0^2 - 62*y*y0 + 961*y)/(x0^2 + y0^2 - 62*y0 + 961))^2,
-r^2 + (r - (r*y0^2 - 62*r*y0 + 961*r + x*x0*y0 - 31*x*x0 + 31*x0^2)/(x0^2 + y0^2 - 62*y0 + 961))^2 + (x - x0*(r*y0 - 31*r + x*x0 - 31*y0 + 961)/(x0^2 + y0^2 - 62*y0 + 961))^2,
-r^2 + (1116*r/2257 + 961*x/2257 - 34596/2257)^2 + (1296*r/2257 + 1116*x/2257 - 40176/2257)^2,
-r^2 + (r - x0*(r*x0 + y*y0)/(x0^2 + y0^2))^2 + (y - y0*(r*x0 + y*y0)/(x0^2 + y0^2))^2,
-r^2 + (r - y0*(r*y0 + x*x0)/(x0^2 + y0^2))^2 + (x - x0*(r*y0 + x*x0)/(x0^2 + y0^2))^2,
]
end;
iniv = [5.0, 20, 10, 12, 10]
res = nls(H, ini=iniv)
println(res)
([5.6146321148429355, 20.875286969374045, 9.74605295956063, 13.244959542108463, 7.680342537201841], true)
直径は 2 × 5.6146321148429355 ≒ 11寸2分2厘9毛である。算額では 9分4厘5毛になっている。
～井於神社奉納「算額」について～ http://io-jinja.org/sangakusetumei.html に引用されている，「算術稽古大全」宅間流五世松岡能一，文化5年(1808) に以下が示されているとあるが，これはまた以上の 2 通りの解とも違う。
(鉤, 股, 弦) = (3.1, 3.6, 4.8)
天 = 鉤 + 股 + 弦
等円径 = (sqrt(2天 * 鉤) + 天) / 2(鉤*股)
0.8935453733438091
弦に正しい（?）数値を代入しても異なる解になる。
(鉤, 股, 弦) = (3.1, 3.6, 4.750789408087881)
天 = 鉤 + 股 + 弦
等円径 = (sqrt(2天 * 鉤) + 天) / 2(鉤*股)
0.8905302961492191
using Plots
using Printf
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, lw=0.5, linestyle=:solid)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
plot!(ox .+ x, oy .+ y, color=color, linewidth=lw, linestyle=linestyle)
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 intersection(x1, y1, x2, y2, x3, y3, x4, y4)
denominator = (x1*y3 - x1*y4 - x2*y3 + x2*y4 - x3*y1 + x3*y2 + x4*y1 - x4*y2)
X = (x1*x3*y2 - x1*x3*y4 - x1*x4*y2 + x1*x4*y3 - x2*x3*y1 + x2*x3*y4 + x2*x4*y1 - x2*x4*y3) / denominator
Y = (x1*y2*y3 - x1*y2*y4 - x2*y1*y3 + x2*y1*y4 - x3*y1*y4 + x3*y2*y4 + x4*y1*y3 - x4*y2*y3) / denominator
(X, Y)
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, x, y, x0, y0) = res[1]
@printf("円の直径 = %.3f\n", 2r)
plot([0, 36, 0, 0], [0, 0, 31, 0])
circle(r, y, r)
circle(x, r, r)
(X1, Y1) = intersection(0, 0, x0, y0, 0, 31, 36, 0)
segment(0, 0, X1, Y1)
(X2, Y2) = intersection(0, 31, x0, y0, 0, 0, 36, 0)
segment(0, 31, X2, Y2)
if more == true
point(r, y, "円1(r,y;r)", :red, :center, :top)
point(x, r, "円2(x,r;r)", :red, :center, :top)
point(x0, y0, "(x0,y0)", :red)
point(0, 0, "O ", :blue, :right, :bottom)
point(36, 0, " A", :blue, :left, :bottom)
point(0, 31, " B", :blue, :left, :bottom)
point(X1, Y1, " C", :blue, :left, :bottom)
point(X2, Y2, "D ", :blue, :right, :bottom)
hline!([0], color=:black, lw=0.5, xlims=(-2, 38))
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
draw(false)
円の直径 = 11.229
