代数計算(微分・積分・行列計算)

このページではGitHubにあげたcalculate_algebra.ipynbの中身をマークダウン形式で書き下します。

Info

以下ではipynbでのoutputをこの中に表示します。

Tip

このようにした理由としてはGitHubのviwerの共同が安定していないためです。これまではみれていたソースが急に見れなくなったりしているようです。

代数計算(微分・積分・行列計算)を行ってみる

代数計算ライブラリsympyを用いて簡単な代数計算を行い、できることの一例を示します。

まずは必要なライブラリをimportしてきます。

from fumi_python_libraries import BasicAlgebraicCalculationClass
import numpy
import sympy

インスタンスを生成します。

alg_calc = BasicAlgebraicCalculationClass()

基本操作

出力形式の確認

まずはデフォルトの値で出力を見てみます。

BasicAlgebraicCalculationClassのクラスではprint_formula()という関数が実装されています。

select_print_mode( bool )で表示形式を判別する変数latex_paste_modeを設定できます。

デフォルトではこの変数はFalseで、IPython.diplay.display()で描画されます（Jupyter Notebookで数式を綺麗に出力してくれる関数です。）

一方でTrueではprint()で出力します。この形式ではLatex(Overleafやcloud latex)、Keynoteの数式モードにペーストして使える形でダンプしてくれます。

alg_calc.print_formula()

Info

$$ f(x) = x $$

latex_paste_mode = Trueにして見てみましょう。

alg_calc.select_print_mode(True)
alg_calc.print_formula()

Info

f(x) = x \\

よく使用するメソッドの使い方

代数の宣言と代数計算を行う。

まずは宣言時ですが、引数として生成する変数リストを渡すことができます。

この際、生成した変数はクラス内の変数リストに格納されるのでget_variable_symbols()を用いて取得することができます。

alg_calc = BasicAlgebraicCalculationClass(["\\alpha","\\beta","x"])

var = alg_calc.get_variable_symbols()

上記のようにして渡された変数は普通のプログラム同様に四則演算ができます。

fun = ( var[0] ) / ( var[2]**3 + var[1] )

このようにして計算した関数をBasicAlgebraicCalculationClassの右辺に渡しましょう。

メソッドはset_right_hand_side(function)です。

alg_calc.set_right_hand_side(fun)
alg_calc.print_formula()

Info

$$ f(x) = \frac{\alpha}{\beta + x^{3}} $$

同様にして、set_left_hand_side(function)左辺にも値を渡すことができます。

ただし、set_left_hand_side(function)の時渡す引数の方は文字列を前提にしています。

これは後述する行列計算の時におかしな描かれ方をすることを防ぐためです。

sympyで宣言したもの四則演算の結果を描画する際は代数を名前の順などにソートして出力します。

行列の代表させた大文字の $A$ などの値をいれて $A^{-1}PA$ などを描くと出力した時に $A^{-1}AP$ などを表示してきます。

そのため少し面倒ですが、左辺は文字列で書くことにしています。

では同様にして左辺を代入してみましょう。

alg_calc.set_left_hand_side('f(\\alpha,\\beta,x)')
alg_calc.print_formula()

Info

$$ f(\alpha,\beta,x) = \frac{\alpha}{\beta + x^{3}} $$

新たに変数を生成する場合はcreate_variables()を使用します。

var = alg_calc.create_variables(["a","b","x"])
alg_calc.set_left_hand_side('f(a,b,x)')
fun = ( var[0] ) / ( var[2]**3 + var[1] )
alg_calc.set_right_hand_side(fun)
alg_calc.print_formula()

Info

$$ f(a,b,x) = \frac{a}{b + x^{3}} $$

またこのクラスが持つ左辺と右辺に変数に値を渡さずにダンプするメソッドprint_external_input(left, right)もあります。

例えば以下のようできます。

alg_calc.print_external_input("x(t)", var[0] * sympy.sin ( var[2] ) + var[1] * sympy.sin ( var[2] ) + abs( var[0] + var[2]**var[1] ) )

Info

$$ x(t) = a \sin{\left(x \right)} + b \sin{\left(x \right)} + \left|{a + x^{b}}\right| $$

念の為もう一度print_formula()を用いてクラス内の左辺と右辺の式を格納している変数が更新されていないことを確認します。

alg_calc.print_formula()

Info

$$ f(a,b,x) = \frac{a}{b + x^{3}} $$

数値を代入する

次に値を代入してみましょう。値を代入して出力する時は複数のメソッドを使用します。まず、set_variable_symbols(　list )にはクラスが持つ変数のうち値を代入したり、微分積分をしたりする変数を指定します。

ここではa, bを定数として、これらの変数に値を代入してみましょう。

alg_calc.set_variable_symbols(["a","b"])

上のメソッド自体は戻り値を持たないので、そのまま代入の操作を行なっていきます。

代入はassign_values(variable list,valule list)を使用します。

今回は $a=1$ , $b=2$ を代入してみます。

alg_calc.assign_values(["a","b"],[1,2])

上記の操作を行うと生成したインスタンス内の右辺ではない変数に一時的値が格納されます。

この一時的に格納された値はget_variable_values()を使用して取得できます。

以下のセルでは取得したのち、print_external_input()でダンプてみます。

alg_calc.print_external_input('\\mathrm{Asigned\\,function}', alg_calc.get_variable_values())

Info

$$ \mathrm{Asigned\,function} = \frac{1}{x^{3} + 2} $$

右辺を更新したい場合は、set_right_hand_side()に値を渡すか、assign_values()の第4引数をTrueにすると更新できます。（第3引数については丸め込みのフラグです。この後説明します。）

このように代入後に変数を更新しない理由はloop処理ないで毎回別の値を代入する時にいちいちloopの中でインスタンスを生成したり、関数を定義しなおしたりしなくても済むようにしているためです。

var = alg_calc.create_variables(["a","b","x"])
alg_calc.set_left_hand_side('f(a,b,x)')
fun = ( var[0] ) / ( var[2]**3 + var[1] )
alg_calc.set_right_hand_side(fun)
alg_calc.print_formula()

Info

$$ f(a,b,x) = \frac{a}{b + x^{3}} $$

alg_calc.set_variable_symbols(["a","b"])
alg_calc.assign_values(["a","b"],[1,2],None,True)
alg_calc.print_formula()

Info

$$ f(a,b,x) = \frac{1}{x^{3} + 2} $$

計算結果/出力結果を丸める

またsympyで出力すると、表示される桁がじどうで調整されないので割り切れない数字を入れると大変なことになったりします。

そのためいくつかの丸め込み用の引数やメソッドが用意されています。

例えば循環小数 $\dfrac{1}{19}=0.\overline{052631578947368421}$ を表示してみましょう。

function = var[0] / ( var[1] + var[2] )
values = [1, 19, 0]
labels = ["a", "b", "x"]
alg_calc.set_right_hand_side( function )
alg_calc.set_variable_symbols( labels )
alg_calc.assign_values( labels, values )
alg_calc.set_left_hand_side( "\\frac{1}{19}" )
alg_calc.set_right_hand_side( alg_calc.get_variable_values() )
alg_calc.print_formula()

Info

$$ \frac{1}{19} = 0.0526315789473684 $$

こんな感じで長ったらしい数字が羅列されても困ってしまうことがあります。

これをなんとかする方法の一つとして、先ほど使用した代入するメソッドassign_values()の第3引数にて丸める桁指定できます。

書き方はassign_values(variable list,valule list, round digit)です。

このメソッドでは代入する前に丸めることになります。

例えば以下のような感じになります。

values = [1, 10.5000000001, 9.4999999999]
labels = ["a", "b", "x"]
alg_calc.set_right_hand_side( function )
alg_calc.set_variable_symbols( labels )
alg_calc.assign_values( labels, values , 6)
alg_calc.set_left_hand_side( "\\mathrm{Round}\\left(\\frac{1}{19},0\\right)" )
alg_calc.set_right_hand_side( alg_calc.get_variable_values() )
alg_calc.print_formula()

Info

$$ \mathrm{Round}\left(\frac{1}{19},0\right) = 0.05 $$

ただしこの場合は代入される前に丸められるので上記のように6桁目で丸めるようにしてもアウトプットは小数点2桁になってしまいます。

それもそのはずで、代入前に四捨五入しているので、 $\frac{1}{10.5+9.5}=\frac{1}{20}$ になっているので想定以上に丸め込まれていることがわかると思います。

このような場合には代入する前ではなく、後に丸めるたくなります。

そのためのメソッドがround_right_hand_side( round_digits )です。

function = var[0] / ( var[1] + var[2] )
values = [1, 19, 0]
labels = ["a", "b", "x"]
alg_calc.set_right_hand_side( function )
alg_calc.set_variable_symbols( labels )
alg_calc.assign_values( labels, values )
alg_calc.set_left_hand_side( "\\mathrm{Round}\\left(\\frac{1}{19},6\\right)" )
alg_calc.set_right_hand_side( alg_calc.get_variable_values() )
alg_calc.round_right_hand_side(6)
alg_calc.print_formula()

Info

$$ \mathrm{Round}\left(\frac{1}{19},6\right) = 0.052632 $$

上記のround_right_hand_side( round_digits )では元々の右辺を更新するので計算を終えて、最後に表示する時のみ使用することが推奨されます。

ObjectをCopyする

さらにクラスをコピーするメソッドもあります。このメソッドはモンテカルロなどにおいてloop処理ないで逐次代入する際にいちいち新しいクラスを宣言する必要がないようにするためです。

copied_alg_calc = alg_calc.copy_class()
copied_alg_calc.print_formula()

Info

$$ \mathrm{Round}\left(\frac{1}{19},6\right) = 0.052632 $$

微分積分をおこなう

sympyが持つ微分積分をおこなうメソッドを用いた微分・積を行うメソッドが用意されています。

その方法を書いておきます。

不定積分をしてみる

不定積分はindefinite_integral(variable, update flag)を使います。

まずは積分する前に新たに積分するための関数をわたします。

一度初期化しても良いですがset_formula（lhs, rhs)で外部から一度に右辺と左辺を登録できます。

alg_calc.select_print_mode(False)
var = alg_calc.create_variables( ["a", "x", "b", "y", "c","\\alpha", "\\beta"] )
fun = var[0] * var[1] + var[2] * var[3]**var[5] + var[4]**var[6]
alg_calc.set_formula( 'f(x,y)', fun )
alg_calc.print_formula()

Info

$$ f(x,y) = a x + b y^{\alpha} + c^{\beta} $$

あとはindefinite_integral(variable, update_flag)で不定積分ができます。

第2引数は積分結果で右辺を更新するかを選択できます。

ind_intg = alg_calc.indefinite_integral( var[3], True )
alg_calc.print_formula()

Info

$$ \int f(x,y)dy = a x y + b \left(\begin{cases} \frac{y^{\alpha + 1}}{\alpha + 1} & \text{for}\: \alpha \neq -1 \\\log{\left(y \right)} & \text{otherwise} \end{cases}\right) + c^{\beta} y $$

定積分をしてみる

次に定積分です。こちらはdefinite_integral( variable , min, max, update flag )でできます。

int_range = [0, 2]
alg_calc.set_formula( 'f(x,y)', fun )
def_intg = alg_calc.definite_integral( var[1] , *int_range, True )
alg_calc.print_formula()

Info

$$ \int_{0}^{2}f(x,y)dx = 2 a + 2 b y^{\alpha} + 2 c^{\beta} $$

微分してみる

次に微分です。こちらはdifferentiate( variable, update flag )です。

alg_calc.set_formula( 'f(x,y)', fun )
dif_func = alg_calc.differentiate( var[4], True )
alg_calc.print_formula()

Info

$$ \frac{d}{dc}f(x,y) = \frac{\beta c^{\beta}}{c} $$

行列計算をしてみる。

最後に行列計算です。

行列の計算では生成するためのメソッドなどが少し異なります。

create_martix_symbols( matrixsymbol, row, col, elementsymbol)で生成します。

第四引数は行列の各要素を代表する記号を指定したい時に使用します。

デフォルトではNoneになっており、この場合は第一引数で指定した文字の小文字が採用されます。

戻り値は行列そのものと各要素の代数を格納したリストです。

n = 10
A, a = alg_calc.create_martix_symbols( "A", n, n )
alg_calc.set_formula( 'A' , A )
alg_calc.print_formula()

Info

$$ A = \left[\begin{matrix}a_{1 1} & a_{1 2} & a_{1 3} & a_{1 4} & a_{1 5} & a_{1 6} & a_{1 7} & a_{1 8} & a_{1 9} & a_{1 10}\\a_{2 1} & a_{2 2} & a_{2 3} & a_{2 4} & a_{2 5} & a_{2 6} & a_{2 7} & a_{2 8} & a_{2 9} & a_{2 10}\\a_{3 1} & a_{3 2} & a_{3 3} & a_{3 4} & a_{3 5} & a_{3 6} & a_{3 7} & a_{3 8} & a_{3 9} & a_{3 10}\\a_{4 1} & a_{4 2} & a_{4 3} & a_{4 4} & a_{4 5} & a_{4 6} & a_{4 7} & a_{4 8} & a_{4 9} & a_{4 10}\\a_{5 1} & a_{5 2} & a_{5 3} & a_{5 4} & a_{5 5} & a_{5 6} & a_{5 7} & a_{5 8} & a_{5 9} & a_{5 10}\\a_{6 1} & a_{6 2} & a_{6 3} & a_{6 4} & a_{6 5} & a_{6 6} & a_{6 7} & a_{6 8} & a_{6 9} & a_{6 10}\\a_{7 1} & a_{7 2} & a_{7 3} & a_{7 4} & a_{7 5} & a_{7 6} & a_{7 7} & a_{7 8} & a_{7 9} & a_{7 10}\\a_{8 1} & a_{8 2} & a_{8 3} & a_{8 4} & a_{8 5} & a_{8 6} & a_{8 7} & a_{8 8} & a_{8 9} & a_{8 10}\\a_{9 1} & a_{9 2} & a_{9 3} & a_{9 4} & a_{9 5} & a_{9 6} & a_{9 7} & a_{9 8} & a_{9 9} & a_{9 10}\\a_{10 1} & a_{10 2} & a_{10 3} & a_{10 4} & a_{10 5} & a_{10 6} & a_{10 7} & a_{10 8} & a_{10 9} & a_{10 10}\end{matrix}\right] $$

第2の戻り値は行列の各要素に値を代入する時に使用されます。

以下では例えば、上記の $10 \times 10$ の行列に対して後からヒルベルト行列を代入してみます。

まずはnumpyを使ってヒルベルト行列を数値で計算します。

def hilbert_matrix( n ):
  H = numpy.zeros( ( n, n ) )
  for i in range( n ):
    for j in range( n ):
      H[ i, j ] = 1 / ( i + j + 1 )
  return H
    
h_np = hilbert_matrix( n ).flatten()

上記の式では最後に行列を平滑化しています。

これはassign_values()で代入する際に渡すのものがリストのためです。

今回は代入した後にsympy.nsimplify( function, rational=True ) )を用いて分数にしてから出力しています。

これはみやすさ重視のためです。（ただし、結構適当に値の近い分数を取ってくることもあるので注意）

alg_calc.set_formula( 'A' , A )
alg_calc.assign_values( a, h_np, None, True )
h_sy = alg_calc.get_right_hand_side()
alg_calc.set_formula( 'H', sympy.nsimplify( h_sy, rational=True ) )
alg_calc.print_formula()

Info

$$ H = \left[\begin{matrix}1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10}\\\frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11}\\\frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12}\\\frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13}\\\frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14}\\\frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15}\\\frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16}\\\frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17}\\\frac{1}{9} & \frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18}\\\frac{1}{10} & \frac{1}{11} & \frac{1}{12} & \frac{1}{13} & \frac{1}{14} & \frac{1}{15} & \frac{1}{16} & \frac{1}{17} & \frac{1}{18} & \frac{1}{19}\end{matrix}\right] $$

上記の処理において肝心の計算部分では数値のみを使っていたので、当然成功していてほしいわけです。

しかしsympyのいいところは代数のまま色々計算できるところです。

そのため最後に、完全に代数のみで対角化をしたのちに値を代入してみます。

今回は $2\times 2$ の行列を用意します。

n = 2
L, l = alg_calc.create_martix_symbols( "Lambda", n, n )
alg_calc.set_formula( "\\Lambda" , L )
alg_calc.print_formula()

Info

$$ \Lambda = \left[\begin{matrix}\lambda_{1 1} & \lambda_{1 2}\\\lambda_{2 1} & \lambda_{2 2}\end{matrix}\right] $$

sympyのdiagonalize()メソッドを使って対角化します。

戻り値は対角行列を $\Lambda=PDP^{-1}$ とする時、 $P$ , $D$の順です。

p, d = alg_calc.right_hand_side.diagonalize()
alg_calc.print_external_input("P",p)
alg_calc.print_external_input("D",d)

Info

$$ P = \left[\begin{matrix}\frac{\lambda_{1 1} - \lambda_{2 2} - \sqrt{\lambda_{1 1}^{2} - 2 \lambda_{1 1} \lambda_{2 2} + 4 \lambda_{1 2} \lambda_{2 1} + \lambda_{2 2}^{2}}}{2 \lambda_{2 1}} & \frac{\lambda_{1 1} - \lambda_{2 2} + \sqrt{\lambda_{1 1}^{2} - 2 \lambda_{1 1} \lambda_{2 2} + 4 \lambda_{1 2} \lambda_{2 1} + \lambda_{2 2}^{2}}}{2 \lambda_{2 1}}\\1 & 1\end{matrix}\right] $$ $$ D = \left[\begin{matrix}\frac{\lambda_{1 1}}{2} + \frac{\lambda_{2 2}}{2} - \frac{\sqrt{\lambda_{1 1}^{2} - 2 \lambda_{1 1} \lambda_{2 2} + 4 \lambda_{1 2} \lambda_{2 1} + \lambda_{2 2}^{2}}}{2} & 0\\0 & \frac{\lambda_{1 1}}{2} + \frac{\lambda_{2 2}}{2} + \frac{\sqrt{\lambda_{1 1}^{2} - 2 \lambda_{1 1} \lambda_{2 2} + 4 \lambda_{1 2} \lambda_{2 1} + \lambda_{2 2}^{2}}}{2}\end{matrix}\right] $$

上記のように対角化を計算できました。

もちろん得られた行列についても代数計算可能です。

まず、対角行列を計算してみましょう。

alg_calc.set_right_hand_side( p**-1 * L * p )
alg_calc.set_left_hand_side( 'P\\Lambda P^{-1}' )
alg_calc.print_formula()

Info

$$ P\Lambda P^{-1} = \Biggl[ \mathrm{omit} \Biggr] $$

このアウトプットは非常に長かったので、このページでは省略します。

結構形が複雑すぎるのでここでは省略

自分で式を書いていたらかなり苦しそうです…

この対角行列に値を代入していきましょう。

l_np = hilbert_matrix( n ).flatten()
alg_calc.set_right_hand_side( p**-1 * L * p )
alg_calc.assign_values( l, l_np)
alg_calc.set_formula( 'PHP^{-1}', alg_calc.get_matrix( alg_calc.get_variable_values() ) )
alg_calc.round_right_hand_side(10)
alg_calc.print_formula()
php = alg_calc.get_right_hand_side()

Info

$$ PHP^{-1} = \left[\begin{matrix}0.0657414541 & 0\\0 & 1.2675918792\end{matrix}\right] $$

描画を綺麗にするために10桁で丸めました。

さてこの値が正しいかを確認する方法ですが、 $P$ の行列（逆行列）をもう一度かけて値を代入した時にヒルベルト行列になっていれば良いはずです。

まずは行列 $P$ を上記の対角行列にかけます。

この際、下のセル内のpは代数なので出力結果には $\lambda_{ij}$ が残ります。

pphpp = p * php * p**-1
alg_calc.print_external_input( 'P^{-1}\\left(PH P^{-1}\\right)P', pphpp )

Info

$$ P^{-1}\left(PH P^{-1}\right)P = \Biggl[ \mathrm{omit} \Biggr] $$

このアウトプットは非常に長かったので、このページでは省略します。

よってもう一度、代入して分数の形で出力してみます。

alg_calc.set_right_hand_side( pphpp )
alg_calc.assign_values( l, l_np, None, True )
h_sy = alg_calc.get_right_hand_side()
alg_calc.set_formula( "H", sympy.nsimplify( h_sy, rational=True ) )
alg_calc.print_formula()

Info

$$ H = \left[\begin{matrix}1 & \frac{1}{2}\\\frac{1}{2} & \frac{1}{3}\end{matrix}\right] $$

非常に良い結果を得ました。