# Plot MATLAB icon

I started to explore this on a whim and haven't succeeded yet… Some introduction for the icon is found here and here, but I can't understand it very well. (I admit that, though playing with NDSolve a lot, I suffered from my unstable foundation of knowledge for PDE.)

A MATLAB code sample for the icon is found at the end of the document here and I have no doubt that it can be translated to Mathematica code mechanically, but it'll be so boring! Can we plot it in a Mathematica style? Has anyone tried it before?

• Take a look at the first "simple example" here. ;-) Unfortunately Mathematica does not have direct support for elliptic differential equations. It doesn't mean you can't do it, but you have to work a bit harder. Take a look here and here. Commented Feb 27, 2014 at 4:48
• I'm not going to do it now because it's midnight here, but what you need to do: 1. discretize the L shaped region using a square grid, and index each cell of the grid with a single index 2. construct the Laplacian as a SparseMatrix lap acting on a vector containing the cells, in order of their indices 3. vec = Eigenvectors[lap, -1] 4. Map this eigenvector back onto the grid and plot it. Commented Feb 27, 2014 at 4:56
• My favourite reference in this site on solving elliptic equations is here and very informative IMO
– gpap
Commented Feb 27, 2014 at 10:32
• For more info on the MATLAB logo by Cleve Moler, see these five blog entries. Commented Sep 3, 2017 at 5:58

While the other answers are nice, the icon deserves a closer look:

Note, in particular, that four of the six edges are not constrained by the ostensible Dirichlet boundary conditions, nor is it clear that they solve a Neumann problem. And indeed, as I noted in the comments this is supported by the OP's first link.

In short, to produce the logo, they took 'pac-man' eigenfunctions around the convex corner, which is the trickiest bit to describe, and used a least-squares approach to impose the boundary conditions. Finally, they discarded all but two terms in the expansion:

After being so careful to satisfy the boundary conditions, the logo uses only the first two terms in the sum. This artistic license gives the edge of the logo a more interesting, curved shape.

In my view, this essentially says it is pointless to try and impose the boundary conditions at all. Since they took a certain artistic license, in replicating it one can simply try and do it by hand. We know that the plot is a function of the form $$v(r,\theta)=\sum_{j=1,2} a_j J_\frac{2}{3}(k_j r)\sin(\tfrac{2}{3}\theta),$$ which has only four free parameters. These are enough that they can, to a good precision, simply be set by hand. I therefore used the Manipulate construct

Manipulate[
Plot3D[(
a BesselJ[2/3, λ Sqrt[x^2 + y^2]] Sin[
2/3 (Arg[y - I x] + π)] +
b BesselJ[2/3, μ Sqrt[x^2 + y^2]] Sin[
2/3 (Arg[y - I x] + π)]
) Boole[x >= 0 || y >= 0], {x, -1, 1}, {y, -1, 1},
BoxRatios -> Automatic, Mesh -> None, Exclusions -> None, PlotPoints -> 50, Axes -> False
]
, {{a, 1.275}, 0, 3}, {{b, 0.805}, 0, 3}, {{λ, 3.18}, 0, 10}, {{μ, 1.96}, 0, 10}]


to find good guesses for the parameters, which are the defaults above. The result is, I think, pretty close to the original:

One could also, if so inclined, attempt to implement their least-squares procedure. This sounds like a hazier problem to me: it is no longer "find an accurate enough solution of this specific problem", but rather "find an accurate solution that will also look like ours when you make this arbitrary restriction". It may be that many of the possible other options for the choices they make ($m$ points on the boundary, $n$ fundamental solutions) will yield similar-looking plots. But that's much too much work, I think.

• Agreed on the boundary constraints. For what it's worth, my image (done quite some time ago) was based on an earlier icon which (as I recall) was constrained at the edges. Commented Mar 4, 2014 at 10:42
• I must admit I ignored the stuff about how they actually created the logo, finding the finite difference solution more interesting! Commented Mar 4, 2014 at 11:49

Here's my attempt. To get the matrix representing the Laplacian I use LaplacianFilter on an array of symbols and CoefficientArrays to extract the coefficients.

n = 200;

shape = ArrayPad[ConstantArray[0, {n/2, n/2}], {{0, n/2}, {0, n/2}}, 1];
shapeVector = Flatten @ Position[Flatten @ shape, 1];

symbolArray = Array[x, {n, n}];
symbolLaplacian = LaplacianFilter[1.0 symbolArray , 1, Padding -> 0];
lapMatrix = Last@CoefficientArrays[
Flatten[symbolLaplacian][[shapeVector]],
Flatten[symbolArray][[shapeVector]]];

{ev} = Eigenvectors[lapMatrix, {-1}];

result = Flatten[shape];
result[[shapeVector]] = ev;
result = Partition[result, n];

ListContourPlot[result]
ListPlot3D[result, BoxRatios -> {1, 1, 0.8}, Mesh -> False,
PlotStyle -> Specularity[White, 30]]


• Very nice! I don't know that LaplacianFilter existed when I first did this. :) Commented Feb 27, 2014 at 16:36
• Er… can you add some instruction for the LaplacianFilter? For me the usage of the 2nd argument isn't so clear. Commented Feb 28, 2014 at 8:54
• @xzczd, LaplacianFilter implements the Discrete Laplace Operator with a convolution. The second argument of LaplacianFilter just determines the size of the convolution kernel. By setting this to 1 the smallest possible kernel is used (3x3) which keeps lapMatrix as sparse as possible. If you use a larger kernel you will get a bit more accuracy but slow down the eigenvector calculation. Commented Feb 28, 2014 at 13:01

I had this laying around from a course in numerical linear algebra I taught a few years ago.

Here's a matrix whose nonzero elements describe the basic shape.

size = 50;
nw = Partition[Table[i, {i, 1, size^2}], size];
sw = Partition[Table[i, {i, size^2 + 1, 2*size^2}], size];
se = Partition[Table[i, {i, 2*size^2 + 1, 3*size^2}], size];
L = ArrayFlatten[{{nw, 0}, {sw, se}}];
{m, n} = Dimensions[L];
L = Join[{Table[0, {n + 2}]},
Map[Join[{0}, #, {0}] &, L],
{Table[0, {n + 2}]}];


If size=4, for example, we get something like so

(* size = 4 *)
L // MatrixForm


Now, we convert this to the Laplacian.

entries[0, {_, _}] = {};
entries[k_, {i_, j_}] := Module[{goodVals},
goodVals = Select[{
{i + 1, j}, {i, j + 1}, {i - 1, j}, {i, j - 1}},
L[[Sequence @@ #]] =!= 0 &];
{k, #} -> 1 & /@ (L[[Sequence @@ #]] & /@ goodVals)];
saRules = Flatten[{Band[{1, 1}] -> -4,
MapIndexed[entries, L, {2}]}];
lap = SparseArray[saRules];


lap is a matrix whose dimension is equal to the number of non-zero entries of L. You're essentially looking for an eigenvector of lap.

evs = Reverse[Eigenvectors[N[lap], 8,
Method -> {"Arnoldi", "Shift" -> 0}]];


Here's a simple visualization

ev = evs[[1]];
vib = Map[If[# > 0, ev[[#]], 0] &, L, {2}];
ListPlot3D[vib, ViewPoint -> {2.3, 2.26, 1}]


OK, let's spruce it up a bit

{m, n} = Dimensions[vib];
delete[0, {i_, j_}] := If[
(i + 1 > m || vib[[i + 1, j]] == 0 ) &&
(j + 1 > n || vib[[i, j + 1]] == 0) &&
(i - 1 < 1 || vib[[i - 1, j]] == 0) &&
(j - 1 < 1 || vib[[i, j - 1]] == 0) &&
((i + 1 > m || j + 1 > n) || vib[[i + 1, j + 1]] == 0) &&
((i - 1 < 1 || j + 1 > n) || vib[[i - 1, j + 1]] == 0) &&
((i + 1 > m || j - 1 < 1) || vib[[i + 1, j - 1]] == 0) &&
((i - 1 < 1 || j - 1 < 1) || vib[[i - 1, j - 1]] == 0),
Null, 0];
delete[x_, {_, _}] := x;
vib2 = MapIndexed[delete, vib, {2}];
ListPlot3D[vib2,
Mesh -> None, PlotStyle -> {RGBColor[0.8, 0.2, 0.2],
Specularity[White, 20]}, Lighting -> {
{"Directional", RGBColor[0, 0.4, 0.4],
{{-40, -100, 2}, {0, 0, 0}}},
{"Directional", RGBColor[0.4, 0.4, 0],
{{80, -50, 20}, {0, 0, 0}}}}, Background -> Black,
Boxed -> False, ViewPoint -> {2.3, 2.26, 1}]


• The link you add is instructive for the understanding of the code, especially this and this part :) Commented Feb 28, 2014 at 8:49
• @xzczd Thanks! Glad you liked it. Commented Feb 28, 2014 at 13:23

I had difficulties in replicating the coloring and lighting in Mathematica from the code in logo.m (otherwise, I'm ignoring their "artistic license" of taking only the first two terms of the series solution); nevertheless, I wanted to demonstrate the use of NDEigensystem[]. Thus:

dom = RegionDifference[Rectangle[{-1, -1}, {1, 1}], Rectangle[{-1, 0}, {0, 1}]];
{{val}, {fun}} = NDEigensystem[{Laplacian[u[x, y], {x, y}],
DirichletCondition[u[x, y] == 0, True]},
u, {x, y} ∈ dom, 1,
Method -> {"Eigensystem" -> "Direct",
"Interpolation" -> {"ExtrapolationHandler" ->
{(0 &), "WarningMessage" -> False}},
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.01}},
"VectorNormalization" -> None}];


Check the eigenvalue:

val
-9.657669292805231


The value is a bit different from the $\approx 9.639724$ obtained by Moler; I'm not sure why.

Now, the plot:

Plot3D[fun[x, y], {x, -1, 1}, {y, -1, 1}, Axes -> None, Background -> Black,
BoundaryStyle -> None, Boxed -> False, BoxRatios -> {1, 1, 0.6},
Lighting -> {{"Ambient", White},
{"Directional", RGBColor[0.8, 0.8, 0.], {.5, -1., .4}},
{"Point", RGBColor[0., 0.8, 0.8], {40, 100, 20}}}, Mesh -> None,
PlotPoints -> 35,
PlotStyle -> Directive[RGBColor[1., 0.1, 0.], Specularity[0.4, 7]],
ViewPoint -> {-2.4, -2.4, 1.3}]


As it turns out, this example is in the documentation for NDEigensystem[].

• You need 2 days of calculations to get this number right, see in comments, still it is bad it is even more off in 13.1, -9.66001. blogs.mathworks.com/cleve/?p=1129#reply_10746 Commented Jul 27, 2022 at 17:29
• Commented Jul 27, 2022 at 17:53
• @Валерий good find! (Actually calculating it using Mathematica seems like a stiff challenge, tho.) Commented Jul 27, 2022 at 17:56