# Comparing analytical solution with numerical solution of Helmholtz equation in a unit square

I am just learning PDE, and I am interested to compare analytical solution with numerical solution of Helmholtz equation in a unit square with zero boundary condition. I am not sure if it possible. Below is the equation I am analyzing, with its eigenvalues and eigenfunctions: $$\begin{cases} \Delta^2u+\lambda u=0,\\ u(0,y)=0,\quad u(x,0)=0,\\ u(1,y)=0,\quad u(x,1)=0, \end{cases}$$ I found out from a textbook that the eigenvalues and eigenfunctions are: $\lambda_{nm}=(n^2+m^2)\pi^2$ and $u_{nm}=\sin n\pi x \sin m\pi y$, where $n,m=1,2,3,\ldots$. I know that the solution to the above equation is trivial, i.e., $u=0$. However, I found from question regarding solution to Helmholtz by other user, that we can solve the equation in term of eigenfunctions and eigenvalues. So, in my case, the first six eigenfunctions are

region = Rectangle[];
{eigenvalues[region], eigenfunctions[region]} = NDEigensystem[{-Laplacian[u[x, y], {x, y}] + u[x, y], DirichletCondition[u[x, y] == 0, True]}, u[x, y], {x, y} \[Element] region, 6];
Grid[Partition[Table[Show[{ContourPlot[eigenfunctions[region][[j]], {x, y} \[Element] region, Frame -> None, ColorFunction -> "BlueGreenYellow", PlotPoints -> 60, PlotRange -> Full, PlotLabel -> eigenvalues[region][[j]]],RegionPlot[region, PlotStyle -> None, BoundaryStyle -> {Black, Thick}]}], {j, 1, Length[eigenvalues[region]]}], 3]]


I wonder if I can do error analysis by suming the eigenfunctions and compare it with $u=0$ in $L_2$ norm.

1) But, what is the exact analytical solution of this equation in terms of eigenfunctions? I know it is the summation of those eigenvalues, but I am not sure how to sum it.

2) And, how to sum those first six images of eigenfunctions with mathematica? Once somebody show me how to sum those images, I think I can do the error analysis myself. Thanks for your time!

• I am not exactly sure I understand what you are trying to do. Maybe you are looking for an eigenfunction expansion? There is an example of this in the documentation of NDEigensystem in the section Applications and then Eigenfunction Expansion. Hope this helps. May 4, 2017 at 3:13
• I think the solution to the above equation is trivial, which is $u=0$. However, if I am not wrong we can also get solution by adding the eigenfunctions? So I plan to do error analysis, for example I want to compare the exact analytical solution, which is $u=0$ with the sum of eigenfunctions. I assume the more eigenfunctions we add, the smaller the error is? I am just learning PDE, so that is what my understanding is. I am not sure if it is correct. So for example, if mathematica gives me these six eigenfunctions, how to add those functions?
– Lila
May 4, 2017 at 3:24
• I think there is some confusion here. $u=0$ is $u_{nm}$ with $n=m=0$. No sum just one term, no error too. However, what is the significance of this solution? The Helmholtz equation describes some physical process and for this matter $u=0$ is to be neglected. May 4, 2017 at 5:03
• The $u$ that I mean is the final analytical solution, i,e,, when the $m,n$ is infinite. Is it possible? Or my understanding is wrong?
– Lila
May 4, 2017 at 5:06
• Note 1 that the solution also depends on initial conditions. You have to supply two. 2 you have chosen a square which due to symmetry will sometimes have identical eigenvalues. See your second and third eigenvectors. These need special treatment when combining them. A rectangle is simpler.
– Hugh
May 4, 2017 at 5:55

Well, I'm not that familiar with this topic, so feel free to point out if I'm wrong.

I think the trivial solution simply cannot be represented by a combination of eigenfunctions, the reason is as follows.

Let's add a constant inhomengeneous term $s$ to the equation

$$\begin{cases} \Delta^2u+\lambda u=s,\\ u(0,y)=0,\quad u(x,0)=0,\\ u(1,y)=0,\quad u(x,1)=0, \end{cases}$$

and solve it with the help of finiteFourierSinTransform. (I don't use DSolve here because it doesn't represent the solution in terms of eigenfunctions and eigenvalues. )

First, interpret the system to Mathematica code:

With[{u = u[x, y]}, eq = Laplacian[u, {x, y}] + lambda u == s;
bc = u == 0 /. List /@ Flatten@Outer[Rule, {x, y}, {0, 1}]]


Then make transform on $x$ and $y$:

ffst = finiteFourierSinTransform;
Format@ffst[f_, __] := Subscript[\[ScriptCapitalF], s][f]

{teq, tbc} = ffst[{eq, bc[[3 ;; 4]]}, {x, 0, 1}, n] /. Rule @@@ bc[[1 ;; 2]]

tteq = ffst[teq, {y, 0, 1}, m] /. ffst[ffst[a_, b__], c__] :> ffst[ffst[a, c], b] /.
Rule @@@ tbc


Finally solve for the transformed solution and transform back:

ttsol = u[x, y] /. First@Solve[tteq /. ffst[ffst[a_, b__], c__] :> a, u[x, y]]

tsol = inverseFiniteFourierSinTransform[ttsol, n, {x, 0, 1}]

sol = HoldForm@Sum[#, {n, C@1}, {m, C@2}] &[
inverseFiniteFourierSinTransform[tsol, m, {y, 0, 1}] //. HoldForm@Sum[a_, __] :> a] Apparently, as long as $s$ is a non-zero number, the particular solution for the problem can be represented by a combination of eigenfunctions and eigenvalues, but when $s=0$, every summand becomes zero. That's the reason why I think eigenfunctions and eigenvalues can't represent this zero solution.

Finally, a image taken at $\lambda=1\,,s=1$:

Plot3D[sol //. {lambda -> 1, s -> 1, C@_ -> 6,
HoldPattern@Sum[a__] :> Total@Flatten@Table@a} // ReleaseHold // Evaluate, {x, 0,
1}, {y, 0, 1}] Since I read somewhere that the solution to the above is infinite sum of eigenfunctions, so add the solutions: starting by adding the first 2 eigenfunctions, then 3, until 6. I notice that after each addition, the image is getting lighter, which means it is closer to the analytical solution $u=0$. So, the more eigenfunctions, the error is smaller. Please correct me if I am wrong. (I don't know how to enter the code using code command, since it involves images).     • Why do you add images, not the functions? Notice also that if $u_{nm}$ is an eigenfunction, so does the rescaled one $\alpha u_{nm}$. They both correspond to the same eigenvalue $\lambda_{mn}$. So why not adding with prefactors? Where does your formula about the summation coming from? May 4, 2017 at 5:09
• I add the images since I thought the image represent the eigenfunctions. I add images since I want to analyze it numerically with the solution from Mathematica, which is the eigenfunctions? So you think it is not possible to use the images from Mathematica? I know the analytical solution, which involves the functions and prefactors. I just want to compare analytical and numerical solution which involves eigenfunctions.
– Lila
May 4, 2017 at 5:20
• This method for adding solutions is wrong, just add the same picture for several times and you'll know what I mean. The correct way is e.g. Plot3D[eigenfunctions[region] // Total // Evaluate, {x, y} \[Element] region], which clearly shows that the naive summation isn't getting closer to 0. Sep 14, 2017 at 17:32