18
$\begingroup$

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.

enter image description here

The result was very good, finding the image below.

enter image description here

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.

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$
    – Sjoerd C. de Vries
    Sep 11, 2014 at 19:03

2 Answers 2

Reset to default
26
$\begingroup$

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]

Mathematica graphics

Improve this answer
$\endgroup$
16
  • $\begingroup$ 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. $\endgroup$
    – Santi Carmesí
    Sep 11, 2014 at 22:17
  • $\begingroup$ 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. >> $\endgroup$
    – Santi Carmesí
    Sep 11, 2014 at 22:19
  • $\begingroup$ What's going on? I'll be even more grateful. $\endgroup$
    – Santi Carmesí
    Sep 11, 2014 at 22:20
  • 1
    $\begingroup$ @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]. $\endgroup$
    – Jens
    Sep 12, 2014 at 3:33
  • $\begingroup$ @SantiCarmesí The gradient field issue is discussed in great detail here. $\endgroup$
    – Mark McClure
    Sep 12, 2014 at 9:15
17
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.