MathWorks Logo, Part Four, Method of Particular Solutions Generates the Logo (2023)

The Method of Particular Solutions computes a highly accurate approximation to the eigenvalue we have been seeking, and guaranteed bounds on the accuracy. It also provides flexibility involving the boundary conditions that leads to the MathWorks logo.

Contents

Eidgenossische Technische Hochschule Zurich

In 1965-66, after finishing grad school, I spent a postdoctoral year at the ETH, the Swiss Federal Institute of Technology in Zurich. At the time, the faculty included three world renowned numerical analysts, Eduard Stiefel, Heinz Rutishauser, and Peter Henrici. Moreover, IBM had a branch of their research division in Ruschlikon, a Zurich suburb, where they held occasional seminars.In October, Leslie Fox, from Oxford, visited the IBM center and gave a talk about the work two of his students, Joan Walsh and John Reid, were doing on the L-shaped membrane. Expansions involving Bessel functions that matched the singularity of the eigenfunction near the reentrant corner were coupled to finite difference grids in the rest of the domain. Shortly after Fox's talk, it occurred to me to use the Bessel function expansions over the entire domain and avoid finite differences altogether.Within a matter of weeks I had done computations on the CDC 1604 mainframe at the ETH that had me very excited. To celebrate, I made an L-shaped ornament out of aluminum foil for the Christmas tree in our apartment.Henrici contributed some theoretical background relating what I was doing to work of Stephan Bergman and a Russian mathematician I. N. Vekua. Together with Fox, we published a paper in the SIAM Journal on Numerical Analysis in 1967.

Method of Particular Solutions

You don't need to know much about Bessel functions to appreciate the method of particular solutions. Using polar coordinates, a particular solution is any function of the form$$ v_\alpha(\lambda,r,\theta) = B_\alpha(\sqrt{\lambda} r) \sin{(\alpha \theta}) $$First of all, notice that the $\sin{(\alpha \theta)}$ portion insures that this solution satisfies the boundary condition$$ v_\alpha(\lambda,r,\theta) = 0 $$along the two rays with $\theta = 0$ and $\theta = {\pi}/{\alpha}$. For our L-shaped domain the reentrant corner has an angle of $\frac{3}{2}\pi$, so we will want to have$$ \alpha = \frac{2}{3}j $$where $j$ is an integer.The Bessel function $B_\alpha$ is defined by an ordinary differential equation that implies any $v_\alpha$ satisfies the partial differential equation$$ \Delta v_\alpha + \lambda v_\alpha = 0 $$If $\alpha$ is not an integer, then $B_\alpha(r)$ has another important property. Its derivatives blow up as $r$ goes to 0. This asymptotic singularity is just what we need to match the singularity in the eigenfunction.We can also take advantage of symmetries. The first eigenfunction is symmetric about the centerline of the L. We can also eliminate components that come from the eigenfunctions of the unit squares that are subdomains. This implies that for the first eigenfunction $\alpha = \frac{2}{3}j$ where $j$ is an odd integer that is not a multiple of 3, that is$$ \alpha_j = \frac{2}{3}j, \ \ j = 1,5,7,11,13, ... $$

Boundary Conditions

Form a linear combination of particular solutions.$$ u(\lambda,r,\theta) = \sum_\alpha c_\alpha v_\alpha(\lambda,r,\theta) $$This function satisfies our differential equation throughout the domain, satisfies the boundary condition on the two edges that meet at the reentrant corner, and has the desired singular behavior there. All we have to do is pick $\lambda$ and the coefficients $c_\alpha$ so that $u$ satisfies the boundary condition, $u = 0$, on the remaining edges of the L.The general idea is to pick a number of points on the outer boundary and investigate the matrix obtained by evaluating the particular solutions in the linear combination at these points. Vary $\lambda$ until that matrix is nearly rank deficient. The resulting null vector provides the desired coefficients.The details of how we carry out this general plan have varied over the years. The best approach is implemented in the function

, available here, described in the PDE chapter of Numerical Computing with MATLAB. It uses the matrix computation workhorse, the SVD, the singular value decomposition.The number of rows in the matrix involved is the number of points on the boundary and the number of columns is the number of terms in the sum. The figures below come from matrices that are 46-by-9. Let $(r_i,\theta_i)$ be the points chosen on the boundary. The matrix elements are$$ a_{i,j}(\lambda) = B_{\alpha_j}(\sqrt{\lambda} r_i) \sin{(\alpha_j \theta_i)} $$The function

is used to vary $\lambda$ and find a minimum of the SVD of this matrix. The resulting matrix is essentially rank deficient and the resulting right singular vector is a null vector that provides the coefficients.The minimum is found at$$ \lambda = 9.639723844 $$Here are the first four particular solutions. The vertical scale reflects the magnitude of the corresponding coefficient.

 mps_figs

All Four Terms

None of these four terms is remotely close to zero on the outer boundary, but when they are added together the sum is zero to graphical accuracy.

 mps_figs(4)

So this is actually the first eigenfunction of the L-shaped membrane.

Just Two Terms

Now we can have some fun. Here is a

plot of the sum of just the first two particular solutions.

 mps_figs(2)

This function is an exact solution of the partial differential equation$$ \Delta u + \lambda u = 0 $$with the value of $\lambda$ reported above. It satisfies the boundary condition $u = 0$ on the two edges meeting at the reentrant corner, but it fails significantly on the outer edges. So it is only with some artistic license that we can refer to it as the first eigenfunction of L-shaped membrane.When I first saw this

plot, it didn't have this gorgeous color map, but I found the shape so attractive that I kept it around. The evolution of the shading, lighting and colors that turned this into our logo is the subject of my next blog.

Upper and lower bounds

The Fox, Henrici and Moler paper was not just about the method of particular solutions. We were also concerned about finding guaranteed upper and lower bounds for the exact eigenvalue of the continuous differential operator. A theorem proved in the paper relates the discrepancy in meeting the boundary condition to the width of a small interval containing the eigenvalue. In 1966 we proved$$ 9.6397238_{05} < \lambda < 9.6397238^{84} $$Today with

on MATLAB, I am finding a value of $\lambda$ in the center of this interval.

 [L,lambda] = membranetx(1); format long lambda

lambda = 9.639723844021997

Reviving the Method of Particular Solutions

In 2005 Timo Betcke and Nick Trefethen published a paper in SIAM Review that pointed out a serious shortcoming in the approach we described in the Fox, Henrici and Moler paper. When the method is attempted on other domains, especially domains with more than one singular corner, the matrix has a tendency to appear nearly rank deficient for all values of $\lambda$. It is difficult to identify local minima. The recommended cure is to add a number of interior points and require the approximate eigenfunction to be nonzero there.Betcke and Trefethen also prove an approximation theoretic result that, when combined with their computations, adds four more digits to our rigorous bounds on the first eigenvalue.$$ 9.63972384402_{16} < \lambda < 9.63972384402^{22} $$Moreover, they are so confident in their numerical results that they are sure the exact eigenvalue lies in the center of this interval, and they are pretty sure they even know the next digit.$$ \lambda \approx 9.63972384402195 $$This agrees with what I get out of

to almost full double precision accuracy. I think we finally have this fellow pinned down.I recommend that anyone reading this blog who wants to know more about the method of particular solutions and computing the eigenvalues of the L-shaped membrane read the Betcke and Trefethen paper.

References

Leslie Fox, Peter Henrici, and Cleve Moler, Approximations and Bounds for Eigenvalues of Elliptic Operators, SIAM Journal Numerical Analysis 4, 89-102, 1967, <https://blogs.mathworks.com/images/cleve/Fox_Henrici_Moler.pdf>Timo Betcke and Lloyd N. Trefethen, Reviving the Method of Particular Solutions SIAM Review 47, 469-491, 2005, https://people.maths.ox.ac.uk/trefethen/publication/PDF/2005_112.pdfCleve Moler, Numerical Computing with MATLAB, Chapter 11, Partial Differential Equations, 2004. <https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/pdes.pdf>

Published with MATLAB® R2014b