# MathWorks Logo, Part Two. Finite Differences 2

Posted by **Cleve Moler**,

After reviewing the state of affairs fifty years ago, I use classic finite difference methods, followed by extrapolation, to find the first eigenvalue of the region underlying the MathWorks logo.

### Contents

#### Fifty Years Ago

My Ph. D. dissertation, submitted to the Stanford Math Department in 1965, was titled "Finite Difference Methods for Eigenvalues of Laplace's Operator". The thesis was not just about the L-shaped region, but the L was my primary example.

My thesis advisor was George Forsythe, who had investigated the L in some of his own research. At the time, many people were interested in methods for obtaining guaranteed upper and lower bounds for eigenvalues of differential operators. Forsythe had shown that certain kinds of finite difference methods could provide lower bounds for convex regions. But the L is not convex, so that made it a region to investigate.

Nobody knew the first eigenvalue very precisely. Forsythe and Wolfgang Wasow published a book in 1960, "Finite Difference Methods for Partial Differential Equations" (John Wiley & Sons), that features the L in several different sections. In section 24.3 Forsythe references some calculations he had done in 1954 that led him to estimate that

$$ \lambda_1 = 9.636 $$

He goes on to report an informal conjecture made in the late '50s from his colleague R. De Vogelaere that

$$ 9.63968 < \lambda_1 < 9.63972 $$

I repeat these numbers here to give you some idea of the accuracy we were working with back then. We will see in a later post that the right answer is just outside De Vogelaere's upper bound.

#### My Calculations in 1964

Here's a page from my thesis. The contour plot at the bottom of the page would eventually evolve into the MathWorks logo. The plot was made on a Calcomp plotter, which involved a computer controlled solenoid holding a ball-point or liquid ink pen moving across a paper rolled on a rotating drum. This was a very effective device and was used in computer centers worldwide for a number of years.

The calculations shown on this page are the first eigenvalue of the finite difference Laplacian for ten values of the mesh size $h$. The finest mesh, $h$ = 1/100, has 29,601 interior points. As I recall, the calculation of the first eigenvalue, which I did with a kind of inverse power method, took about half an hour on Stanford's campus main frame, an IBM 7090.

Here is a plot of these values I computed fifty years ago. I've added a black line at what we now know to be the limiting value, the first eigenvalue of the continuous differential operator. You can see that the finite difference values are converging quite slowly. Furthermore, they are converging from above. Forsythe was not going to get his lower bound for this nonconvex region.

thesis_graph

#### Slow Convergence

The show rate of convergence of the finite difference eigenvalue with decreasing mesh size is caused by the fact that the corresponding eigenfunction is singular at the reentrant corner in the region. The gradient becomes infinite at that corner. You can see the contour lines crowding together near the corner. If you tried to make a drum head or tambourine by stretching a membrane over an L-shaped frame, it would rip.

To be more precise, use polar coordinates, $r$ and $\theta$, centered at the reentrant corner. Note that to cover the interior of the region $\theta$ will run from $0$ to $\frac{3}{2}\pi$. It turns out that approaching the corner the first eigenfunction has the asymptotic behavior

$$ u(r,\theta) \approx r^\alpha \sin{\alpha \theta} $$

with $\alpha = \frac{2}{3}$. This implies that the $k$ -th derivatives of the eigenfunction have singularities that behave like

$$ \partial^k u \approx r^{\alpha-k} $$

When this fact is combined with the Taylor series analysis of the discretization error of the five-point finite difference operator we find that the first eigenvalue converges at the rate

$$ \lambda_{1,h} - \lambda_1 \approx O(h^{2 \alpha}) = O(h^\frac{4}{3}) $$

This is significantly slower than the $O(h^2)$ convergence rate obtained when the eigenfunction is not singular.

#### Difference Methods in MATLAB

Now let's use my laptop and the sparse capabilities in MATLAB. The functions `numgrid`, `delsq`, `spy`, and `eigs` have been included in the `sparfun` directory since its introduction. We begin by generating a small L-shaped numbering scheme.

```
n = 10;
L = numgrid('L',n+1)
```

L = 0 0 0 0 0 0 0 0 0 0 0 0 1 5 9 13 17 21 30 39 48 0 0 2 6 10 14 18 22 31 40 49 0 0 3 7 11 15 19 23 32 41 50 0 0 4 8 12 16 20 24 33 42 51 0 0 0 0 0 0 0 25 34 43 52 0 0 0 0 0 0 0 26 35 44 53 0 0 0 0 0 0 0 27 36 45 54 0 0 0 0 0 0 0 28 37 46 55 0 0 0 0 0 0 0 29 38 47 56 0 0 0 0 0 0 0 0 0 0 0 0

To illustrate how this scheme works, superimpose it on a finite difference grid.

Lgrid(L)

The statement

A = delsq(L);

constructs matrix representation of the five-point discrete Laplacian for the grid `L`.

```
whos A
issym = isequal(A,A')
nnzA = nnz(A)
```

Name Size Bytes Class Attributes A 56x56 3156 double sparse issym = 1 nnzA = 244

So for this grid, the discrete Laplacian `A` is a 56-by-56 sparse double precision symmetric matrix with 244 nonzero elements. That's an average of

ratio = nnz(A)/size(A,1)

ratio = 4.3571

nonzero elements per row. Each interior grid point is connected to its four nearest neighbors. For example, point number 44 is connected to points 35, 43, 45, and 53. So the nonzero elements in column 44 are

A(:,44)

ans = (35,1) -1 (43,1) -1 (44,1) 4 (45,1) -1 (53,1) -1

There are 4's on the diagonal and -1's is positions on the off-diagonals corresponding to connections between neighbors. They can be seen in the `spy` plot which shows the location of the nonzeros.

spy(A)

Scale `A` by the square of the mesh size and then use `eigs` to find its five smallest eigenvalues.

```
h = 2/n
A = A/h^2;
eigvals = eigs(A,5,'sm')
```

h = 0.2000 eigvals = 30.3140 27.7177 19.0983 14.6708 9.6786

#### My Calculations in 2014

Run for decreasing mesh size, up to the memory capacity of my laptop.

```
type L_finite_diffs
```

% L_finite_diffs % Script for MathWorks Logo, Part Two % Finite Differences opts.tol = 1.0e-12; eigvals = zeros(13,5); tic for n = 50:50:650 L = numgrid('L',n+1); h = 2/n; A = delsq(L)/h^2; e = eigs(A,5,'sm',opts); eigvals(n/50,:) = fliplr(e'); fprintf('%4.0f %10.6f %9.0f %10.0f %16.12f\n', ... n,h,size(A,1),nnz(A),e(5)) end toc save eigvals eigvals

```
type L_finite_diffs_results
```

n h size(A) nnz(A) lambda_1,h 50 0.040000 1776 8684 9.661331286788 100 0.020000 7301 36109 9.649547111146 150 0.013333 16576 82284 9.645736333124 200 0.010000 29601 147209 9.643932372382 250 0.008000 46376 230884 9.642903412412 300 0.006667 66901 333309 9.642247498956 350 0.005714 91176 454484 9.641797238078 400 0.005000 119201 594409 9.641471331746 450 0.004444 150976 753084 9.641225852818 500 0.004000 186501 930509 9.641035124947 550 0.003636 225776 1126684 9.640883203502 600 0.003333 268801 1341609 9.640759699387 650 0.003077 315576 1575284 9.640657572983 Elapsed time is 29.729023 seconds.

#### Extrapolation

I know the Taylor series analysis of the discretization error involves terms with both $h^{4/3}$ and $h^2$, so let's use these as the basis for extrapolation. Use backslash to do a least squares fit.

load eigvals eigvals = eigvals(:,1); n = (50:50:650)'; h = 2./n; e = ones(size(h)); X = [e h.^(4/3) h.^2]; c = X(4:end,:)\eigvals(4:end); lambda1 = c(1); n = (50:10:650)'; h = 2./n; fit = c(1) + c(2)*h.^(4/3) + c(3)*h.^2; plot(n(1:5:end),eigvals,'o',n,fit,'-',[0 700],[lambda1 lambda1],'k-') set(gca,'xtick',100:100:600,'xticklabel',... num2str(sprintf('2/%3.0f\n',100:100:600))) xlabel('h')

Frankly, I was surprised when the extrapolation of the finite difference eigenvalues produced this result.

```
format long
lambda1
```

lambda1 = 9.639723753239963

I know from techniques that I will describe in later posts that this agrees with the exact eigenvalue of the continuous differential operator to seven decimal places. Before I wrote this blog I never tried this computation on computers with enough memory to handle meshes this fine and was never able to obtain this kind of accuracy from extrapolated finite difference eigenvalues.

#### Eigenfunction

Let's do a contour plot of the eigenfunction. The tricky part is using the grid numbering in `L` to insert the eigenvector back onto the grid.

close n = 100; h = 2/n; L = numgrid('L',n+1); A = delsq(L)/h^2; [V,E] = eigs(A,5,'sm'); L(L>0) = -V(:,5); contourf(fliplr(L),12); axis square axis off

The plot shows off `parula`, the MATLAB default colormap with Handle Graphics II and Release 2014b. The colormap derives its name from the parula, which is a small warbler found throughout eastern North America that exhibits these colors. As our contour plot demonstrates, the color intensity varies smooth and uniformly as the contour levels increase. This is not true of `jet`, the old default. Steve Eddins has more about `parula` in his blog.

Get the MATLAB code

Published with MATLAB® R2014b

**Category:**- Differential Equations,
- Eigenvalues,
- History,
- Logo

### Note

Comments are closed.

## 2 CommentsOldest to Newest

**1**of 2

Love the eigenfunction contour plot, Cleve!

The parula colormap is named for the tropical parula. I’m planning to post more about the name later.

**2**of 2

Thanks for the clarification, Steve. I’m afraid that I don’t know much about birds.

## Recent Comments