Error en el algoritmo RK4 en Python

Estoy resolviendo una ecuación de Schrodinger no lineal (NLS):

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

después de aplicar la Transformada de Fourier, se convierte en:

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

donde uhat es la Transformada de Fourier de u . La ecuación (2) anterior es un IVP bastante bien establecido, que se puede resolver con el método de Runge-Kutta del cuarto oder. Aquí está mi código para resolver la ecuación (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 ht = 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]+hw[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,:] 

Imprimí los primeros seis pasos de la iteración, el error ocurrió en el sexto paso, no entiendo por qué sucedió esto. Los resultados de los 6 pasos son:

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

En el paso 6, los valores de la iteración son una locura. Aslo, el error de desbordamiento se produjo aquí.

¿¿Alguna ayuda?? ¡¡¡¡Gracias!!!!

Dos errores diferentes fueron evidentes en el primer análisis.

  1. (Se encontró que no es válido para el número de python) Como se dijo varias veces, las implementaciones estándar de fft no contienen la escala por dimensión, esto es responsabilidad del usuario. Así, para un vector u de n componentes,

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

    Si desea utilizar uhat = fft(u) , entonces la reconstrucción debe ser u=ifft(uhat)/n .

  2. En el paso RK4, tienes que decidir en un lugar para el factor h . O bien (como ejemplo, los otros de forma análoga)

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

    o

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

Sin embargo, corregir estos puntos solo retrasa la explosión. Que exista la posibilidad de una explosión dinámica no es una maravilla, se debe esperar del término cúbico. En general, solo se puede esperar un crecimiento exponencial “lento” si todos los términos son lineales o sub-lineales.

Para evitar las singularidades “no físicas”, uno tiene que escalar el tamaño del paso inversamente proporcional a la constante de Lipschitz. Como la constante de Lipschitz aquí es de tamaño u^2 , uno tiene que adaptarse dinámicamente. Encontré que el uso de 1000 pasos en el intervalo [0,1], es decir, h=0.001 , procede sin singularidad. Esto sigue siendo válido para 10 000 pasos en el intervalo [0,10].


Actualización La parte de la ecuación original sin la derivada de tiempo es autoadjuntiva, lo que significa que el cuadrado de la norma de la función (integral sobre el cuadrado del valor absoluto) se conserva en la solución exacta. El cuadro general es así una rotación. El problema ahora es que las partes de la función podrían “girar” con un radio tan pequeño o una velocidad tan alta que un paso de tiempo representa una gran parte de una rotación o incluso múltiples rotaciones. Esto es difícil de capturar con métodos numéricos, lo que requiere una reducción en el paso del tiempo. El nombre general de este fenómeno es “ecuación diferencial rígida”: los métodos explícitos de Runge-Kutta no son adecuados para problemas rígidos.


Update2: Empleando métodos utilizados anteriormente , se puede resolver la parte lineal en el dominio de frecuencia desacoplada usando

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

lo que permite una solución estable con pasos más grandes. Al igual que en el tratamiento de la ecuación de KdV, la parte lineal i*u_t+0.5*u_xx=0 desacopla bajo la DFT para

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

y por lo tanto puede ser fácilmente resuelto en los exponenciales correspondientes

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

La ecuación completa se aborda utilizando la variación de las constantes estableciendo

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

Esto levanta parte de la carga de la rigidez para los componentes más grandes de kx , pero aún así, la tercera potencia permanece. Por lo tanto, si el tamaño de paso aumenta, la solución numérica explota en muy pocos pasos.

Código de trabajo a continuación

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