Blow-up from PDE or aliasing error? How to know and check?

81 Views Asked by At

I have this PDE: $u_{t}= u_{xx}+g(t)\cdot u_{x}^{2}$ (with periodic boundary conditions in the interval $(-\pi,\pi)$)

that I don't have an exact solution for or know if it blows up. However, by using the pseudospectral I get a blow up, : enter image description here

In general, can I check numerically if a blow-up comes from the PDE itself or the cumulative error from picking a very small timestep? Maybe by checking the power-spectra for various time-steps (both small and big)

       function [u,x,t,kk,amp] = pseudo_spec_ppde(u0,t_0,t_f,M,N)
       % Mesh in space
     dx = 2*pi/N;
     x = -pi:dx:pi-dx;
    x = x';
   kk = 0:1:(N/2-1);

   %we set high beta
     beta=100;

   % define the mesh in time
    dt = (t_f-t_0)/M;
    t = t_0:dt:t_f;

     u(:,1) = u0;
  [kk,amp(:,1)] = find_spec(u0);

        U_now = fft(u(:,1));
 for j=1:M
 %approximate the convolution sum
 G = convolve2(U_now);
 % initialize the zero frequency mode
 U_new(1) = U_now(1) + dt*G(1)*normcdf((beta-t(j))/beta,0.5,0.2);
 % the nonzero modes
 for k = 1:floor(N/2)-1
    U_new(k+1) = (1 - k^2*dt)*U_now(k+1) + dt*G(k+1)*normcdf((beta-   t(j))/beta,0.5,0.2);
    U_new(N-k+1) = (1 - k^2*dt)*U_now(N-k+1) + dt*G(N-k+1)*normcdf((beta-    t(j))/beta,0.5,0.2);
  end
  % update U_now for the next solution column
  U_now = U_new;
 %return to physical space
  u(:,j+1) = real(ifft(U_new));

  [kk,amp(:,j+1)] = find_spec(u(:,j+1));
 end

With the subroutine

        function V = convolve2(u)
       N = length(u);

       U = zeros(1,N);
     %we zero-pad the middle to improve accuracy of the frequencies and thus         of the IFFT

     U(1:floor(N/2)) = u(1:floor(N/2));
    U(N-floor(N/2)+1:N) = u(floor(N/2)+1:N);

   V(1) = 0;
   for k=1:floor(N/2)-1
 V(k+1) = i*k*U(k+1);
  V(N-k+1) = -i*k*U(N-k+1);
   end
 %we bring each u_x term to physical space
  V1 = ifft(V);
 V2 = ifft(V);
 v = V1.*V2;
  %we bring their product to the frequency space 
 v = fft(v);
  clear V
 V(1) = v(1);
 for k=1:floor(N/2)-1
 V(k+1) = v(k+1);
V(N-k+1) = v(N-k+1);
end
V(floor(N/2)+1) = 0;