RK4程序c语言,Python中RK4算法出错

在第一次解析中,两个不同的错误是显而易见的.

>(发现对python numpy无效)正如多次告诉的那样,fft的标准实现不包含维度的缩放,这是用户的责任.因此,对于n个分量的向量u,

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

如果你想使用uhat = fft(u),那么重建必须是u = ifft(uhat)/ n.

>在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步,这仍然适用.

更新原始方程中没有时间导数的部分是自伴的,这意味着函数的范数平方(绝对值的平方上的积分)保留在精确解中.因此,一般情况是轮换.现在的问题是功能的某些部分可能以如此小的半径或如此高的速度“旋转”,使得时间步长代表旋转的大部分或甚至多次旋转.这很难用数值方法捕获,因此需要减少时间步长.这种现象的一般名称是“刚性微分方程”：显式Runge-Kutta方法不适合刚性问题.

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