Using Kutta Merson on NLS

Question

I'm trying to use the Kutta-Merson to get the same results as in the book Solitons, Nonlinear Evolution Equations and Inverse Scattering - M. J. Ablowitz - pg 140

The author propose using the Kutta-Merson in the following scheme:

$$ i\frac{du_j}{dt} + (u_{j-1} - 2u_j + u_{j+1})/h^2 + |u_j|^2(u_{j-1} + u_{j+1}) = 0 $$

given condition of periodicity and the following initial condition

$$ u(x,0) = 0.5 + 0.05 \cos(\mu x) + 10^{-5} i \sin(\mu x) $$ where $L = 2\pi\sqrt{2}$ is the length of the interval and $\mu = 2\pi/L$

I'm not very familiar with this method for PDE, but here is what i try to do, using scilab:

I define the laplacian and the nonlinear using scilab functions, and put the scheme in a form of Kutta's method

$$\frac{du_j}{dt} = F(u_j)$$

I use the condition of periodicity to define F. Here is the Laplacian:

function l =L(u)
    A = zeros(M+1,M+1)
    // Condição periódica imposta ao laplaciano
    A(1,1) = -2; A(1,2) = 1; A(1,M+1) = 1
    for j=2:M
        A(j,j-1) = 1; A(j,j) = -2; A(j,j+1) = 1
    end
    A(M+1,1) = 1; A(M+1,M) = 1; A(M+1,M+1) = -2
    A = h^(-2)*A

    l = A*u
endfunction

And here is the nonlinear therm

function n = N(u)
    n(1) = abs(u(1))^2 * (u(M+1) + u(2))
    for j=2:M
        n(j) = abs(u(j))^2 *(u(j-1) + u(j+1))
    end
    n(M+1) = abs(u(M+1))^2 * (u(M) + u(1))
endfunction

Then F becomes:

function f = F(u)
    f = %i*(L(u) + N(u)) 
endfunction

Finally I use this to calculate the solution

for n=1:Nt

    k1 = k*F(U(:,n))
    k2 = k*F(U(:,n) + 1/3*k1)
    k3 = k*F(U(:,n) + 1/6*k1 + 1/6*k2)
    k4 = k*F(U(:,n) + 1/8*k1 + 3/8*k3)
    k5 = k*F(U(:,n) + 1/2*k1 - 3/2*k3 + 2*k4)

    U(:,n+1) = U(:,n) + 1/6*k1 + 2/3*k4 + 1/6*k5
end

I'm using 5000 points at time grid and 20 for x-grid. I really don't know nothing about so if I'm doing some weird thing let me know and recommend some material about it please

Answer 1

It makes no sense to compute the matrix A in every call of F, indeed you need not construct this matrix at all. Use vector operations

function l =Lin(u)
  l = -2*u;
  l(1:M) += u(2:M+1); l(M+1) += u(1);
  l(2:M+1) += u(1:M); l(1) += u(M+1);
  l *= h^(-2);
endfunction

Similarly you can vectorize the non-linear part

function n = NonLin(u)
    n = 0*u;
    n(1)=u(M+1)+u(2);
    n(2:M)= u(1:M-1)+u(3:M+1);
    n(M+1)=u(M)+u(1);
    n = n.*abs(u).^2
endfunction

function f = F(u)
  f = 1i.*(Lin(u) + NonLin(u));
endfunction 

Using M=50 and k=0.02 and the other data as described, I get as output

There appears to be an oscillation back to the original curve from time 5 to 10