SymPyの活用法

背景

近頃、ある一定の境界条件を満たす多項式の解析解を必要とすることがありました。
解析解を求めるのは、苦労しそうだったので、SymPyを使って楽にできないか検討してみました。

SymPy

SymPyはシンボリック演算をするためのpythonライブラリです。
以下のようなことができます。

たとえば以下のうような連立方程式の解析解を求めたいとします。
$ bx + ay = ab $
$ x - cy = c$
この連立方程式の解は、

import sympy
a = sympy.symbols("a")
b = sympy.symbols("b")
c = sympy.symbols("c")
x = sympy.symbols("x")
y = sympy.symbols("y")
sympy.solve([b*x+a*y-b*a,x-c*y-c],(x,y))

この結果は、$\frac{a c \left(b + 1\right)}{a + b c},\frac{b \left(a - c\right)}{a + b c}$となって正しい結果が得られます。

多項式の係数

ここでは、多項式

$f(x) = A_{0}+A_{1} x+A_{1} x^{2} + A_{3} x^{3} + A_{4} x^{4} + A_{5} x^{5}$

の中で、

$f(0) = 0$
$f(1) = S_{1}$
$f'(0) = 0$
$f'(1) = 0$
$f''(0) = 0$
$f''(1) = 0$
を満たす係数を求めたいとします。

for i in range(6):
    exec(f"A{i} = sympy.Symbol(\"A{i}\")")
S1 = sympy.Symbol("S1")
x = sympy.Symbol("x")

f = A0 + A1*(x) + A2*(x)**2 + A3*(x)**3 + A4*(x)**4 + A5*(x)**5
f1 = sympy.diff(f, x)
f2 = sympy.diff(f1, x)

まず上記のように変数と求めたい関数を用意して、

$f(0) = 0$の条件は

cond0 = f.subs(x, 0) - 0

のように表現します。
同様に残りすべての条件を

cond1 = f.subs(x, 1) - S1
cond2 = f1.subs(x, 0) - 0
cond3 = f1.subs(x, 1) - 0
cond4 = f2.subs(x, 0) - 0
cond5 = f2.subs(x, 1) - 0

で定義します。

このすべての条件を満たす$A_{0},A_{1},A_{2},A_{3},A_{4},A_{5}$は

solution = sympy.solve([cond0, cond1, cond2, cond3, cond4, cond5],(A0,A1,A2,A3,A4,A5))

とすると、答えは
$(A_{0},A_{1},A_{2},A_{3},A_{4},A_{5}) = (0, 0, 0, 10 S_{1}, - 15 S_{1}, 6 S_{1})$
と求まりました。

この結果を$f(x)$に再代入すると、

for i in range(6):
    exec(f"f = f.subs(A{i}, solution[A{i}])")

$f(x) = 10S_{1}x^{3} -15S_{1}x^{4} + 6S_{1}x^{5}$
$f'(x) = 30S_{1}x^{2} -60S_{1}x^{3} + 30S_{1}x^{4}$
$f''(x) = 60S_{1}x -180S_{1}x^{2} + 120S_{1}x^{3}$

となり、求めたい関数$f(x)$が分かりました。