ode - Error in RK4 algorithm in Python -

i solving nonlinear schrodinger (nls) equation:

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

after applying fourier transform, becomes:

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

where uhat fourier transform of u. equation (2) above pretty stated ivp, can solved 4th oder runge-kutta method. here code solving equation (2):

import numpy np import math matplotlib import pyplot plt 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 nx matrix     t[0]   = timespan[0]         w[0,:] = uhat0          # enter initial conditions in w     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 grid , 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 ode ----- timespan = [0,10.] nt       = 100 uhatsol  = rk4(timespan,uhat0,nt)  print np.shape(uhatsol) print uhatsol[:6,:] 

i printed out fist 6 steps of iteration, error occurred @ 6th step, don't understand why happened. results of 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]] 

at 6th step, values of iteration crazy. aslo, overflow error occurred here.

any help?? thank you!!!!

two different errors obvious in first parsing.

1. (found not valid python numpy) told several times, standard implementations of fft not contain scaling dimension, responsibility of user. vector u of n components,

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

   if want use uhat = fft(u), reconstruction has u=ifft(uhat)/n.

2. in rk4 step, have decide on 1 place factor h. either (as example, others analogously)

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

   or

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

---

however, correcting these points delays blowup. there possibility dynamical blow-up no wonder, expected cubic term. in general 1 can expect "slow" exponential growth if terms linear or sub-linear.

to avoid "unphysical" singularities, 1 has scale step size inversely proportional lipschitz constant. since lipschitz constant here of size u^2, 1 has dynamically adapt. found using 1000 steps in interval [0,1], i.e., h=0.001, proceeds without singularity. still holds true 10 000 steps on interval [0,10].

---

update part of original equation without time derivative self-adjoint, means norm square of function (integral on square of absolute value) preserved in exact solution. general picture rotation. problem parts of function might "rotate" such small radius or such high velocity time step represents large part of rotation or multiple rotations. hard capture numerical methods, requires reduction in time step. general name phenomenon "stiff differential equation": explicit runge-kutta methods ill-suited stiff problems.

---

update2: employing methods used before, 1 can solve linear part in decoupled frequency domain using

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

which allows stable solution larger step sizes. in treatment of kdv equation, linear part i*u_t+0.5*u_xx=0 decouples under dft

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

and can solved corresponding exponentials

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

the full equation tackled using variation of constants setting

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

this lifts of burden of stiffness larger components of kx, still, third power remains. if step size gets large, numerical solution explodes in few steps.

working code below

import numpy np import math matplotlib import pyplot plt 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 grid , 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 @ wrap-around uhat0   = np.fft.fft(u0)  #------ solving ode ----- t0 = 0; tf = 10.0; timespan = [t0, tf] # nt       = 500 # limit case, barely stable, visible spurious bumps in phase nt       = 1000 # boring  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()