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

[John B. Thoo]

On May 16, 2009, at 2:35 AM, Francesco Potorti` wrote:

> > Now, my modified script (appended below) has several "for loops," and
> > I've read here on more than one occasion that it helps to vectorize
> > "for loops" to speed things up.  However, I can't see how to
> > vectorize the "for loops" in my script.  I would appreciate help with
> > this from anyone.
>
> By having a quick look, I'd say that you should be able to vectorise
> much of the code.
>
> [... snip, snip ...]

Thank you again.

Now my problem is that my code does not seem to be producing what I  
expect.  Because I'm new to Octave and coding and numerics, I don't  
know if my problem lies with my code, my math, or both.  I would  
appreciate any help from anyone.

Here is a recap of what I'm trying to accomplish.  My code is  
appended at the bottom.

For practice, I'm trying to solve the inviscid Burgers's equation,   
u_t + u*u_x = 0,  by first Fourier transforming the equation and  
initial conditions.  Writing the equation as

  (BE)  u_t = - (u^2/2)_x,  u(x,0) = u0,

I apply the F.T. and then solve the system of ODEs

  (d/dt)U(k,t) = -ik/2 * sum_m {U(m,t)*U(k-m,t)},  U(k,0) = U0,

where  U = F.T. of u  and the sum is over  m = -infinity..infinity.   
Finally, the solution of (BE) should be

  u(x,t) = sum_k {U(k,t)*exp(ikx)},

where the sum is over  k = -infinity..infinity.

Now, when I use the initial condition

  u(x,0) = 1 if x <= 0, 0 otherwise,

for example, I don't recover the initial condition when I plot the  
solution with  t = 0, and I don't understand why.  In my code, I use   
U(-k) = U(k)',  the complex conjugate.

Again, I don't know if my problem lies with my code, my math, or  
both, and I would appreciate any help from anyone.  TIA.

---John.

------------begin code------------
t0 = 0;   % initial time
tf = 10;   % final time
a = -1;   % left end point
b = 1;   % right end point
M = 110;   % number of time intervals
N = 100;   % number of space intervals
dt = (tf - t0)/M;
t = t0:dt:tf;
dx = (b - a)/N;
xx = a:dx:b;


global K;   % truncation index
K = 3;

% define initial condition  zz = u(x,0)
zz = (xx <= 0);   % step function: 1 on L, 0 on R
%zz = exp (-(xx.^2)./(0.1^2));   % Gaussian hump
fftzz = fftshift (fft (zz));
c = floor (length (fftzz)/2) + 1;
z0 = fftzz (c-K:c+K);
T = linspace (t0, tf, M+1)';


function w = cplxfcn(z,t)
global K;

  U1 = zeros (K, 2*K+1);
  U0 = zeros (1, 2*K+1);
  U2 = zeros (K, 2*K+1);
  U = zeros (2*K+1, 2*K+1);
  V = zeros (2*K+1, 1);

  for k = 1:K
    U1(k, 1:k) = z(k:-1:1);
    U1(k, k+1:k+K) = z(2:K+1)';
    U2(k, k+1:k+K+1) = z(K+1:-1:1);
  endfor

  for k = 1:K-1
    U2(k, k+K+2:2*K+1) = z(2:K-k+1)';
  endfor

  U = [U1; U0; U2];

  V(1:K) = z(K+1:-1:2)';
  V(K+1) = z(1);
  V(K+2:2*K+1) = z(2:K+1);

  UV = U*V;
  w(1:K) = i*(K:-1:1)./2.*UV.'(1:K);
  w(K+1) = UV.'(K+1);
  w(K+2:2*K+1) = -i*(1:K)./2.*UV.'(K+2:2*K+1);
endfunction


global n;
n = length(z0);

function y = realfcn(x,t)
global n;

 z = x(1:n) + i*x(n+1:2*n);
 w = cplxfcn(z,t);
 y = [real(w),imag(w)];
endfunction


% Each row of Z is a time slice of solution with time values from T.
[X, istate, msg] = lsode ('realfcn',[real(z0),imag(z0)], T);
Z = X(:,1:n) + i*X(:,n+1:2*n);


% recover  u
u = zeros (N+1, M+1);
EXP = zeros (2*K+1, N+1);

for k = 1:2*K+1
  EXP(k, 1:N+1) = exp (i*(k-K-1)*xx(1:N+1));
endfor

u = (Z*EXP).';


% Plot the solution
plot(xx, real (u(:,1))) %, 'r', xx, real (u(:,2)), 'g', xx, real (u(:, 
3)), 'b')
------------end code------------