Re: Trying to solve du/dt = i*u

[John B. Thoo]

Hello, Thomas and everyone.

On May 13, 2009, at 8:53 PM, Thomas Shores wrote:

> I'm not going to dissect your code, but I will include a template  
> script that reduces your problem to making a well-formed definition  
> of the complex rhs f(z,t).  It's at the bottom of this email, and  
> if you remark/deremark the appropriate lines, you'll get each of  
> the first two examples you proposed.
>
> However, I have a few comments about your method, which is a  
> separation of variables type argument with complex exponentials  
> whose time-dependent coefficients are just the Fourier coefficients  
> of the solution:
>
> 1. There is an immediate problem with truncation m=-M..M of the  
> double series in the indices k, m which comes originally from a  
> double series in k and ell, with k+ell=m. The problem is that for  
> all indices k except k=0, you will be missing some terms with index  
> k-m, which compounds the mathematical truncation error.  It also  
> means you have to be careful about limits when you calculate the  
> sums involved in the rhs.
>
> 2.  There is a question as to why you are using a Fourier method at  
> all for Burgers equation.  This nonlinear hyperbolic problem will  
> generate shocks in the solution at a minimum time T = -min u_0(x),  
> where u_0(x) = u(x,0) is the initial condition and the derivative   
> u_0(x)'  is somewhere negative.  At this time, the problem becomes  
> ill-posed and Fourier expansions are not helpful.  It is this very  
> difficulty that is the starting point for the modern theory of  
> "numerical methods for conservation laws" (see, e.g., LeVeque's  
> monograph of the same name.)

I don't know if this is of any continuing interest, but I want to  
bring it (almost) to a close.

After longer than it should have taken, I've finally written a short  
code to solve  u_t + f'(u)*u_x = 0  using a "pseudospectral" method  
avoiding the immediate problem with truncation m=-M..M of the double  
series in the indices k, m.  (I understand your objections to using a  
Fourier method, but that's what I've been directed to do, so I'm  
giving it my best.)

Here is the my code.  I would appreciate any comments you may have  
about it as I'm still very much learning.  A few notes:

1) I had to use  .'  in a couple of places where I hadn't expected,  
namely, in the lines

  u_x = ifft (ifftshift (i*k.*fftshift (fft (u)).'));

  u_t = real (-fprime (u).*u_x.');


Obviously I still don't understand when something is outputted as a  
column vector or a row vector.  Is there a rule of thumb that can  
guide me?

2) The solution does not seem to travel a distance one in time one  
when I set  f'(u) = 1.  Should this be a red flag?

%%------------begin code---------------

% This script solves the pde  u_t + f'(u)*u_x = 0  using a  
"pseudospectral"
% method.


% set time and space parameters
t0 = 0;   % initial time
tf = 1;   % final time
L = -1;   % left end point
R = 1;   % right end point
M = 10;   % number of time intervals
N = 100;   % number of space intervals

dt = (tf - t0)/M;   % time mesh
t = t0:dt:tf;   % time interval
dx = (R - L)/N;   % space mesh
x = a:dx:b;   % space interval


% define initial condition  u0 = u(x,0)
%u0 = (x <= 0);   % step function: 1 on L, 0 on R
u0 = exp (-25*x.^2);   % Gaussian hump

global lenu0;
lenu0 = length (u0);


% define the equation
function speed = fprime (u)
  speed = u;   % set  fprime(u)  in  u_t + fprime(u)*u_x = 0
endfunction

function u_t = F (u,t)
global lenu0;

  speed = fprime (u);

  c = floor (lenu0/2);
  k = -c:c;

  u_x = ifft (ifftshift (i*k.*fftshift (fft (u)).'));

  u_t = real (-fprime (u).*u_x.');
endfunction


% solve the equation
u = lsode("F", u0, t).';


% plot the solution
clf

%% static plot
%plot(x, u(:,1), 'r', x, u(:,2), 'g', x, u(:,3), 'b')

%% animation
wait = max (1/(tf - t0 + 1), 0.05);   % pause time between frames

for j = 1:M+1
  plot (x, u(:, j), '-o');
  title (strcat ('time =', num2str (t(j))));
%  axis ([-1, 1, -1, 1]);
  pause(wait);
endfor

%%------------end code---------------

Thanks for your help.

---John.