2017-12-24 10 views
0

は、私が遭遇した問題を提起してみましょう:考える
は、常微分方程式DX/DT = F(xは、t)は位相肖像次のコードで初期値の配列のために位相軌跡をプロットすることによりプロットされる:(DX/DT = xで見て 複数のデータセットを計算し、プロットするときmatplotlibの/ scipyのダウンロードは、(のODEの位相肖像画を)奇妙な結果を生成

import matplotlib.pyplot as plt 
import numpy as np 
from scipy.integrate import odeint 

def phase_trajectories(func, inits, xbound, ybound): 
    """func: RHS of ODE to be integrated by odeint 
    inits: list or array of initial-value-tuples 
    xbound, ybound: Boundaries for plotting 
    """ 
    t = np.linspace(0, 10, 100) # solution is calculated for these values 
    axis = plt.gca() 
    def f_neg(x, t): 
     return -func(x, t) 
    for i in inits: 
     sol_fwd = odeint(func, i, t) # calc solution forward 
     sol_bwd = odeint(f_neg, i, t) # and backward in time 
     sol = np.vstack((np.flipud(sol_bwd), sol_fwd)) # put both together 
     sol_x = sol[:,0] 
     sol_y = sol[:,1] 
     sol_x_masked = np.ma.masked_outside(sol_x, *xbound) # masking the data because 
     sol_y_masked = np.ma.masked_outside(sol_y, *ybound) # of blow-up of the solution 
     axis.plot(sol_x_masked, sol_y_masked, c='b') 
    plt.plot() 

具体例としてt)²初期値はx(t0)= x0

# inits contains tuples (t0, x0) 
inits = np.array([(i, j) for i in range(-2, 3) for j in (-1, 1)]) 

def func(x, t): 
    # receives x as array with shape (n,) from odeint 
    # returns [dt/dt, dx/dt] 
    return np.array([1, x[1]**2]) 

phase_trajectories(func, inits, (-2, 2), (-2, 2))プロット予想される軌道を呼び出すだけでなく、関数が呼び出されるたびに異なるように思われる不要なノイズ追加:、この底面を得るために example with noise
another example

を私はphase_trajectoriesを一度に全体のinits配列ではなく、最初の値のペアごとに別々に、nextメソッドを手動で呼び出すジェネレータを使用して消耗するまで:

def generator_pt(): 
    for i in inits: 
     yield phase_trajectories(func, (i,), (-2, 2), (-2, 2)) 

g = generator_pt() 
next(g) 

これは、しかし、考え、時には(妙にないたびに)私に望ましい結果を与える:
It's supposed to look like that

を私は与えられたODEのソリューションは、それが唯一の特定のintervall上に存在する意味、吹くことを知っていますその限界に近づくと分岐するので、私はデータを保持している配列をマスクしてプロットしています。
それにもかかわらず、データが計算され、プロットされるたびに、不要なノイズが異なる理由は、私には謎です。

ここに述べられた所見に基づいて、私はこの不正行為の原因がmatplotlibの中で求められているか、むしろ私の側でこのライブラリの恐ろしい誤用を求めていると思う。 私はすでに異なるmatplotlibバックエンドを使い、必要なpackegesが更新された別の環境(Python 2.7と3.5)でコードを実行しようとしましたが、これまでのところ苦労していませんでした。 。

おそらく、誰かが私にヒントを与えたり、これらの結果を説明したり、少なくともこのようなことが起こっている理由を理解するのに役立つ洞察を提供することができます。

StackOverflowの、あなたは私の唯一の希望...

答えて

0

ある問題は数値的なODEを解決するために、古い実装ですscipys odeint Methodeの、内にあるようです。 scipy v1.0.0以来、妨害ノイズが現れないようにsolve_ivpの新しい実装が利用可能です。参考

import matplotlib.pyplot as plt 
import numpy as np 
from scipy.integrate import solve_ivp 

def phase_trajectories(func, inits, xbound, ybound): 
    axis = plt.gca() 
    t = (0, 10) 
    def f_neg(t, x): 
     return -func(t, x) 
    for f in (func, f_neg): 
     for i in inits: 
      x, y = solve_ivp(f, t, i).y 
      x_ma = np.ma.masked_outside(x, *xbound) 
      y_ma = np.ma.masked_outside(y, *ybound) 
      axis.plot(x_ma, y_ma, c='b') 
    plt.plot() 

関数と初期値looks as aspectedに対する結果。しかし、なぜodeintがその奇妙な振る舞いをしているのかは未だに未だにわかっておらず、根底にあるFORTRANの実装に見出されるべきである。

関連する問題