Using Kutta Merson on NLS

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

• If you write the equations of your method it would be easier to follow you. – nicoguaro Jul 28 '19 at 14:43
• The equations of kutta-merson are in the hyperlink, in the first paragraph, but any way here it is encyclopediaofmath.org/index.php/Kutta-Merson_method – Matheus Fachini Jul 28 '19 at 20:27
• It is more probable that somebody help you if you include the details in your question instead of needing to go read somewhere else. – nicoguaro Jul 28 '19 at 20:34
• What is your time interval, and what is the expected shape of the solution vs. the shape that you obtained? – Lutz Lehmann Aug 6 '19 at 14:29

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