Trying to solve u_tt - u_xx = 0

[John B. Thoo]

Hello, everyone.

I'm trying to solve the wave equation

(E)  u_tt - u_xx = 0,

(I)  u(x,0) = u0,  u_t(x,0) = u0_t.

Since lsode can solve a system of first-order, I define  v = u, w =  
u_t  and write (E) as

(E-sys)  v_t = w,  w_t = v_xx,

(I-sys)  v0(x) = u0,  w0(x) = u0_t.

Now comes my difficulty: putting all this into Octave.  Actually, I  
think my code (given below) works correctly.

My question is this: Because I'm new to Octave and numerics in  
general, can you tell me if my code is decent or how I can improve on  
it?

%------------------begin here--------------------
%% Set time parameters
%%
t0 = 0;   % initial time
tf = 2;   % final time
M = 20;   % number of time intervals (at least 10 for animation)
dt = (tf - t0)/M;   % time mesh width for plotting (see next line)
t = t0:dt:tf;   % time interval for plotting; does not affect the  
time mesh
                % used in the solver (either "lsode" or "ode45")  
which uses
                % its own adaptive time stepping

%% Set space parameters
%%
ord = 8;   % n = 2^ord  space grid points --> n  Fourier modes(?)
x = 0:(2^ord-1);   % }
x = x/2^ord;       % } space interval
x = 2*pi*x;        % }
dx = x(2)-x(1);   % space mesh width

global lx;
lx = length (x);


%% Define the initial conditions

%u0 = exp (-25*(x - 1).^2);
%u0 = exp (-25*(x-1).^2) + 0.6*exp (-9*(x-3).^2);
u0 = cos (x);

%u0_t = zeros (1, lx);
u0_t = -sin (x);


%% "make_k" is to multiply  k  by the fft  correctly
function k = make_k (n)
  k = 1:n;
  k = heaviside(n/2 + 1 - k).*(-1*(k - 1)) + heaviside(k - n/2).*(n  
- k + 1);
  k = -k;
endfunction

global k;
k = make_k (lx);

%% "dderiv" calculates the second derivative in Fourier space
function d2ydx2 = dderiv (u, k)
  d2ydx2 = real (ifft (-k.*k.*fft (u).'));
endfunction

%% Define  v_t  and  w_t
function v_t = f (w, t)
  v_t = w;
endfunction

function w_t = g (v, t)
global k;

  w_t = dderiv (v, k);
endfunction

%% Define the system of ODEs to be solved
function U_t = F (U, t)  % use this for "lsode"
global k;
global lx;

  v = U(1:lx);
  w = U(lx+1:2*lx);

  v_t = f (w, t);
  w_t = g (v, t);

  U_t = [v_t.', w_t];
endfunction


%% Solve the system of ODEs

%lsode_options ("absolute tolerance", 1.0e-2*ones (1, lx));
%lsode_options ("relative tolerance", 1.0e-2);

u = lsode ("F", [u0, u0_t], t);



%% Plot the solution
clf

u = u.';
u = u(1:lx,:);

%% 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

N = 10;   % number of time intervals to plot
time = ones (1, N);
time(2:N+1) = floor (M/N)*(1:N) + 1;

for j = 1:N+1
  plot (x, u(:,time(j)));
  title (strcat ('time =', num2str ((tf-t0)*(time(j) - 1)/M)));
%  axis ([0 2*pi -1.5 1.5]);
  axis([0 2*pi 2*min(real(u0)) 2*max(real(u0))]);
  pause(wait);
endfor
%-------------------end here---------------------


Thanks for your patience and your help.

---John.