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.

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

Given a Butcher array (A, b, c) of a Runge-Kutta method, one can deduce it’s order using graph theory. For a RK method to be order s, it must satisfy all its graph theory conditions for orders 1 to s. One could find these order conditions from a graduate-level differential equations textbook such as A First Course in Numerical Analysis of Differential Equations by Arieh Iserles (Cambridge Texts in Applied Mathematics, 2009).

You can find the MATLAB code here.

I had trouble finding information on the recursive LU decomposition algorithm for Gaussian Elimination and even harder time finding MatLab code for my graduate course in Matrix Computations at Berkeley. After reading various documents I finally figured out the code and decided to show it share it with you all.

The purpose of a recursive algorithm (a.k.a. “divide and conquer”) is to minimize data movement between fast and slow memory without actually specifying the size of the cache.

For more information, you can view my Professor’s notes. He sites a source written by Toledo in 1997 which I am not aware of. You can also view notes I found online written by Dongarra, Eijkhout, and Luszczek of the University of Tennessee titled Recursive Approach in Sparse Matrix LU Factorization.

For my code, we assume that A is an n by m matrix with n > m and m = 2^k for some integer k. My code is not optimized and is meant for educational purposes. Passing through pointers of A instead of the submatrices and not creating new variables to name parts of the matrix would speed up the program. The output P is not the permutation matrix, but you can create permutation matrix using two lines of code.

The results of this code is that the solution are accurate within the error bounds specified in page 49 of Demmel’s Applied Numerical Linear Algebra. However, a significant increase in error occurs in the existence of columns of equal values.

.M file