Cleaning Up Mac Ports

This is a good tip, as stuff tends to get installed for one reason or another, and then never removed.  I also don't entirely understand why macports doesn't remove the inactive ports when upgrades are installed.

Checking Progress on Asynchronous Jobs

I recently discovered the wc command in *nix environments, which is a handy solution for when I run asynchronously processed jobs.  I often need to run ensembles of simulations.  If I process them serially, and put an output command somewhere in the script, I can immediately see how far along I am by checking a log file.  But if the data is processed asynchronously, it can be a bit tougher to see what's going on.

Running the following on my log file tells me the number of lines so far:

wc -l log.out

Hence, if the log file has n lines, and I have a total of N cases, then, assuming each case produces a single output statement, n/N % of my job is done.

Graphics Grievances

When refereeing papers with plots of numerical simulations, I very often have to write back to the authors telling them that figures, while readable on the screen, where they can be blown up, are completely illegible when printed out in black and white.

When preparing images (MATLAB and matplotlib images anyways), unless you intend the figure to take up all of a letter sized paper, fonts and line widths must be checked for readability.  For the typical size figure that you'll see in a publication, a font size of at least 16, and often 18, is needed for the axes labels, tick marks, and other elements.  Line widths should generally be at least 2.  These things can easily be set in MATLAB using the exportfig script.  For matplotlib, consider editing your ~/.matplotlib/matplotlibrc file. You can set the default values for all of these graphical items as desired.

Pseudo-Spectral Methods

I was previously affiliated with the University of Toronto, where I wrote some pages for the math department's wiki on how to get pseudo-spectral algorithms up and running for nonlinear wave equations such as Nonlinear Schrödinger (NLS) and Korteweg - de Vries (KdV).  Much of this was based on Trefethen's Spectral Methods in MATLAB.

While that's still an excellent book that I recommend to people who are starting out, it has a very serious flaw in its discussion of nonlinear problems.  Trefethen's p27.m example, which simulates KdV solitons does not include a de-aliasing step.  Early in my career, this led me to think that de-aliasing was unnecessary.  That was, until, I started looking at problems with cubic nonlinearities; the p27.m example is only quadratic.

Consider focusing, cubic NLS on the torus, \(\mathbb{T}\),
\[
i u_t + u_{xx} + |u|^2 u = 0, \quad u: \mathbb{T} \times \mathbb{R} \to \mathbb{C}.
\]
When posed on the whole real line, instead of the torus, it admits soliton solutions,
\[
\Phi(x,t) = e^{i \lambda^2 t} \sqrt{2} \lambda {\rm sech}(\lambda x), \quad \lambda \neq 0.
\]
Since they are exponentially localized, if we take \( \lambda \) sufficiently large, we can approximately simulate it on the finite domain size.

Modifying p27.m for NLS, we have:

% Set discretization parameters
N = 512;
dt = .1/N;
x = (2*pi/N)*(-N/2:N/2-1)';
% Set initial condition
lam = 10;
u = sqrt(2) * lam * sech(lam * x);
v = fft(u);

% Construct integration operators
k = [0:N/2-1 0 -N/2+1:-1]';
E = exp(-.5i * dt * k.^2);
E2 = E.^2;
g = 1i * dt;

% Set simulation parameters
tmax = 5; nplt = floor((tmax/50)/dt); nmax = round(tmax/dt);
udata = u; tdata = 0;

% Integrate using the Integrating Factor method
for n = 1:nmax
  t = n*dt;
  a = g.*fft(abs( ifft(     v   ) ).^2 .* ifft(v));
  b = g.*fft(abs( ifft(E.*(v+a/2)) ).^2 .* ifft(E.*(v+a/2)) );
  c = g.*fft(abs( ifft(E.*v + b/2) ).^2 .* ifft(E.*v + b/2) );
  d = g.*fft(abs( ifft(E2.*v+E.*c) ).^2 .* ifft(E2.*v+E.*c));
  v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
  if mod(n,nplt) == 0
    u = ifft(v);
    udata = [udata u]; tdata = [tdata t];
  end
end

% Plot and save
figure;
for n = 1:length(tdata)
  plot(x, abs(udata(:,n)),...
    x, real(udata(:,n)), x, imag(udata(:,n)),...
    x,abs(udata(:,1)), 'linewidth',2);
  xlim([-pi, pi]);
  ylim([-16,16]);

  set(gca,'fontsize',16);
  xlabel('x');
  legend('|u|', 'Re. u', 'Im. u', '|u_0|');
  title(sprintf(' t = %g', tdata(n)));
  drawnow;
  saveas(gcf,sprintf('frame%i.png', n));
end

This generates the following output:

Notice the appearance of very high frequency oscillations.  These are aliasing errors. De-aliasing them, by zeroing out the upper half of the frequencies in Fourier, produces the following, desired, output:

This is the result of the following code.  Notice that I doubled the number of grid points to compensate for the de-aliasing.

% Set discretization parameters
N = 1024;
dt = .1/N;
x = (2*pi/N)*(-N/2:N/2-1)';
% Set initial condition
lam = 10;
u = sqrt(2) * lam * sech(lam * x);
v = fft(u);
v(N/4+2:end - N/4+1) = 0;

% Construct integration operators
k = [0:N/2-1 0 -N/2+1:-1]';
E = exp(-.5i * dt * k.^2);
E2 = E.^2;
g = 1i * dt;

% Set simulation parameters
tmax = 5; nplt = floor((tmax/50)/dt); nmax = round(tmax/dt);
udata = u; tdata = 0;

% Integrate using the Integrating Factor method
for n = 1:nmax
  t = n*dt;
  a = g.*fft(abs( ifft(     v   ) ).^2 .* ifft(v));
  a(N/4+2:end - N/4+1) = 0;
  b = g.*fft(abs( ifft(E.*(v+a/2)) ).^2 .* ifft(E.*(v+a/2)) );
  b(N/4+2:end - N/4+1) = 0;
  c = g.*fft(abs( ifft(E.*v + b/2) ).^2 .* ifft(E.*v + b/2) );
  c(N/4+2:end - N/4+1) = 0;
  d = g.*fft(abs( ifft(E2.*v+E.*c) ).^2 .* ifft(E2.*v+E.*c));
  d(N/4+2:end - N/4+1) = 0;
  v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
  if mod(n,nplt) == 0
    u = ifft(v);
    udata = [udata u]; tdata = [tdata t];
  end
end

% Plot and save
figure;
for n = 1:length(tdata)
  plot(x, abs(udata(:,n)),...
    x, real(udata(:,n)), x, imag(udata(:,n))...
    , 'linewidth',2);
  xlim([-pi, pi]);
  ylim([-16,16]);

  set(gca,'fontsize',16);
  xlabel('x');
  legend('|u|', 'Re. u', 'Im. u');
  title(sprintf(' t = %g', tdata(n)));
  drawnow;
  saveas(gcf,sprintf('frame%i.png', n));
end

While I still recommend Trefethen's book to people, I would supplement it with another text that addresses the aliasing problem. Two books I like are Boyd's, Chebyshev and Fourier Spectral Methods, and Canuto et al's, Spectral Methods: Fundamentals in Single Domains.

(Partial) Apologia

While I had hoped to make this blog free from any particular programming ideology (only FOSS, or only GNU/Linux software, or only Apple, or only Satan's licensed software), I have rapidly realized that there is one platform which I will never write about, Windows. I haven't had a Windows machine in almost a decade. My personal machines have been OS X, and the machines at the institutions I've been at have generally, but not always, been Linux. This has led me to ask many a summer research student to do one of the following:

1. Install a Linux distribution on their own machine;
2. Use the department's Linux machines;
3. Figure everything out in Windows on his own;
4. Find a new summer research advisor.

I have felt somewhat guilty about this, but I'm just not knowledgeable enough to know how to set up different software libraries within a Windows programming environment, how to connect that machine to a remote server, or much else beyond web browsing and email.

Despite this, I recognize that Windows, as desktop environment, has come a long way in ten years, and it is no longer synonymous with the blue screen of death.

nohup and MATLAB

Occasionally, I need to run a large ensemble of MATLAB simulations and would prefer to offload this to a server. This can easily be accomplished with the following script, called matbg.sh (MATLAB background shell script):

#! /bin/bash
# Clear the DISPLAY.
unset DISPLAY  # unset DISPLAY for some shells
# Call MATLAB with the appropriate input and output,
# make it immune to hangups and quits using ''nohup'',
# and run it in the background.
nohup nice -10 matlab -nodisplay -nodesktop -nojvm -nosplash <$1 >$2 &

Log on to your machine, place this script in the same location as the m-file you wish to run and then execute:

./matbg.sh myscript1.m mylogfile.out

You're then welcome to log out and wait for the job to finish, but make sure that your m-file saves your output for post processing. If for some reason, such as if you are using some aspect of the parallel computing toolbox, you need Java, just remove the -nojvm flag from the above script. The next you may want is to be alerted to the completion of the job. This can be accomplished by tacking on the following to the end of your m-file:

% Email message.  Feel free to add in details about how long the job took,
% if it completed successfully, etc.:
mesg = sprintf('finished job');
% The command which will send the email:
cmd =...
sprintf(['echo "sending email notification" | mail -s "%s" ''%s'],...
mesg,'my.email@domain.com');
system(cmd);

Making Animations from Image Sequences

Despite making animations out of image sequences all the time, I often forget the right command line arguments, so I'm going to post them here for posterity.  Given a sequence of images in an acceptable format, say PNG, labeled frame001.png, frame002.png, etc., we can create an animation using FFMPEG with the following command:

ffmpeg -r 10 -i frame%d.png movie.mp4

which will generate a .mp4 formatted animation with 10 fps.  I use ffmpeg on the command line using the one I got from MacPorts.

I should note that I was a big fan of Quicktime Pro, which would nicely import image sequences and spit out animations in all sorts of formats.  It was an affordable piece of software, about $30 as I recall, but Apple seems to have abandoned it in favor pushing people to buy Final Cut, which is out of my budget just for making animations.

Something else that people often want to do is include animations in their Beamer presentations.  I've had success following these instructions with .mp4 formatted animations.  However, I have only gotten the actual presentation to work properly with Acrobat; Skim never seems to work for me.

Greetings

Welcome to Numerical Flux, a blog on scientific computing, numerical analysis, applied math and related issues.  The intention of this blog is twofold.  First, I aim to offer information and examples of all the work that goes into creating journal publications (codes, scripts and assorted tricks), but is never really part of the publication.  Second, I intend to survey a variety of topics and tools without adhering to any particular ideology.  I hope to have posts on MATLAB, SciPy, Fortran, MPI, etc.