One of my goals when I started my website just after I finished grad school was to start sharing little snippets of code I had used or developed over the course of my undergraduate and graduate degree. One thing I did use, and still use frequently, is the Fourier transform. The Fourier transform is used in countlessly many algorithms in seismic imaging because of it’s relationship to sampling theory, ability to identify frequency content in signals, and provide solutions to differential equations.

One day I stumbled upon a neat series of lectures by my undergraduate NSERC supervisor, Dr. Michael Lamoureux, for a Seismic Imaging Summer School in 2006. Hidden in a footnote of this set of lecture notes contains one of my favorite quotes about a mathematical theorem, namely the divergence theorem:

And, as I like to explain to my kids, the Divergence Theorem simply says you can measure how much methane gas a cow produces by either measuring how much is produced in each cubic centimeter of the cow, or simply by measuring how much leaks out its mouth and rear end, and other places on the surface of the cow. – Michael Lamoureux

However, I digress! In these lecture notes Michael describes a variety of theoretical considerations when investigating partial differential equations, and in particular the famous wave equation with initial conditions:

In the third lecture of these notes, the Fourier transform is used to find a solution of the equation. This is also given as an example in Elias Stein’s famous book on harmonic analysis (though I have it written here as Page 395, I’m not sure which volume. It’s likely volume 1). In this case, the Fourier transform of the solution can be written as the sum of sines and cosines:

In the future, I might include the derivation though for now this will be sufficient. We simply apply the inverse Fourier transform to obtain a solution. In 2009 or so, I wrote a little Matlab code to demonstrate this solution. There are some scaling factors to match that of the Fourier transform implementation of Matlab. Luckily I’m a bit of a hoarder when it comes to old school work so here it is:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 |
% Solution of wave equation with initial % conditions and constant velocity, u(x,0) = f0(x);, (d/dt)f(x,0) = f1(x); % (d^2u/dx^2) = (d^2u/dt^2) function [u] = waveqnsteinsol(f0,f1,t,dx) N = length(f0); ft_f0=fft2(f0); ft_f1=fft2(f1); % Both the positive and negative frequencies corresponding to the frequency domain. W=[0:floor((N-1)/2) -ceil((N-1)/2):-1]/(dx*N); [w,z] = meshgrid(W,W); % Compute the amplitude functions amp1 = cos(2*pi*(sqrt(w.^2+z.^2))*t); amp2 = 2*pi*sinc(2*sqrt(w.^2+z.^2)*t); % Finish the operations inside the integrands t1ft = ft_f0.*amp1; t2ft = ft_f1.*amp2; % Inverse Fourier transform t1 = ifft2(t1ft); % t2 = ifft2(t2ft); % % Sum of two FIO's u = t1 + (1/(2*pi))*t2; |

To be honest, there are a lot of poor programming practices going on here. For example when I wrote this, prior to my introduction into software development, I probably gave myself a nice little pat on the back for putting in those comments. Those who employ clean coding practices see where I’m going with this. If you take a look through the code, every time there is a comment it could be naturally refactored into a function with a descriptive name. Moreover, who am I trying to kid with variable names like “W”! No one would quickly realize that I’m referring to a grid of frequencies. Though it is academic code in nature, it’s still embarrassing!

Moving forward, we can then use this function to plot the solution for a couple different initial conditions:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 |
%% 2-d Wave equation % Sample frequency Fs = 10; % Corresponding sample rate dx = 1/Fs; % Domain of interest x = -5:dx:5; [X,Y] = meshgrid(x,x); N = length(x); %% Gaussian Initial condition f1 = zeros(size(X)); f0 = 1*exp(-(X.^2+Y.^2)); %% Square box initial condition f0 = ones(size(X)); f1 = zeros(size(X)); f0(1:40,:) = 0; f0(end-40:end,:) = 0; f0(:,1:40) = 0; f0(:,end-40:end) = 0; %% View the initial condition imagesc(f0) %% Plot of the solution for k seconds for k=0:.1:10 u = waveqnsteinsol(f0,f1,k,dx); figure(1) %imagesc(x,x,real(u)), axis image, colormap hot, colorbar; mesh(X,Y,real(u)); zlim([-1 1]) title(['T = ',num2str(k) ]) end |

Here are two sample videos:

Notice that we see some funky stuff happening at the boundaries. The mathematical solution given is valid on the entire plane, so it should just disappear off our screen. The fact that there is reflections at the boundary is an artifact of the discrete nature of computational mathematics, and these are typically handled using things like absorbing boundary conditions or perfectly matched layers. The easiest, but most computational inefficient way of handling it, would be to make an even larger grid and only record for the amount of time it takes for the reflection to reach the area you are recording.