# Plot Matlab icon

I started to explore this on a whim and hasn't succeeded yet… Some introduction for the icon is found here but I can't understand it very well. (I admit that, though playing with NDSolve a lot, I suffered 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. –  Szabolcs Feb 27 '14 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. –  Szabolcs Feb 27 '14 at 4:56
My favourite reference in this site on solving elliptic equations is here and very informative IMO –  gpap Feb 27 '14 at 10:32

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. –  Mark McClure Mar 4 '14 at 10:42
I must admit I ignored the stuff about how they actually created the logo, finding the finite difference solution more interesting! –  Simon Woods Mar 4 '14 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. :) –  Mark McClure Feb 27 '14 at 16:36
Er… can you add some instruction for the LaplacianFilter? For me the usage of the 2nd argument isn't so clear. –  xzczd Feb 28 '14 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. –  Simon Woods Feb 28 '14 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 :) –  xzczd Feb 28 '14 at 8:49
@xzczd Thanks! Glad you liked it. –  Mark McClure Feb 28 '14 at 13:23