# Solve Laplace equation using NDSolve

I am new to Mathematica, a friend recommended this software and started using it, in fact download the trial version to know.

I recently did a program in C to calculate numerically the solution to the Laplace equation in two dimensions for a set of points as in the figure.

The result was very good, finding the image below.

This program took me about 100 lines in C, my friend told me that Mathematica could do it in a couple of lines, which seemed quite interesting. I studied a bit and found that Mathematica can solve the Laplace and Poisson equations using NDSolve command.

However, this command requires to be given to the specific boundary conditions. The boundary condition in which $\phi = 0$, it is quite easy to introduce. But on the inside border, where $\phi = 100$, I failed to get the condition. My idea is to make this problem as follows:

            uval = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] ==
0,
u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0,
u[40, y] == u[60, y] == u[x, 40] == u[x, 60] == 100},
u, {x, 0, 100}, {y, 0, 100}]


But this does not work, and I could not give the boundary conditions on the inner square. Any help would be the most grateful.

• Aren't there conflicting boundary conditions? For instance, u[x,40] equals u[100,y] at {x,y}=={100,40} making you state 0==100 at that point. – Sjoerd C. de Vries Sep 11 '14 at 19:03

You need a DirichletCondition (new in V10) here. Using regions (also from V10):

Ω = RegionDifference[Rectangle[{0, 0}, {100, 100}], Rectangle[{40, 40}, {60, 60}]];
sol = NDSolveValue[{
D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 100.,
x == 40 && 40 <= y <= 60 ||
x == 60 && 40 <= y <= 60 ||
40 <= x <= 60 && y == 40 ||
40 <= x <= 60 && y == 60],
u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0
}, u, {x, y} ∈ Ω]

DensityPlot[sol[x, y], {x, y} ∈ Ω,  Mesh -> None, ColorFunction -> "Rainbow",
PlotRange -> All, PlotLegends->Automatic]


• Thank you very much. I wondered, with a couple of lines !!! Just one more question friend, if not a lot to ask. How could plot the vector field phi? That is, graphic $-\nabla sol(x,y)$. Mi idea es hacer: field = -Grad[sol[x, y], {x, y}] VectorPlot[campo[x, y], {x, 0, 100}, {y, 0, 100}] But I shows a blank graph. – Santi Carmesí Sep 11 '14 at 22:17
• And I get the messages: > InterpolatingFunction::dmval: Input value {42.8643,42.8643} lies outside the range of data in the interpolating function. Extrapolation will be used. >> General::stop: Further output of InterpolatingFunction::dmval will be suppressed during this calculation. >> – Santi Carmesí Sep 11 '14 at 22:19
• What's going on? I'll be even more grateful. – Santi Carmesí Sep 11 '14 at 22:20
• @SantiCarmesí For the gradient and field strength, you can use my answer here. The error message in the potential plot appears because the inner boundary isn't sampled quite accurately. One can get rid of this by first defining rm=RegionMember[\[CapitalOmega]] and then replacing the plot argument sol[x,y] by If[rm[{x,y}],sol[x,y],100]. – Jens Sep 12 '14 at 3:33
• @SantiCarmesí The gradient field issue is discussed in great detail here. – Mark McClure Sep 12 '14 at 9:15

Building on @SjoerdC.deVries answer you can use:

sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 100.,
x == 40 && 40 <= y <= 60 || x == 60 && 40 <= y <= 60 ||
40 <= x <= 60 && y == 40 || 40 <= x <= 60 && y == 60],
u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0},
u, {x, y} \[Element] \[CapitalOmega],
Method -> {"FiniteElement",
"MeshOptions" -> {"BoundaryMeshGenerator" -> "Continuation"}}]


Adding this option results in a better boundary approximation and you will not get the warning messages about extrapolation at the boundary.

• The results I'm getting seem slightly discontinuous and jagged. Even when I decrease the MaxCellMeasure it doesn't seem to improve – Young Jun 26 '16 at 19:09
• @Young, looks good to me. – user21 Jun 26 '16 at 23:29
• Is there some global settings that need to be changed from the default? – Young Jun 26 '16 at 23:30
• @Young, no. But since I do not know/see what the issue is I can not really comment. – user21 Jun 26 '16 at 23:33