Jonathan Ong's Blog

This is MATLAB code for the Level Set and Fast Marching Methods since I can’t find any online.  This is for Sethian‘s Math 228B course at UC Berkeley.

Inputs:

• H => mesh spacing = 2^-H
• K => time step = 2^-K
• OBS => if OBS == 1, obstacles are turned on
• PLOT => if PLOT is defined, the contour is plotted at every time step

fastmarching.m

levelsetmethod.m

The following is sample MATLAB code for solutions to hyperbolic equations for UC Berkeley course Math 228B – Numerical Solutions of Partial Differential Equations. These are codes for the following schemes: Lax-Friedrich, Lax-Wendroff with a viscosity term, Godunov, Glimm, and the exact solution to a simple Reimann problem.  Note that Lax-Wendroff doesn’t work very well without a viscosity term; it may be due to my own coding error.

These samples solve Burger’s equation u_t + uu_x = 0 with initial condition either u = 0 for x < 0, u = 1 for x >= 0 or u = 1 for x < 0, u = 0 for x >= 0. Glimm’s, Godunov’s, and Reimann’s  assume that u >= 0, otherwise solving the exact Reimann solution would be a little more tricky.

For all the following codes, these inputs are the same:
-  dx = 2^-H
-  dt = 2^-K
-  POS = 1 => u = 0 for x < 0 & u = 1 for x >= 0
-  POS != 1 => u = 1 for x < 0 & u = 0 for x >= 0
-  PLOT = 1 => plots the solutions from time 0 to 1

Functions:
-  Lax-Friedrich: [x,u] = LF(H,K,POS,PLOT)
-  Lax-Wendroff: [x,u] = LW(H,K,POS,e,PLOT)
- e is viscosity term for u_xx, doesn’t work very well without it
-  Glimm: [x,u] = Glimm(H,K,POS,PLOT)
-  Godunov: [x,u] = Godunov(H,K,POS,PLOT)
-  Reimann: y = Reimann(x,ul,ur)
- Solves the exact solution for a point between two states, ul = left state, ur = right state, x = location with 0 being the center. Assumes that ul & ur >= 0, otherwise the solution would be incorrect.

I had a phone interview this morning for a Strategist position at Goldman Sachs which I regrettably failed. As you may have known, my car’s engine blew, causing the ride back from Berkeley to take an extra three hours and I ended up not sleeping before the phone interview to take care of my car ASAP. Thus, I was incapable of any profound thought process.

The interviewer asked me these three technical questions, all of which I did not solve until after I slept. These are actually good questions, much better than questions given to me in interviews pertaining to similar positions. Let’s see if you can solve them.

1. Which is larger, the expectation of the exponential of x or the exponential of the expectation of x? In other terms, which one is larger: E[e^x] or e^(E[x])?
2. What is the expectation of the exponential of a normally distributed random variable? Find E(e^x) where x ~ N(0,1).
3. What is the expected number of flips until you flip two heads in a row?

Let me know if you would like the answers.

Suppose you have  a nonuniform (or uniform) domain x and range y for which you’d like to approximate the kth derivative at xbar. This program calculates this approximation using Taylor expansions and solving a Vandermonde matrix (or the transpose of one) to calculate constants c. You then dot-product c with y to obtain the approximate kth derivative at xbar (which depends if your input is a row or column vector, so I did not include that part of the program). To understand the program, read these lecture notes from UC Berkeley’s Numerical Solution to Differential Equations course.

Notes:

• The dimension or length(x) has to be greater than k otherwise the method will fail. This make sense since you can’t approximate the first derivative using only one point.
• x and xbar are standardized normal to avoid a high condition number when solving the Vandermonde matrix.
• The estimated error is not included, though it is approximately to the order of length(x)-k assuming no symmetry.
• Remove the last 5 lines from the code if you don’t care for the iterative refinement (will make the constants accurate to almost machine precision).

Code (copy and paste):

function c=mkfdstencil(x,xbar,k)
% Jonathan Ong - http://jonon.gs/blog/?p=411
% Calculates the finite difference approximation weights given a nonuniform
% (or uniform) domain x, point of interest xbar, and kth derivative.
% The inputs are (x,xbar,k) and the output is c.
% c is the constants s.t. c.*y (or y*c if y is a row vector) ~ f^k(xbar)
%
% Example:
% x = [1:5]; % Domain of integers 1 to 5
% y = exp(x); % Range of the function exp(x)
% c = mkfdstencil(x,3,2); % Find the finite difference weights for the kth derivative using all 5 points
% % c =
% %   -0.083333333333333
% %    1.333333333333333
% %   -2.500000000000000
% %    1.333333333333334
% %   -0.083333333333334
% KthDerivAt3 = y*c; % The approximate solution = 19.841479123877672
% ActualSolution = exp(3) % The actual solution = 20.085536923187668
%
% Notes:
% Calculated by solving a system of linear equations of taylor expansions
% centered at xbar and calculate at the nodal points of x.
% x and xbar are standardized normal to avoid a high condition number when
% solving the Vandermonde matrix
% Estimated error is not included but is on the order of length(x)-k without
% symmetry factored.
% Iterative refinement is included so the constants are accurate to almost
% machine precision.

if round(k) ~= k
    error('k must be an integer');
end
LX = length(x);
if LX < k+1
    error('not enough sample points in x to approximate the kth derivative');
end
A = ones(LX, LX);
b = zeros(LX, 1);
b(k+1) = 1;
AVG = mean(x);
STD = std(x);
x = (x-AVG)/STD;
xbar = (xbar-AVG)/STD;
x2 = x-xbar;
for i = 1:LX-1
    A(i+1,:) = (x2).^i/factorial(i);
end
c = (A\b)/STD^k;
for i = 1:3
    r = b - A*c*STD^k;
    r2 = (A\r)/STD^k;
    c = c + r2;
end

Today in my graduate Numerical Linear Algebra course, Professor Kahan, who substituted for Professor Demmel, discussed the differences between the nonsymmetrical and symmetrical eigenvalue problem. He then asked, “What do matrices and potato chips have in common?”

Attempting to be an active listener, I thought about the possibilities. Since we were discussing symmetric matrices, I thought about a saddle point which has a shape of a Pringles potato chip if it is within a certain boundary. This of course a logical conjecture since the Hessian of a matrix, which tells you if a critical point is a saddle point, is a symmetric matrix.


I then thought about particular shapes which could be present in topology or geometry that I’m not aware about since I prefer numerics, but I wasn’t sure how that could relate to a matrix. I only thought of this because the words “tensor” and “torus” popped into my mind when I heard “potato chips”, the former having to do with matrices (among other things) and the latter with shapes.

Alas, after a good 15 seconds of thinking, he told us the answer:

You can’t just have one.

…ARE YOU SERIOUS? THIS WAS A JOKE? I couldn’t tell because even after he gave his punch line, he still wasn’t laughing, but boy was I laughing inside. I couldn’t believe that he was joking and that I took him so seriously. However, he was talking about dual spaces and how all matrices have one. Too bad I fell asleep in my undergrad linear algebra course so I don’t remember any of that.

The whole lecture was made of these subtle jokes, though this was the only explicit one. I remember his saying, “But it’s miraculous that I remember any of this at my age.” He’s 76 but he sure doesn’t act like it. This was probably one of the single greatest math lecture’s I’ve ever had.

If you don’t know who Kahan is, he’s a very notable mathematician and computer scientist; he even has his own nontrivial wiki! He was the father of the IEEE standard and won a Turing award. Around campus he’s known as the man who wears overalls and despises tech companies for not supporting processors with higher than 64-bit architectures for its numerical accuracy. Who knew he was such a character?