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

[John B. Thoo]

Hello, Thomas.

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

> %-------------begin script------------------
> % Script for solving dz/dt = f(z,t), z(t0) = z0
> % Here t0 is first coordinate of time node vector T
> % at which solution is evaluated.
> % User defines f as cplxfcn below as well as
> % initial condition and time node vector.
> % Vectors z, z0, f(z,t) must be row vectors.
>
> %z0 = 1;
> z0 = [1+i,1-i,-2]; % initial condition as row vector
> %T = linspace(0,2*pi,20);
> T = linspace(0,2,100); % time node vector
>
> function w = cplxfcn(z,t)
>  %w = i*z;
>  w = i*[z(3)*z(2)'-z(2)*z(1)',-2*z(3)*z(1)'-z(1)*z(1),-3*z(2)*z(1)];
> end
>
> 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)];
> end
>
> % 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);
>
> % Write your own output/plot routines for the solution matrix Z, e.g.,
> %plot(Z(:,1),'r')
> plot(Z(:,1), 'r', Z(:,2), 'g', Z(:,3), 'b')

Thanks very much.  Your example was very helpful, for I think I  
managed to generalize it to compute more than 3 Fourier modes (is  
that the right terminology?) in solving Burgers' equation.  However,  
it's pretty slow.

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.

TIA.

---John.

%-------------begin script------------------
% Script for solving dz/dt = f(z,t), z(t0) = z0
% Here t0 is first coordinate of time node vector T
% at which solution is evaluated.
% User defines f as cplxfcn below as well as
% initial condition and time node vector.
% Vectors z, z0, f(z,t) must be row vectors.

%%z0 = 1;
%z0 = [1+i,1-i,-2]; % initial condition as row vector
%%T = linspace(0,2*pi,20);
%T = linspace(0,2,100); % time node vector

t0 = 0;   % initial time
tf = 10;   % final time
a = -1;   % left end point
b = 1;   % right end point
M = 100;   % 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;
K = 100;

% 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 = fft (zz);
z0 = fftzz (1:K);
T = linspace (t0, tf, M+1)';

function w = cplxfcn(z,t)
 %w = i*z;
 %w = i*[-z(3)*z(2)'-z(2)*z(1)',-2*z(3)*z(1)'-z(1)*z(1),-3*z(2)*z(1)];

global K;

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

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

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

  for j=1:K
    w(j) = -i*j/2*(U*V).'(j);
  endfor
end

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)];
end


% 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  uhat
uhat = zeros (K, M+1);
for j = 1:K
  uhat(j, :) = Z(:, j);
endfor

% define  u
u = zeros (N+1, M+1);
EXP = zeros (N+1, K);

for j = 1:K
  EXP(:, j) = exp (i.*j.*xx);
endfor

for j = 1:M+1
  u(:, j) = EXP*uhat(:, j);
endfor


% Write your own output/plot routines for the solution matrix Z, e.g.,
%plot(Z(:,1),'r')
%plot(Z(:,1), 'r', Z(:,2), 'g', Z(:,3), 'b')
plot(xx, real (u(:,1)), 'r', xx, real (u(:,2)), 'g', xx, real (u(:, 
3)), 'b')