错误在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()