錯誤在Python

Question

我求解非線性Schrodinger（NLS）方程RK4算法：

(1): i*u_t + 0.5*u_xx + abs(u)^2 * u = 0 

施加傅立葉變換後，就變成：

(2): uhat_t = -0.5*i*k^2 * uhat + i * fft(abs(u)^2 * u) 

其中uhat是傅立葉變換的u 。上面的公式（2）是一個相當規範的IVP，可以用第四或者Runge-Kutta方法求解。這裏是我的解決方程（2）代碼：

import numpy as np 
import math 
from matplotlib import pyplot as plt 
from matplotlib import animation 


#----- Numerical integration of ODE via fixed-step classical Runge-Kutta ----- 

def RK4(TimeSpan,uhat0,nt): 
    h = float(TimeSpan[1]-TimeSpan[0])/nt 
    print h 
    t = np.empty(nt+1) 
    print np.size(t)  # nt+1 vector 
    w = np.empty(t.shape+uhat0.shape,dtype=uhat0.dtype) 
    print np.shape(w)  # nt+1 by nx matrix 
    t[0] = TimeSpan[0]  
    w[0,:] = uhat0   # enter initial conditions in w 
    for i in range(nt): 
     t[i+1] = t[i]+h 
     w[i+1,:] = RK4Step(t[i], w[i,:],h) 
    return w 

def RK4Step(t,w,h): 
    k1 = h * uhatprime(t,w) 
    k2 = h * uhatprime(t+0.5*h, w+0.5*k1*h) 
    k3 = h * uhatprime(t+0.5*h, w+0.5*k2*h) 
    k4 = h * uhatprime(t+h,  w+k3*h) 
    return w + (k1+2*k2+2*k3+k4)/6. 

#----- Constructing the grid and kernel functions ----- 
L = 40 
nx = 512 
x = np.linspace(-L/2,L/2, nx+1) 
x = x[:nx] 

kx1 = np.linspace(0,nx/2-1,nx/2) 
kx2 = np.linspace(1,nx/2, nx/2) 
kx2 = -1*kx2[::-1] 
kx = (2.* np.pi/L)*np.concatenate((kx1,kx2)) 

#----- Define RHS ----- 
def uhatprime(t, uhat): 
    u = np.fft.ifft(uhat) 
    z = -(1j/2.) * (kx**2) * uhat + 1j * np.fft.fft((abs(u)**2) * u) 
    return z 

#------ Initial Conditions ----- 
u0  = 1./np.cosh(x)#+1./np.cosh(x-0.4*L) 
uhat0 = np.fft.fft(u0) 

#------ Solving for ODE ----- 
TimeSpan = [0,10.] 
nt  = 100 
uhatsol = RK4(TimeSpan,uhat0,nt) 
print np.shape(uhatsol) 
print uhatsol[:6,:] 

我打印出來的拳頭6個步驟的迭代，在第6步發生的錯誤，我不明白爲什麼會這樣。 6個步驟的結果是：

nls.py:44: RuntimeWarning: overflow encountered in square 
    z = -(1j/2.) * (kx**2) * uhat + 1j * np.fft.fft((abs(u)**2) * u) 
(101, 512) 
[[ 4.02123859e+01 +0.00000000e+00j -3.90186082e+01 +3.16101312e-14j 
    3.57681095e+01 -1.43322854e-14j ..., -3.12522653e+01 +1.18074871e-13j 
    3.57681095e+01 -1.20028987e-13j -3.90186082e+01 +1.62245217e-13j] 
[ 4.02073593e+01 +2.01061092e+00j -3.90137309e+01 -1.95092228e+00j 
    3.57636385e+01 +1.78839803e+00j ..., -3.12483587e+01 -1.56260675e+00j 
    3.57636385e+01 +1.78839803e+00j -3.90137309e+01 -1.95092228e+00j] 
[ 4.01015488e+01 +4.02524105e+00j -3.89110557e+01 -3.90585271e+00j 
    3.56695007e+01 +3.58076808e+00j ..., -3.11660830e+01 -3.12911766e+00j 
    3.56695007e+01 +3.58076808e+00j -3.89110557e+01 -3.90585271e+00j] 
[ 3.98941946e+01 +6.03886019e+00j -3.87098310e+01 -5.85991079e+00j 
    3.54849686e+01 +5.37263725e+00j ..., -3.10047495e+01 -4.69562640e+00j 
    3.54849686e+01 +5.37263725e+00j -3.87098310e+01 -5.85991079e+00j] 
[ 3.95847537e+01 +8.04663227e+00j -3.84095149e+01 -7.80840256e+00j 
    3.52095058e+01 +7.15970026e+00j ..., -3.07638375e+01 -6.25837011e+00j 
    3.52095070e+01 +7.15970040e+00j -3.84095155e+01 -7.80840264e+00j] 
[ 1.47696187e+22 -7.55759947e+22j 1.47709575e+22 -7.55843420e+22j 
    1.47749677e+22 -7.56093844e+22j ..., 1.47816312e+22 -7.56511230e+22j 
    1.47749559e+22 -7.56093867e+22j 1.47709516e+22 -7.55843432e+22j]] 

在第6步，迭代的值是瘋狂的。 Aslo，溢出錯誤發生在這裏。

任何幫助？謝謝！！！！

Answer 1

第一次解析中兩個不同的錯誤是顯而易見的。

1. （發現無效蟒蛇numpy的）至於說了好幾次，的fft標準的實現不包含由尺寸縮放，這是用戶的責任。因此，對於一個矢量n組件u，

   fft(ifft(uhat)) == n*uhat and ifft(fft(u))==n*u 
   

   如果你想使用uhat = fft(u)，然後重建必須是u=ifft(uhat)/n。

2. 在RK4步驟中，您必須決定一個因子爲h的地方。或者（如例如，類似的其他人）

   k2 = f(t+0.5*h, y+0.5*h*k1) 
   

   或

   k2 = h*f(t+0.5*h, y+0.5*k1) 
   

---

然而，校正這些點僅延遲爆破。難怪有動力爆炸的可能性，這從三次方來看是可以預料的。一般來說，如果所有的項都是線性或亞線性的，那麼只能預期「緩慢」的指數增長。

爲了避免「非物理」奇點，人們必須按照Lipschitz常數成反比地縮放步長。由於這裏的Lipschitz常數的大小是u^2，所以必須動態適應。我發現在區間[0,1]中使用1000步，即h=0.001，沒有奇點。在[0,10]區間內，這仍然適用於10 000步。

---

更新沒有時間導數的原始等式的部分是自伴隨，這意味着該函數的模的平方（超過絕對值的平方積分）的確切溶液被保留。因此總體情況是一個旋轉。現在的問題是，功能的某些部分可能會以如此小的半徑「旋轉」，或者如此高的速度以至於時間步長代表了旋轉的大部分甚至多次旋轉。這很難用數值方法來捕捉，因此需要減少時間步長。這個現象的一般名稱是「僵硬的微分方程」：顯式龍格 - 庫塔方法不適用於僵硬問題。

---

UPDATE2：用人methods used before，人們可以使用

vhat = exp(0.5j * kx**2 * t) * uhat 

，其允許具有更大的步長大小的穩定溶液解決解耦頻域中的線性部分。正如在KdV方程的治療中，線性部分i*u_t+0.5*u_xx=0的DFT下解耦至

i*uhat_t-0.5*kx**2*uhat=0 

，因此可以容易地解決到相應的指數

exp(-0.5j * kx**2 * t). 

然後將該充滿方程使用變異解決常量通過設定

uhat = exp（-0.5j * kx ** 2 * t）* vhat。

這提升了kx較大部件的剛度負擔，但第三個功率仍然存在。因此，如果步長變大，數值解決方案將在很少的步驟中爆炸。下面

import numpy as np 
import math 
from matplotlib import pyplot as plt 
from matplotlib import animation 


#----- Numerical integration of ODE via fixed-step classical Runge-Kutta ----- 

def RK4Stream(odefunc,TimeSpan,uhat0,nt): 
    h = float(TimeSpan[1]-TimeSpan[0])/nt 
    print h 
    w = uhat0 
    t = TimeSpan[0] 
    while True: 
     w = RK4Step(odefunc, t, w, h) 
     t = t+h 
     yield t,w 

def RK4Step(odefunc, t,w,h): 
    k1 = odefunc(t,w) 
    k2 = odefunc(t+0.5*h, w+0.5*k1*h) 
    k3 = odefunc(t+0.5*h, w+0.5*k2*h) 
    k4 = odefunc(t+h,  w+k3*h) 
    return w + (k1+2*k2+2*k3+k4)*(h/6.) 

#----- Constructing the grid and kernel functions ----- 
L = 40 
nx = 512 
x = np.linspace(-L/2,L/2, nx+1) 
x = x[:nx] 

kx1 = np.linspace(0,nx/2-1,nx/2) 
kx2 = np.linspace(1,nx/2, nx/2) 
kx2 = -1*kx2[::-1] 
kx = (2.* np.pi/L)*np.concatenate((kx1,kx2)) 

def uhat2vhat(t,uhat): 
    return np.exp(0.5j * (kx**2) *t) * uhat 

def vhat2uhat(t,vhat): 
    return np.exp(- 0.5j * (kx**2) *t) * vhat 

#----- Define RHS ----- 
def uhatprime(t, uhat): 
    u = np.fft.ifft(uhat) 
    return - 0.5j * (kx**2) * uhat + 1j * np.fft.fft((abs(u)**2) * u) 

def vhatprime(t, vhat): 
    u = np.fft.ifft(vhat2uhat(t,vhat)) 
    return 1j * uhat2vhat(t, np.fft.fft((abs(u)**2) * u)) 

#------ Initial Conditions ----- 
u0  = 1./np.cosh(x) #+ 1./np.cosh(x+0.4*L)+1./np.cosh(x-0.4*L) #symmetric or remove jump at wrap-around 
uhat0 = np.fft.fft(u0) 

#------ Solving for ODE ----- 
t0 = 0; tf = 10.0; 
TimeSpan = [t0, tf] 
# nt  = 500 # limit case, barely stable, visible spurious bumps in phase 
nt  = 1000 # boring but stable. smaller step sizes give same picture 
vhat0 = uhat2vhat(t0,uhat0) 

fig = plt.figure() 
ax1 = plt.subplot(211,ylim=(-0.1,2)) 
ax2 = plt.subplot(212,ylim=(-3.2,3.2)) 
line1, = ax1.plot(x,u0) 
line2, = ax2.plot(x,u0*0) 

vhatstream = RK4Stream(vhatprime,[t0,tf],vhat0,nt) 

def animate(i): 
    t,vhat = vhatstream.next() 
    print t 
    u = np.fft.ifft(vhat2uhat(t,vhat)) 
    line1.set_ydata(np.real(np.abs(u))) 
    line2.set_ydata(np.real(np.angle(u))) 
    return line1,line2 

anim = animation.FuncAnimation(fig, animate, interval=15000/nt+10, blit=False) 

plt.show()