技術ドキュメント: sympy_with_test.py

プログラムの動作

本プログラム sympy_with_test.py は、Pythonのシンボリック計算ライブラリ SymPy を用いて、ユーザーがコマンドラインから指定した数式に対して微分、積分、またはテイラー展開の処理を実行し、その結果を数値計算による近似値と比較して標準出力に表示するスクリプトです。

目的: SymPy によるシンボリックな数式処理の結果が、NumPy を用いた数値計算による近似値とどの程度一致するかを視覚的に検証し、理解を深めることを目的としています。

主な機能:

  • 数式解析: コマンドライン引数として与えられた数式文字列を SymPy が扱える形式に変換します。

  • シンボリック計算:

    • 指定された数式の1階微分を計算します。

    • 指定された数式の不定積分を計算します。

    • 指定された数式を特定の基準点の周りで指定された次数までテイラー展開します。

  • 数値評価: SymPy の結果式および元の数式を、NumPy を用いて特定の \(x\) の範囲で数値的に評価します。

  • 数値比較:

    • 微分モードでは、中央差分近似による数値微分結果と SymPy の微分結果を比較します。

    • テイラー展開モードでは、テイラー展開の結果と元の関数値を比較します。

    • 積分モードでは、現在のところ数値的な比較は行われません。

  • 結果表示: 計算結果、元の関数値、数値近似値、および両者の相対誤差を表形式で標準出力に表示します。

解決する課題: 抽象的なシンボリック計算の結果と、具体的な数値計算の挙動との関連性を、具体的な例を通して理解する手助けとなります。特に、数値微分の近似精度や、テイラー展開の近似範囲における誤差の振る舞いを観察できます。

原理

本プログラムは主に SymPyNumPy という二つの強力なライブラリの機能を利用して、シンボリック計算と数値計算を比較します。

  • 数式解析とシンボリック処理: プログラムはまず、コマンドライン引数で与えられた数式文字列(例: "exp(x)")を sympy.parsing.sympy_parser.parse_expr 関数を使用して SymPy の内部表現である式オブジェクトに変換します。これにより、SymPy が提供するシンボリックな操作が可能になります。

    • 微分: sym.diff(expr, x) を呼び出すことで、式 expr を変数 x について1階微分します。

    • 積分: sym.integrate(expr, x) を呼び出すことで、式 expr を変数 x について不定積分します。

    • テイラー展開: sym.series(expr, x, x0, norder).removeO() を使用して、式 expr を基準点 x0 の周りで norder 次までテイラー展開します。.removeO() は、ランダウの記号 \(O((x-x_0)^n)\) を除去し、テイラー多項式部分のみを抽出します。テイラー展開の原理は以下の級数展開に基づきます。 $\(f(x) = \sum_{n=0}^{\infty} \frac{f^{(n)}(a)}{n!}(x-a)^n = f(a) + f'(a)(x-a) + \frac{f''(a)}{2!}(x-a)^2 + \dots\)$

  • 数値計算と評価: SymPy で得られたシンボリックな式(元の式、微分結果、積分結果、テイラー展開結果など)は、そのままでは具体的な数値として評価できません。そこで、sympy.lambdify(x, sympy_expr, modules=["numpy"]) 関数を使用して、SymPy の式を NumPy の関数として評価できるPython関数に変換します。これにより、指定された \(x\) の値の範囲で、高速に数値計算を実行できます。

  • 数値微分(中央差分近似): 微分モードでは、SymPy によるシンボリックな微分結果と比較するために、数値微分を計算します。これは、点 \(x\) における微分を、ステップサイズ \(h\) を用いて以下の中央差分近似の式で計算します。 $\(\frac{df(x)}{dx} \approx \frac{f(x+h) - f(x-h)}{2h}\)$ この方法は、前進差分や後退差分よりも一般的に高い精度が得られます。

  • 誤差計算: 各 \(x\) の点において、SymPy による計算結果と、数値的な手法(数値微分、元の関数値など)によって得られた値との間の相対誤差が計算されます。相対誤差 \(\delta\) は以下の式で定義されます。 $\(\delta = \frac{|V_{SymPy} - V_{Numeric}|}{|V_{Numeric}|}\)\( ここで \)V_{SymPy}\( は `SymPy` の結果、 \)V_{Numeric}$ は数値計算または元の関数値です。分母が0の場合は0として扱われます。これにより、両者の結果がどの程度一致しているかを定量的に評価できます。

必要な非標準ライブラリとインストール方法

本プログラムを実行するためには、以下のPythonライブラリが必要です。

  • sympy: シンボリックな数式処理を行うためのライブラリ。

  • numpy: 数値計算(配列操作、数学関数など)を行うためのライブラリ。

これらのライブラリは、Pythonのパッケージ管理ツール pip を使用してインストールできます。

pip install sympy numpy

必要な入力ファイル

本プログラムは、コマンドライン引数として数式文字列を直接受け取るため、特別な入力ファイルは必要ありません。

生成される出力ファイル

本プログラムは、現在、ファイルを生成しません。全ての計算結果は標準出力 (stdout) に表形式で表示されます。 なお、--export オプションがコマンドライン引数として定義されていますが、現在の実装ではこのオプションが指定されてもファイルは出力されません。

コマンドラインでの使用例 (Usage)

基本的なコマンドラインの書式は以下の通りです。

python sympy_with_test.py "formula_string" [--mode MODE] [--xmin XMIN] [--xmax XMAX] [--nx NX] [--norder NORDER] [--x0 X0] [-s STEP] [--export]

各引数の説明:

  • formula_string: 処理対象となる数式を文字列として指定します。例えば、"1/(exp(x)+1)", "exp(x)", "x**2 + cos(x)" のように記述します。

  • --mode MODE: 実行する計算モードを指定します。以下のいずれかの値を指定します。

    • diff (デフォルト): 数式の1階微分を計算します。

    • int: 数式の不定積分を計算します。

    • taylor: 数式のテイラー展開を計算します。

  • --xmin XMIN: 数値評価を行う \(x\) の最小値を浮動小数点数で指定します (デフォルト: -1.0)。

  • --xmax XMAX: 数値評価を行う \(x\) の最大値を浮動小数点数で指定します (デフォルト: 1.0)。

  • --nx NX: xmin から xmax の間で等間隔に生成する \(x\) の点の数を整数で指定します (デフォルト: 11)。

  • --norder NORDER: テイラー展開の次数を整数で指定します (デフォルト: 5)。--mode taylor の場合のみ有効です。

  • --x0 X0: テイラー展開の基準点(展開の中心)を浮動小数点数で指定します (デフォルト: 0.0)。--mode taylor の場合のみ有効です。

  • -s STEP, --step STEP: 数値微分における差分計算のステップサイズを浮動小数点数で指定します (デフォルト: 1e-5)。--mode diff の場合のみ有効です。

  • --export: このオプションを指定しても、現在のところファイルは出力されません。

コマンドラインでの具体的な使用例

例1: 関数の微分とその数値微分との比較

数式 \(e^{-x^2}\) の1階微分を計算し、x の範囲を -2 から 2 まで、5 点で評価します。

python sympy_with_test.py "exp(-x**2)" --mode diff --xmin -2 --xmax 2 --nx 5

実行結果の説明: このコマンドは、元の関数 \(f(x) = e^{-x^2}\)\(x\) に関する1階微分をSymPyでシンボリックに計算し、その結果を各 \(x\) 点での数値微分(中央差分近似)と比較します。

====================================================================================================
 Mode: Differentiation (1st)
 Result Formula: -2*x*exp(-x**2)
====================================================================================================
       x |   Original f(x) |   Numeric Diff |   SymPy Result |   Rel.Err
----------------------------------------------------------------------------------------------------
  -2.000 |        0.018316 |       0.073264 |       0.073264 |     1.45e-11
  -1.000 |        0.367879 |       0.735759 |       0.735759 |     6.35e-12
   0.000 |        1.000000 |       0.000000 |       0.000000 |     0.00e+00
   1.000 |        0.367879 |      -0.735759 |      -0.735759 |     6.35e-12
   2.000 |        0.018316 |      -0.073264 |      -0.073264 |     1.45e-11
  • Mode: 実行モードが「Differentiation (1st)」であることが示されます。

  • Result Formula: SymPyによって計算された微分結果の数式 \(-2xe^{-x^2}\) が表示されます。

  • 各行は異なる \(x\) の値に対応し、Original f(x) は元の関数値、\(f(x)\)

  • Numeric Diff は中央差分近似によって計算された数値微分値。

  • SymPy Result はSymPyで得られた微分式の値。

  • Rel.ErrSymPy ResultNumeric Diff の間の相対誤差です。この例では、非常に小さい誤差であり、SymPyのシンボリック微分結果が数値微分と極めてよく一致していることを示しています。

例2: テイラー展開とその元の関数値との比較

数式 \(\cos(x)\)\(x_0=0\) の周りで8次までテイラー展開し、x の範囲を \(-\pi\) から \(\pi\) まで、10 点で評価します。

python sympy_with_test.py "cos(x)" --mode taylor --norder 8 --xmin -3.14159 --xmax 3.14159 --nx 10 --x0 0

実行結果の説明: このコマンドは、元の関数 \(f(x) = \cos(x)\)\(x_0=0\) の周りで8次までテイラー展開し、その展開結果を元の関数 \(\cos(x)\) の値と比較します。

====================================================================================================
 Mode: Taylor Expansion (n=8, x0=0.0)
 Result Formula: x**8/40320 - x**6/720 + x**4/24 - x**2/2 + 1
====================================================================================================
       x |   Original f(x) |   Original f(x) |   SymPy Result |   Rel.Err
----------------------------------------------------------------------------------------------------
  -3.142 |       -1.000000 |       -1.000000 |      -0.976020 |     2.40e-02
  -2.464 |       -0.776662 |       -0.776662 |      -0.778848 |     2.81e-03
  -1.787 |       -0.211475 |       -0.211475 |      -0.211425 |     2.37e-04
  -1.110 |        0.444988 |        0.444988 |       0.444988 |     1.19e-09
  -0.432 |        0.906951 |        0.906951 |       0.906951 |     1.13e-13
   0.245 |        0.969850 |        0.969850 |       0.969850 |     7.51e-15
   0.922 |        0.603378 |        0.603378 |       0.603378 |     3.33e-10
   1.599 |       -0.021028 |       -0.021028 |      -0.021028 |     1.51e-06
   2.276 |       -0.638424 |       -0.638424 |      -0.637292 |     1.77e-03
   2.953 |       -0.977464 |       -0.977464 |      -0.957548 |     2.04e-02
  • Mode: 実行モードと、テイラー展開の次数(n=8)および基準点(x0=0.0)が表示されます。

  • Result Formula: SymPyによって計算された8次テイラー展開の多項式 \(1 - \frac{x^2}{2} + \frac{x^4}{24} - \frac{x^6}{720} + \frac{x^8}{40320}\) が表示されます。

  • Original f(x) は元の関数 \(\cos(x)\) の値。

  • SymPy Result はSymPyで得られたテイラー展開式の値。

  • Rel.ErrSymPy ResultOriginal f(x) の間の相対誤差。基準点 \(x_0=0\) に近いほど誤差は非常に小さく、基準点から離れるにつれて誤差が徐々に大きくなっていることが観察できます。

例3: 不定積分の計算と簡易的な表示

数式 \(x^2 + \cos(x)\) の不定積分を計算します。現在、数値的な比較は実装されていません。

python sympy_with_test.py "x**2 + cos(x)" --mode int

実行結果の説明: このコマンドは、数式 \(f(x) = x^2 + \cos(x)\)\(x\) に関する不定積分をSymPyでシンボリックに計算します。数値積分との比較は行われないため、Numeric Integ* の列は nan (Not a Number) と表示されます。

====================================================================================================
 Mode: Integration (Indefinite)
 Result Formula: x**3/3 + sin(x)
====================================================================================================
       x |   Original f(x) |   Numeric Integ* |   SymPy Result |   Rel.Err
----------------------------------------------------------------------------------------------------
  -1.000 |        1.000000 |              nan |      -0.651377 |     0.00e+00
  -0.800 |        0.896707 |              nan |      -0.457805 |     0.00e+00
  -0.600 |        0.815980 |              nan |      -0.274399 |     0.00e+00
  -0.400 |        0.772592 |              nan |      -0.106547 |     0.00e+00
  -0.200 |        0.767355 |              nan |      -0.007604 |     0.00e+00
   0.000 |        1.000000 |              nan |       0.000000 |     0.00e+00
   0.200 |        1.026756 |              nan |       0.007604 |     0.00e+00
   0.400 |        1.050635 |              nan |       0.106547 |     0.00e+00
   0.600 |        1.088661 |              nan |       0.274399 |     0.00e+00
   0.800 |        1.144901 |              nan |       0.457805 |     0.00e+00
   1.000 |        1.301169 |              nan |       0.651377 |     0.00e+00
  • Mode: 実行モードが「Integration (Indefinite)」であることが示されます。

  • Result Formula: SymPyによって計算された不定積分 \(\frac{x^3}{3} + \sin(x)\) が表示されます。

  • Numeric Integ*: 不定積分の数値的な比較は現在の実装では行われていないため、nan と表示されます。それに伴い相対誤差も0となります。