Problem with solving the linear mKdV equation using Spectral methood for spartial derivatives and Backward euler for time derivatives

Question

The mKdV equation: $$ u_{t}=-\varepsilon c^{2}u_{x}-u_{xxx} \\ u(0,x)=u_{0}(x) $$ and P.C. $[-2\pi,2\pi]$

I use spectral method for the derivatives $u_{x}$ and $u_{xxx}$: \begin{align} u_{t}&=-\varepsilon c^{2}u_{x}-u_{xxx} \\ \hat{u}_{t}&=-\varepsilon c^{2}\widehat{u_{x}}-\widehat{u_{xxx}} \\ \hat{u}_{t}(t,j)&=-i\varepsilon c^{2}j\hat{u}(t,j)+ij^{3}\hat{u}(t,j). \end{align} And now with Backward Euler for time derivative: $$ \frac{\hat{u}^{n+1}-\hat{u}^{n}}{h} =-i\varepsilon c^{2}j\hat{u}^{n+1}+ij^{3}\hat{u}^{n+1} $$ After algebraic operation I get: $$ \hat{u}^{n+1}=\frac{\hat{u}^{n}}{1-hij^{3}+hi\varepsilon c^{2}j} $$ If i set $\varepsilon=1,c=1$, it's seems the Matlab code is only good for solutions that are independent on time like $\sin(x)$ but for $\sin(\frac{1}{2}x-\frac{3}{8}t)$ the numerical solution is not good- my numerical solution seems to be correct in space but the change in time almost not exists in it (compering to the analytic solution). The implementation is

k = 0.1; %space step size
h = 0.1; %time step size
beta = 4*k^3/h; 

N = round(4*pi/k); % Number of space steps
x = linspace(-2*pi, 2*pi, N+1)';
x(end)=[];

tF=3;
M=tF/h;
NT=M; 
m = N/2;
J=fftshift(x/(4*pi));
%J =k*fftshift(-m:m-1)';
t = h*(0:NT)';

epsilon = 1; %Parameter
c = 1; %Parameter

%initial conditions
f = sin(0.5*x);
Uhat(:,1) = fft(f);
I = diag(ones(N,1));
g=1-(h*1j)*(J.^3)+(h*epsilon*(c^2)*1j*J); %BE
for i = 1:NT
    Uhat(:,i+1) = Uhat(:,i)./g;
end
for i = 1:NT+1
    U(:,i) = real(ifft(Uhat(:,i)));
end
for i=1:NT+1;
    for j=1:N;
      Uanalitic(j,i) = sin(0.5*x(j)-(3/8)*t(i));
    end
end
Err = abs(U - Uanalitic);

Thanks for the help :-)

Answer 1

For a traveling solution $u(t,x)=f(x+bt)$ you get the equation $$ (b+εc^2)f'(s)+f'''(s)=0. $$ Now in your example you set $ε=c=1$, so that $f(s)=A+B\sin\left((\sqrt{b+1})\,s+C\right)$ is the solution formula. To get that $4\pi$ periodic, you can set $\sqrt{b+1}=\frac12$, then indeed $b=-\frac34$.

---

The general $4\pi$ periodic solution has the Fourier series $$ u(t,x)=\sum_{j\in\Bbb Z} \hat u_j(t)e^{ijx/2} $$ Then $$ u_x(t,x)=i\sum_{j\in\Bbb Z} \frac{j}2\hat u_j(t)e^{ijx/2} \\ u_{xxx}(t,x)=-i\sum_{j\in\Bbb Z} \frac{j^3}8\hat u_j(t)e^{ijx/2} $$

To get the $J$ array that actually corresponds to these multiplicands none of your two attempts are correct, you need

m = round(2*pi/k); 
N = 2*m;         % Number of space steps
k = 4*pi / N;
x = linspace(-2*pi, 2*pi, N+1)';
x(end)=[];
J = 0.5*fftshift(-m:m-1)'

As the initial solution is solely in the basis frequency, one can expect that all the steps in the spectral method will replicate this behavior up to floating point rounding errors, so that $$u^n(x)=A^n\sin(x/2)+B^n\cos(x/2).$$ Inserting into the implicit Euler equation gives $$ A^n\sin(x/2)+B^n\cos(x/2)=(A^{n+1}-\tfrac38hB^{N+1})\sin(x/2)+(B^{n+1}+\tfrac38hA^{n+1})\cos(x/2) $$ so that \begin{align} A^{n+1}+iB^{n+1}&=(1+i\tfrac38h)^{-1}(A^n+iB^n) \\ &=e^{-i\tfrac38h-\frac{9}{2^7}h^2+i\frac{9}{2^9}h^3+\frac{3^4}{2^{14}}h^4+O(h^4)}(A^n+iB^n) \end{align} This means that the phase shift is correct up to $O(h^2)$, but the solution decreases like $e^{-\frac{9}{128}ht}$ which would have to be taken into account. Thus the expectes solution is closer to $$ e^{-\frac{9}{128}ht}\sin\left(\frac{x}2-\frac{3t}8+\frac{9h^2t}{512}\right) $$ Correspondingly the error is reduced dramatically by comparing to

for i=1:NT+1;
  Uanalitic(:,i) = exp(-9*t(i)*h/128).*sin(0.5*x-(3/8-9*h^2/512)*t(i));
end