📂ジュリア

ジュリアのシンボリック演算パッケージSymbolics.jlの紹介

概要

Juliaのシンボリック代数システムであるSymbolics.jlについて紹介します¹。このパッケージは、特にJuliaの基本文法と一緒に、非常に直感的で強力なインターフェースを提供します。

SymEngine.jlとの違い

Symbolics.jlはJuliaでネイティブに実装されており、パフォーマンスやインターフェースの面で多くの優れた点があります。Juliaでシンボリック演算をする方法を紹介したポストで紹介されたSymEngine.jlは、元々C++で書かれたライブラリをJuliaでラッピングしたものです。ラッピングされているため、当時は特に機能開発なしで使いやすかったかもしれませんが、今では元々SymEngineを使い慣れていない限り、わざわざSymEngine.jlを使う理由はありません。

コード

変数宣言 @variables

julia> using Symbolics

julia> @variables t x y z
4-element Vector{Num}:
 t
 x
 y
 z

julia> w = log(y) + cos(x) + 2t^3
log(y) + cos(x) + 2(t^3)

シンボリックに使用する変数は、@variablesの後に列挙する形で宣言できます。これらの変数は自然に超越関数の引数などにも使用できます。

展開 expand

julia> (x + y)^3
(x + y)^3

julia> expand((x + y)^3)
x^3 + 3(x^2)*y + 3x*(y^2) + y^3

expand関数を通じて展開ができます。関数を取らない場合、パフォーマンス上の理由などで、特に指示されていない展開は行われません。

簡略化 simplify

julia> simplify(2x + (cos(x))^2 + 3x + (sin(x))^2)
1 + 5x

simplify関数を通じて、0になる項を取り除いたり、係数同士を合わせて簡単にすることができます。例で見るように、三角関数の二乗の和が1になることを自動で計算してくれます。

置換と代入 substitute

substitute関数を通じて、変数を置換したり代入することができます。どのように置換が行われるかは、辞書を作って入れてあげればいいです。 $$ w = \log y + \cos x + 2t^3 $$ この$w$に置換と代入をしてみましょう。

julia> substitute(w, Dict(t => x))
log(y) + cos(x) + 2(x^3)

シンボリック変数同士の置換はもちろん可能です。

julia> substitute(w, Dict(t => 2, x => π, y => exp(1)))
16.0

$t = 2$、$x = \pi$、$y = e$に代入した結果は数値的に計算してリターンされます。

微分 Differential

julia> Dt = Differential(t)
(::Differential) (generic function with 3 methods)

julia> Dt(w)
Differential(t)(log(y) + cos(x) + 2(t^3))

julia> expand_derivatives(Dt(w))
6(t^2)

特にSymbolicsでは、直接微分する関数に対して微分する変数を引数として渡すのではなく、特定の変数に対する微分を行う微分作用素を定義して、ここにシンボルを渡す方式で微分を行います。一見するとコードが汚く複雑に見えるかもしれませんが、数学者の視点から見ると非常に直感的で妥当なインターフェースだと言えます。微分した結果は、展開と同じように直ちに計算されず、expand_derivatives関数を通じて結果がリターンされます。

julia> Dx = Differential(x)
(::Differential) (generic function with 3 methods)

julia> Dy = Differential(y)
(::Differential) (generic function with 3 methods)

julia> Dx(w) |> expand_derivatives
-sin(x)

julia> Dy(w) |> expand_derivatives
1 / y

作用素を定義するという観点から見れば、偏微分も非常に自然に実装されます。

ベクトルと行列

julia> @variables σ ρ β
3-element Vector{Num}:
 σ
 ρ
 β

julia> Lorenz = [-σ*x + σ*y
                 -x*y + ρ*x - y
                  x*y - β*z]
3-element Vector{Num}:
     -x*σ + y*σ
 -y - x*y + x*ρ
      x*y - z*β

Juliaの文法そのままにベクトルを作ればいいです。例としてローレンツアトラクタの右辺も上記のように簡単にシンボルで表現できます。

julia> M = [x*y y^2
             2x   9]
2×2 Matrix{Num}:
 x*y  y^2
  2x    9

行列も非常に直感的に定義されます。

ヤコビアン

julia> Symbolics.jacobian(Lorenz, [x, y, z])
3×3 Matrix{Num}:
     -σ       σ   0
 -y + ρ  -1 - x   0
      y       x  -β

ヤコビ行列はSymbolics.jacobianを通じて簡単に求めることができます。

全体のコード

using Symbolics

@variables t x y z
w = log(y) + cos(x) + 2t^3

(x + y)^3
expand((x + y)^3)

simplify(2x + (cos(x))^2 + 3x + (sin(x))^2)

substitute(w, Dict(t => x))
substitute(w, Dict(t => 2, x => π, y => exp(1)))

Dt = Differential(t)
Dt(w)
expand_derivatives(Dt(w))

Dx = Differential(x)
Dy = Differential(y)
Dx(w) |> expand_derivatives
Dy(w) |> expand_derivatives

@variables σ ρ β
Lorenz = [-σ*x + σ*y
          -x*y + ρ*x - y
           x*y - β*z]

M = [x*y y^2
      2x   9]

Symbolics.jacobian(Lorenz, [x, y, z])

環境

• OS: Windows
• julia: v1.10.0
• Symbolics v5.16.1