# 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. Sep 11, 2014 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. Sep 11, 2014 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. >> Sep 11, 2014 at 22:19
• What's going on? I'll be even more grateful. Sep 11, 2014 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, 2014 at 3:33
• @SantiCarmesí The gradient field issue is discussed in great detail here. Sep 12, 2014 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 Jun 26, 2016 at 19:09
• @Young, looks good to me. Jun 26, 2016 at 23:29
• Is there some global settings that need to be changed from the default? Jun 26, 2016 at 23:30
• @Young, no. But since I do not know/see what the issue is I can not really comment. Jun 26, 2016 at 23:33