1
$\begingroup$

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

Improve this question
$\endgroup$
  • $\begingroup$ If you write the equations of your method it would be easier to follow you. $\endgroup$ – nicoguaro♦ Jul 28 '19 at 14:43
  • $\begingroup$ 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 $\endgroup$ – Matheus Fachini Jul 28 '19 at 20:27
  • $\begingroup$ It is more probable that somebody help you if you include the details in your question instead of needing to go read somewhere else. $\endgroup$ – nicoguaro♦ Jul 28 '19 at 20:34
  • $\begingroup$ What is your time interval, and what is the expected shape of the solution vs. the shape that you obtained? $\endgroup$ – Lutz Lehmann Aug 6 '19 at 14:29
1
$\begingroup$

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

enter image description here

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.