Python中RK4算法中的错误 [英] Error in RK4 algorithm in Python

问题描述

我正在求解非线性薛定inger(NLS)方程:

I am solving a nonlinear Schrodinger (NLS) equation:

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

应用傅立叶变换后，它变为:

after applying Fourier Transform, it becomes:

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

其中，uhat是u的傅立叶变换.上面的等式(2)是一个相当明确的IVP，可以通过第四阶Runge-Kutta方法求解.这是我用于求解方程式(2)的代码:

where uhat is the Fourier Transform of u. The equation (2) above is a pretty stated IVP, which can be solved by the 4th oder Runge-Kutta method. Here are my code for solving equation (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步发生了，我不明白为什么会这样.这六个步骤的结果是:

I printed out the fist 6 steps of the iteration, the error occurred at the 6th step, I don't understand why this happened. The results of the 6 steps are:

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]]

在第六步，迭代的值是疯狂的. Aslo，这里发生了溢出错误.

At the 6th step, the values of the iteration are crazy. Aslo, the overflow error occurred here.

有帮助吗??谢谢！！！

Any help?? Thank you!!!!

推荐答案

在第一次解析中，明显出现了两个不同的错误.

Two different errors were obvious in the first parsing.

1. (发现对于python numpy无效)多次告诉我们，fft的标准实现不包含按维度进行缩放，这是用户.因此，对于具有n个分量的向量u，

1. (found to be not valid for python numpy) As told several times, the standard implementations of fft do not contain the scaling by the dimension, this is the responsibility of the user. Thus for a vector u of n components,

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

如果要使用uhat = fft(u)，则重构必须为u=ifft(uhat)/n.

If you want to use uhat = fft(u), then the reconstruction has to be u=ifft(uhat)/n.

在RK4步骤中，您必须为系数h确定一个位置.任何一个(例如，其他类似地)

In the RK4 step, you have to decide on one place for the factor h. Either (as example, the others analogously)

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

或

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

---

但是，更正这些点只会延迟爆炸.发生动态爆炸的可能性也就不足为奇了，这是可以预料的.通常，只有当所有项都是线性或次线性时，才可以期望缓慢的"指数增长.

---

However, correcting these points only delays the blowup. That there is the possibility for dynamical blow-up is no wonder, it is to be expected from the cubic term. In general one can only expect "slow" exponential growth if all terms are linear or sub-linear.

为避免非自然的"奇异性，必须将步长的大小与Lipschitz常数成反比.由于此处的Lipschitz常数的大小为u^2，因此必须动态调整.我发现在间隔[0,1]中使用1000步，即h=0.001时，没有奇异之处.对于间隔[0,10]的10000步仍然适用.

To avoid "unphysical" singularities, one has to scale the step size inversely proportional to the Lipschitz constant. Since the Lipschitz constant here is of size u^2, one has to dynamically adapt. I found that using 1000 steps in the interval [0,1], i.e., h=0.001, proceeds without singularity. This still holds true for 10 000 steps on the interval [0,10].

更新原始方程式中没有时间导数的部分是自伴随的，这意味着函数的范数平方(绝对值平方的整数)被精确保留.解决方案.因此，一般情况是旋转.现在的问题是，函数的某些部分可能以很小的半径或很高的速度旋转"，以至于时间步长代表了旋转的大部分，甚至是多次旋转.这很难用数值方法来捕获，因此需要减少时间步长.这种现象的总称是刚性微分方程":明确的Runge-Kutta方法不适用于刚性问题.

Update The part of the original equation without the time derivative is self-adjoint, which means that the norm square of the function (integral over the square of the absolute value) is preserved in the exact solution. The general picture is thus a rotation. The problem now is that parts of the function might "rotate" with such a small radius or such a high velocity that a time step represents a large part of a rotation or even multiple rotations. This is hard to capture with numerical methods, which requires thus a reduction in the time step. The general name for this phenomenon is "stiff differential equation": Explicit Runge-Kutta methods are ill-suited for stiff problems.

Update2 :使用以前使用的方法，可以解决解耦中的线性部分频域使用

Update2: Employing methods used before, one can solve the linear part in the decoupled frequency domain using

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

允许使用较大步长的稳定解决方案.与处理KdV方程一样，线性部分i*u_t+0.5*u_xx=0在DFT下解耦至

which allows for a stable solution with larger step sizes. As in the treatment of the KdV equation, the linear part i*u_t+0.5*u_xx=0 decouples under the DFT to

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

因此可以轻松地求解成相应的指数

and can thus be easily solved into the corresponding exponentials

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

然后通过设置

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

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

这为kx的较大组件减轻了一些刚度负担，但仍然保留了三次方.因此，如果步长变大，则数值解将以很少的步数爆炸.

This lifts some of the burden of the stiffness for the larger components of kx, but still, the third power remains. So if the step size gets to large, the numerical solution explodes in very few steps.

下面的工作代码

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()

这篇关于Python中RK4算法中的错误的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！