1

I am trying to solve the solution of a 1D KdV equation (ut+uux+uxxx=0) starting from a two solitons initial condition using Fourier spectral method. My code blows up the solution for a reason that I cannot figure out. I would appreciate any help. This is my main code file:

%Spatial variable on x direction
L=4; %domain on x
delta=0.05; %spatial step size
xmin=-L; %minimum boundary
xmax=L; %maximum boundary 
N=(xmax-xmin)/delta; %number of spatial points
x=linspace(xmin,xmax,N); %spatial vector

% 1D Initial state
sigma1 = 2; sigma2 = 4;
c1     = 2;     c2 = 1;
h1     = 1;     h2 = 0.3;
U = h1*sech(sigma2*(x+c1)).^2+h2*sech(sigma1*(x-c2)).^2;

%Fast Fourier Transform to the initial condition                 
Ut = fftshift(fft(U));
Ut = reshape(Ut,N,1);
%1D Wave vector disretisation 
k = (2*pi/L)*[0:(N/2-1) (-N/2):-1];  
k(1) = 10^(-6);
k = fftshift(k);
k = reshape(k,N,1);

%first derivative (advection)
duhat = 1i *k .*Ut;
du = real(ifft(ifftshift(duhat))); %inverse of FT
%third derivative (diffusion)
ddduhat = -1i * (k.^3) .*Ut;
dddu = real(ifft(ifftshift(ddduhat))); %inverse of FT

% Time variable
dt = 0.1; %time step
tspan = [0 4];

%solve
Time = 50; 
for TimeIteration = 1:2:Time
t= TimeIteration * dt; 
[Time,Sol] = ode45('FFT_rhs_1D',tspan,U,[],du,dddu);
Sol = Sol(TimeIteration,:);
%plotting
plot(x,abs(Sol),'b','LineWidth',2);
end

And this is the function that solves the equation:

function rhs = FFT_rhs_1D(tspan,U,dummy,du,dddu)
%solve the right hand side
rhs = - U .* du - dddu;
end

Thank you.

User123
  • 29
  • 9
  • What is the purpose of the for loop? Does this code even run? (`FFT_rhs_1D` has too many input arguments for `ode45`) – David Jun 30 '20 at 23:40
  • The loop shows you how the solution changes with time. You can surely remove it but then we need to choose a specific time to see the resulting solution. The function file does not run, it is being called by the main file. The ode45 file has the main inputs; initial condition and the time span, in addition to the two derivatives. All of these inputs are necessary for the right hand side of the kdv equation to be evaluated. – User123 Jul 01 '20 at 00:12
  • I am thinking that the reason behind why the code is not giving a stable solution is because of the higher order numerical stiffness, but I don't know what should be done in this situation or what line in the code that needs to be modified/removed! – User123 Jul 01 '20 at 00:18
  • @Lama Did you ever resolve this? I also have trouble with blowups, even with just a single soliton. – user196574 Feb 07 '22 at 08:59
  • The (first?) problem here is that the space derivative are only computed once for the initial state and then kept constant. So essentially you are solving a system of decoupled scalar equations `dU(i)/dt = a(i)*U(i)+b(i)` and convergence depends on the sign of `a(i)=du(i)`. Of course, this solution is completely unrelated to the original PDE. /// Move the Fourier transform into the ODE function. Perhaps even employ operator splitting, that is, apply an integrating factor that combines the linear terms, see https://stackoverflow.com/questions/29617089 for a local discussion. – Lutz Lehmann Feb 28 '22 at 18:29

0 Answers0