11

Pythonを使用して数値的に積分を解いています:

ここに画像の説明を入力

ここで、a(x)は任意の値を取ることができます。正、負、[-1;1] の内側または外側、および eta は極小の正の量です。a(x) の値を変更する 2 番目の外積分があります。

Sokhotski-Plemelj の定理を使用してこれを解決しようとしています。 ここに画像の説明を入力

ただし、これには原則の値を決定することが含まれますが、これは Python でメソッドを見つけることができません。私はそれがMatlabで実装されていることを知っていますが、Pythonでプリンシパル値を決定するライブラリまたは他の方法を知っている人はいますか(プリンシパル値が存在する場合)?

4

1 に答える 1

7

sympy を使用して積分を直接評価できます。eta->0 の実部が主要な値です。

from sympy import *
x, y, eta = symbols('x y eta', real=True)
re(integrate(1/(x - y + I*eta), (x, -1, 1))).simplify().subs({eta: 0})
# -> log(Abs(-y + 1)/Abs(y + 1))

もちろん、 Matlab のシンボリック ツールボックスintでも同じ結果が得られます (これに関連する Matlab の他のツールについては知りません --- 特定のツールを知っている場合は指定してください)。

主値の数値計算について質問されました。f(y)その答えは、解析形式や動作がわからない関数しかない場合、それらを数値的に計算することは一般に不可能だということです。被積分関数の極がどこにあり、どのような順序であるかなどを知る必要があります。

一方、積分が の形式f(y) / (y - y_0)であることがわかっている場合scipy.integrate.quadは、主値を計算できます。たとえば、次のようになります。

import numpy as np
from scipy import integrate, special

# P \int_{-1}^1 dx 1/(x - wvar) * (1 + sin(x))
print(integrate.quad(lambda x: 1 + np.sin(x), -1, 1, weight='cauchy', wvar=0))
# -> (1.8921661407343657, 2.426947531830592e-13)

# Check against known result
print(2*special.sici(1)[0])
# -> 1.89216614073

詳しくはこちらをご覧ください。

于 2013-10-09T18:08:57.993 に答える