All geek questions in one place

Error in RK4 algorithm in Python

I solve the nonlinear Schrödinger equation (NLS):

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

after applying the Fourier transform, it becomes:

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

where uhatis the Fourier transform u. Equation (2) above is a fairly stated IVP, which can be solved by the fourth Runge-Kutta method. Here is 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,:]

I printed 6 iteration steps, the error in the sixth step, I do not understand why this happened. Results 6 steps:

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

In step 6, the iteration values ​​are crazy. In addition, an overflow error occurred here.

Any help ?? Thank!!!!

+1
source share
1 answer

In the first analysis, two different errors were apparent.

  • (, python numpy) , fft , . , u n,

    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)

. , , , . "" , .

"" , , . u^2, . , 1000 [0,1], .. h=0.001, . 10 000 [0,10].

. , , ( ) . , . , "" , . , , , . - " ": - .

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

Source: https://habr.com/ru/post/1672558/

More articles:


All Articles