0

私はPythonに比較的慣れていないので、それを使用して二次非線形微分方程式、特に電解質のポアソン-ボルツマン方程式を解こうとしています。

phi''(r) + (2/r)*phi'(r) = (k^2)*sinh(phi(r))

本質的には、パラメーター k によって支配される減衰率で、電解質内の帯電した表面から離れる静電ポテンシャル (ファイ) の減衰を表します。

  • phi(r) - r でのポテンシャル
  • dphi(r) - r におけるポテンシャルの導関数
  • r - 表面からの距離 (この場合、r = 1 から r = R まで解いています)

と境界条件

  • ファイ(1) = 5
  • dphi(R) = 0

コードの問題ビットは次のとおりです

from scipy.integrate import odeint
from scipy.optimize import root
from pylab import * # for plotting commands

k = 0.5 
phi = 5 
dphi = -10
R = 21

def deriv(z,r): # return derivatives of the array z   (where z = [phi, phi'])
    return np.array([
    (z[1]),
    ((k**2)*sinh(z[0]))-((2/r)*z[1])
    ])


result = odeint(deriv,np.array([phi,dphi]),np.linspace(1,R,1017), full_output = 1)

一般に、k の値が小さい場合、積分は正常に機能し、scipy.optimize の root を使用して境界条件に従って解決できますが、k の値が比較的大きい (0.5 以上) を使用しようとすると、積分に問題が発生し、次のエラーを出力します

Excess work done on this call (perhaps wrong Dfun type).

full_output = 1 で実行し、tcurパラメーターを調べたところ、特定のポイントまでスムーズにカウントアップし、その後、非常に大きな値から非常に小さな値まで激しく振動しているようです。同じ時点nfeで、関数評価の数はゼロに減少します。正しく動作している場合、tcur パラメーターは 1 から R までスムーズに実行されます。私の実装に問題がありますか、それとも方程式に不安定性がありますか?

助けてくれてありがとう、

デイブ

4

1 に答える 1

0

ODE はおそらく不安定です。関連する式

phi''(r) = (k^2)*sinh(phi(r))

k=1(対応する r=1 の初期条件)の解があります。

phi(r) = 2 arctanh(sin(r))

解は で特異点を持っていr=pi/2ます。数値 ODE ソルバーは、この点を超えて方程式を積分することはできません。一次微分項 (とにかく特異点の近くでは無視できるはずです) を持つ同様の方程式が同様に動作することはもっともらしいです。

実際の問題は、ODE ソルバーを使用した射撃法は、境界値問題を解決するのに適していないということです。かなり安定している選点法を使用する必要があります。たとえば、 http ://www.scholarpedia.org/article/Boundary_value_problemおよびその中の参照を参照してください。

Python ソフトウェアについては、https://pypi.python.org/pypi?%3Aaction=search&term=boundary+value+problem&submit=searchを参照してください。

通常、このようなソルバーを自分で作成することも非常に簡単です。必要なステップは、問題を (多数の) 方程式のセットに離散化することだけであり、その後でrootそれらを解くことができるからです。

于 2014-12-14T13:58:13.407 に答える