Julia/SymPy: 数値代入法による係数の決定
https://blog.goo.ne.jp/r-de-r/e/4076953bde76502ffa5504a5d415c99d
であるが，一番簡単に解くには以下のようにする。
using SymPy
@syms a, b, c, d, x
eq = 9/(x-1)^2/(x+2)^2 - a/(x-1) - b/(x-1)^2 - c/(x+2) - d/(x+2)^2
solve(eq, a, b, c, d)
これで，解が以下のように示される。
Dict{Any, Any} with 4 entries:
b => 1
d => 1
c => 2/3
a => -2/3
算額（その1）
算額--2022年
https://blog.goo.ne.jp/r-de-r/e/12c8f8f08dbedbd16d5b8f54e32a17a8
が解けた。自力でやったのではなく，SymPy を使って。それでも相当骨が折れた。
三重県三重郡三重郡菰野町 廣幡神社 文化9年(1812)
http://www.wasan.jp/mie/hirohata2.html
山形県鶴尾市大山 椙尾神社 文政元年(1818年8月)
http://www.wasan.jp/yamagata/sugio.html
一〇一 大宮市高鼻町 氷川神社 明治31年(1898)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
岩手県一関市花泉町 大門神社・大門観世音菩薩 明治33年(1900)
http://www.wasan.jp/iwate/daimon5.html
答えは，「外側の円の半径を 1 とすると，大円，中円，小円の半径はそれぞれ (9 + √17)/32，(23 - √17)/64，(5 - √17)/4 である。」
三種類の円を含む外側の円が，原点を中心とする半径 1 の円とする。
図と記号の対応を示す。
using SymPy
@syms x, y
@syms α::real, β::real, γ::real # 大円，中円，小円の半径
@syms φ::real # 上の大円の中心の y 座標
@syms χ::real; # 右上の中円の中心の x 座標
## 1. それぞれの円の方程式
αeq = x^2 + (y - φ)^2 - α^2 # 大円の方程式
βeq = (x - χ)^2 + (y - β)^2 - β^2 # 中円の方程式
γeq = (x - γ)^2 + y^2 - γ^2 # 小円の方程式
## 2. 通過点による制約
eq1 = αeq(x => 0, y => 1) # (0, 1)
eq2 = αeq(x => 0, y => 1-2α) # (0, 1-2*α)
## 3. 重心間の距離による制約
eq5 = χ^2 + (1 - α - β)^2 - (α + β)^2 # α to β
eq6 = (χ - γ)^2 + β^2 - (γ + β)^2 # γ to β
eq7 = γ^2 + (1 - α)^2 - (α + γ)^2 # α to γ
## 4. 大円と中円の接点で成り立つこと
eq8 = sqrt(χ^2 + β^2) + β - 1
answer = solve([eq1, eq2, eq5, eq6, eq7, eq8], (φ, α, χ, β, γ))
3-element Vector{NTuple{5, Sym}}:
(0, 1, -1, 0, -1/2)
(23/32 - sqrt(17)/32, sqrt(17)/32 + 9/32, sqrt(2)/sqrt(9 - sqrt(17)), 23/64 - sqrt(17)/64, 5/4 - sqrt(17)/4)
(sqrt(17)/32 + 23/32, 9/32 - sqrt(17)/32, -sqrt(2)/sqrt(sqrt(17) + 9), sqrt(17)/64 + 23/64, sqrt(17)/4 + 5/4)
二組の解が得られるが，最初のものは明らかに不適である。
φ2, α2, χ2, β2, γ2 = [answer[2][i] for i in 1:5]
eq1, eq2 を検算してみる。
eq1(α => α2, φ => φ2).evalf() |> println
0
eq2(α => α2, φ => φ2).evalf() |> println
-0.e-124
eq5, eq6, eq7 を検算してみる。
eq5(χ => χ2, α => α2, β => β2).evalf() |> println
eq6(β => β2, γ => γ2, χ => χ2).evalf() |> println
eq7(γ => γ2, α => α2).evalf() |> println
-0.e-125
-0.e-125
-0.e-125
eq8 を検算してみる。
eq8(χ => χ2, β => β2).evalf() |> println
-0.e-125
## 5. 大円の半径
α2 |> println
sqrt(17)/32 + 9/32
α2 |> println
sqrt(17)/32 + 9/32
α2 |> N
0.410097050800551892181919057999189907035849975792925638574957299159217323323052
## 6. 中円の半径
β2 |> println
23/64 - sqrt(17)/64
β2 |> N
0.2949514745997240539090404710004050464820750121035371807125213504203913383384761
## 7. 小円の半径
γ2 |> println
5/4 - sqrt(17)/4
γ2 |> N
0.2192235935955848625446475360064807437132001936565948914003416067262614134155944
## 8. 大円の中心の y 座標
φ2 |> println
23/32 - sqrt(17)/32
φ2 |> N
0.5899029491994481078180809420008100929641500242070743614250427008407826766769523
## 9. 中円の中心の x 座標
χ2 |> println
sqrt(2)/sqrt(9 - sqrt(17))
χ2 |> N
0.6403882032022075687276762319967596281433999031717025542998291966368692932921995
二重根号を外して，簡約化する。
χ22 = sympy.sqrtdenest(χ2) |> simplify |> factor
χ22 |> println
(1 + sqrt(17))/8
χ22 |> N
0.6403882032022075687276762319967596281433999031717025542998291966368692932921995
## 10. 確認のため図を描いてみる
using Plots
function circle2(ox, oy, r, color)
theta = 0:0.01:2pi
x = r.*cos.(theta)
y = r.*sin.(theta)
plot!(ox .+ x, oy .+ y, color=color)
end;
gr(size=(500, 500), aspectratio=1, label="")
plot()
circle2(0, φ2, α2, :red) # 大円
circle2(0, -φ2, α2, :red)
circle2(χ2, β2, β2, :green) # 中円
circle2(χ2, -β2, β2, :green)
circle2(-χ2, β2, β2, :green)
circle2(-χ2, -β2, β2, :green)
circle2(γ2, 0.0, γ2, :black) # 小円
circle2(-γ2, 0.0, γ2, :black)
circle2(0, 0, 1, :blue); # 外円
savefig("fig1.png");
冒頭に示した，図と記号の対応を示すもの
3 つの円と，7 つのパラメータ
(cax, cay) = (0, convert(Float64, φ2))
(cbx, cby) = (convert(Float64, χ2), convert(Float64, β2))
(ccx, ccy) = (convert(Float64, γ2), 0)
(ra, rb, rc) = (convert(Float64, α2), convert(Float64, β2), convert(Float64, γ2))
annotate!(cax, cay, text("(0, φ)", 10, :center, :top))
annotate!(cbx,cby, text("(χ, β)", 10, :center, :top))
annotate!(ccx, ccy, text("γ, 0)", 10, :center, :top));
sx = [0.0, convert(Float64, χ2), convert(Float64, γ2)]
sy = [convert(Float64, φ2), convert(Float64, β2), 0.0]
scatter!(sx,sy);
ex = [cax, cbx, ccx] .+ [ra, rb, rc]
ey = [convert(Float64, φ2), convert(Float64, β2), 0.0]
scatter!(ex, ey);
for i in 1:3
annotate!((sx[i] + ex[i])/2 , (sy[i] + ey[i])/2, text(["α", "β", "γ"][i], 10, :bottom))
plot!([sx[i], ex[i]], [sy[i], ey[i]])
end
plot!()
circle2(-γ2, 0.0, γ2, :black)
circle2(0, 0, 1, :blue) # 外円
savefig("fig0.png");
整式を整式で割ったときの商と剰余
using SymPy
@syms x
(x,)
1. 商 quo()
quo(f, g, *gens, **args)
Compute polynomial quotient of f and g.
quo はエクスポートされていないので，sympy.quo() として使用する
例 x^2+1 を 2x - 4 で割ったときの商を求めよ。
sympy.quo(x^2 + 1, 2x - 4) |> println
x/2 + 1
2. 剰余 rem()
rem(f, g, *gens, **args)
Compute polynomial remainder of f and g.
rem はエクスポートされていないので，sympy.rem() として使用する
例 x^2+1 を 2x - 4 で割ったときの剰余を求めよ。
sympy.rem(x^2 + 1, 2x - 4) |> println
5
3. 商と剰余 div()
div(f, g, *gens, **args)
Compute polynomial division of f and g.
div はエクスポートされていないので，sympy.div() として使用する
例 x^2+1 を 2x - 4 で割ったときの商と剰余を求めよ。
sympy.div(x^2 + 1, 2x - 4) |> println
(x/2 + 1, 5)
検算
(2x - 4) * (x/2 + 1) + 5 |> expand |> println
x^2 + 1
4. 応用例
x^4 + a * x^2 - 2 * x + 3 を x^2 + x + b で割ったとき余りが 0 になるような実数の組 (a, b) をすべて求めよ。
@syms a, b, x
remainder = sympy.rem(x^4 + a * x^2 - 2 * x + 3, x^2 + x + b)
remainder |> println
-a*b - a*x + b^2 + 2*b*x - b - 3*x + 3
(remainder.coeff(x, 1), remainder.coeff(x, 0)) |> println # 若干わかりにくい
(-a + 2*b - 3, -a*b + b^2 - b + 3)
(remainder.coeff(x), remainder(x => 0)) |> println # Julia 流
(-a + 2*b - 3, -a*b + b^2 - b + 3)
solve([remainder.coeff(x), remainder(x => 0)], (a, b)) |> println
Tuple{Sym, Sym}[(-5, -1), (3, 3)]
三点を通る円の方程式を求める
以下を参照した
3点を通る円「ますいしいからの挑戦状(16)!!」をwolframAlphaとsympyでやってみた。
https://qiita.com/mrrclb48z/items/62d3502529d09b0ef843
注: Julia の SymPy では，Point(), Circle(), Point2D() ともにエクスポートされていないので，sympy. をつけて使用しなければならない。
using SymPy
@syms x y a b c
A = sympy.Point(8, 5)
B = sympy.Point(1,-2)
C = sympy.Point(9, 2)
eq = x^2+y^2+a*x+b*y+c
result = solve([eq(x => A.x, y => A.y),
eq(x => B.x, y => B.y),
eq(x => C.x, y => C.y)],[a,b,c])
println(result) # 解 a, b, c の値
Dict{Any, Any}(b => -4, c => -5, a => -8)
円の方程式は以下の通り。
eq2 = eq(a => result[a], b => result[b], c => result[c])
println(eq2) # eq の a, b, c に解を代入
x^2 - 8*x + y^2 - 4*y - 5
ans2 = sympy.Circle(eq2)
println(ans2) # 円の中心座標と半径
Circle(Point2D(4, 2), 5)
以下に示すのはものすごく簡単（SymPy は大変な計算をしているのだろうか）
eq3 = sympy.Circle(A, B, C).equation().expand()
println(eq3) # 「三点を通る方程式を展開する」という単純な命令で得られる
x^2 - 8*x + y^2 - 4*y - 5
陰関数 f(x, y) = 0 となる曲線を描くのに，ImplicitEquations パッケージを使う。
注: f ⩵ 0 は Eq(f, 0) と書くのと同じであるが，トップレベルでは ⩵ も Eq もそのままでは使えない ImplicitEquations.Eq としないとエラーになる。
using Plots
using ImplicitEquations
f(x,y) = x^2 - 8*x + y^2 - 4*y - 5
plot(ImplicitEquations.Eq(f, 0), aspect_ratio=:equal, xlims=(-2, 10), ylims=(-4, 8))
scatter!([A.x, B.x, C.x, 4], [A.y, B.y, C.y, 2])
Julia Version 1.8.2 (2022-09-29) が公開されていました
UTF-8 で，ASCII 文字と漢字かな混じりの文字列から，m番目からn番目までの文字列を取り出すとき，Julia では直接的な関数（指定法）がないのかなあ？
index = collect(eachindex(s))
で何文字目か？という情報を取り出し，それに基づいて文字列を抽出する。
s = "123あいうえおabcde漢字仮名ひらがな"
index = collect(eachindex(s))
21-element Vector{Int64}:
1
2
3
4
7
10
13
16
19
20
21
22
23
24
27
30
33
36
39
42
45
s[index[3:5]]
"3あい"
s[index[12:15]]
"de漢字"
関数定義しておく
substring(s, range) = s[collect(eachindex(s))[range]]
そうすれば，
substring(s, 3:5)
"3あい"
substring(s, 12:15)
"de漢字"
0^0 はいくつ？
インターネット上にも多くのページがある。
多くのプログラミング言語では 0^0 = 1 とされている。Julia も同じである。
0^0
1
図を描いてみよう。
using Plots
x = 0:0.001:2
plot(x, x.^x, label="", size=(300, 200))
正の方向から 0 に近づけていくと確かに 1 に近づいていくようだ。
x = 1e-15
x^x
0.9999999999999655
x = big"1" / big"100000000000000000000000000000000000000000000000000000000000000"
x^x
0.999999999999999999999999999999999999999999999999999999999998572397242343691675
SymPy の limit() で見てみよう。
using SymPy
@syms x
limit(x^x, x, 0, dir="+") # dir="+" はデフォルト
1
確かに 1 になる。
ちなみに，x^x が最小値を取るときの x はなんだろう。
f = diff(x^x) # 微分して
x^x(log(x) + 1)
x0 = solve(f)[1] # 傾きが 0 になる x0 を求める。
e^(-1)
x0.evalf()
0.367879441171442
x0^x0.evalf()
0.692200627555346
x が exp(-1) ≒ 0.36788 のとき最小値 0.69220 となる。
Julia で二重根号を外す
二重根号を外すには `sympy.sqrtdenest()` を使う。
using SymPy
x = √(5 + 2√Sym(6))
x |> N
3.1462643699419723423291350657155704455124771291873287012324867174426654953709
y = sympy.sqrtdenest(√(Sym(5) + 2√Sym(6)))
y |> N
3.1462643699419723423291350657155704455124771291873287012324867174426654953709
x - y |> N
0
1. 工夫が必要な例
以下はいずれも単に simplify() するだけでは解けない。
数式の項を取り出す args も役に立つ。
@syms x
y = 1 + 2(x+3) - 3x^2
y.args
(7, -3*x^2, 2*x)
1.1. 例 1
各項を別々に sqrtdenest し，簡約化する。
using SymPy
@syms x, f, s
x = √(2 + √Sym(15)/2) - √(2 - √Sym(15)/2)
f, s = x.args
(sqrt(sqrt(15)/2 + 2), -sqrt(2 - sqrt(15)/2))
f = sympy.sqrtdenest(f)
s = sympy.sqrtdenest(s)
y = f + s |> simplify
x |> N
1.732050807568877293527446341505872366942805253810380628055806979451933016908798
y |> N
1.732050807568877293527446341505872366942805253810380628055806979451933016908798
または，式を2乗して簡約化し，平方根を取る。
x^2 |> simplify |> sqrt
1.2. 例 2
x = 2√Sym(3)/√(Sym(6) -√Sym(27)) - √Sym(2)/√(Sym(9) - 4√Sym(5))
f, s = x.args
(-sqrt(2)/sqrt(9 - 4*sqrt(5)), 2*sqrt(3)/sqrt(6 - 3*sqrt(3)))
f = sympy.sqrtdenest(f) |> simplify
s = sympy.sqrtdenest(s) |> simplify
y = f + s |> simplify
x |> N
-2.127001479758296282603298193936525220323279534045494771601492023532366539644017
y |> N
-2.127001479758296282603298193936525220323279534045494771601492023532366539644017
1.3. 例 3
二重根号を外してから簡約化する。
x = 2/(1 + √Sym(2) + √Sym(3)) + √(2 - √Sym(3))
sympy.sqrtdenest(x) |> simplify
1
x |> N
1.0
