Python中RK4算法中的错误

Python中RK4算法中的错误

在第一次解析中，显然有两个不同的错误。

1. （发现对于python numpy无效） 多次告知，的标准实现

   fft

   不包含按维度进行缩放，这是用户的责任。因此，对于矢量

   u

   的

   n

   分量，

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

如果要使用

uhat = fft(u)

，则必须进行重建

u=ifft(uhat)/n

。

1. 在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]的10000步仍然适用。

---

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

---

Update2：
使用以前
使用的方法，可以使用以下方法解耦去耦频域中的线性部分

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 npimport mathfrom matplotlib import pyplot as pltfrom 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,wdef 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   = 40nx  = 512x   = 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) * uhatdef 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-arounduhat0   = 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 phasent       = 1000 # boring  but stable. smaller step sizes give same picturevhat0 = 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,line2anim = animation.FuncAnimation(fig, animate, interval=15000/nt+10, blit=False)plt.show()