I have a PDE with mixed boundaries (Neumann and Dirichlet on some sides) in the region

$(t,x,y) \in \left( 0, T\right) \times\left\{ -L \leq x \leq L, 0 \leq y \leq h(x) \right\}$

where $h(x)$ is something like $\exp\left\{ -(x-x^*)^2\right\}$, doesn't matter. And I have Neumann boundary condition on the curve $\left(x, h(x)\right)$.

How can I provide such region to NDSolve?

  • 1
    $\begingroup$ Perhaps you may find a mapping that transforms that into a rectangle $\endgroup$ – Dr. belisarius Dec 25 '13 at 18:30
  • 1
    $\begingroup$ This related question may help. Basically, go from (t,x,y) to (t,x,u) where u=H corresponds to the boundary y=h(x). $\endgroup$ – Timothy Wofford Dec 26 '13 at 8:29

This problem can be easily solved using V10's new FEM functionality. For concreteness, let's suppose we want to solve the heat equation $$u_t - \Delta u = 0$$ over the region $$\left\{(x,y): -1 \leq x \leq 1, \; 0 \leq y \leq e^{-x^2}\right\}.$$ We'll take the initial temperature distribution to be identically 1, i.e. $u(x,y,0)=1$; we'll suppose the bottom edge is held at 1 while the left and right edges are held at 0, i.e. $u(x,0,t)=1$ and $u(-1,y,t)=u(1,y,t)=0$; and we'll suppose the curved top is insulated, i.e. the normal derivate of $u$ is zero along the curve $(x,e^{-x^2})$.

omega = ImplicitRegion[-1 <= x <= 1 && 0 <= y <= Exp[-x^2], {x, y}];
mesh = ToElementMesh[omega];
gamma1 = DirichletCondition[u[t, x, y] == 0, x == 1 || x == -1];
gamma2 = DirichletCondition[u[t, x, y] == 1, y == 0];
u = NDSolveValue[{D[u[t, x, y], t] - Laplacian[u[t, x, y], {x, y}] == 
  NeumannValue[0, y == Exp[-x^2]], gamma1, gamma2,
  u[0, x, y] == 1}, u, Element[{x, y}, mesh], {t, 0, 3}];

Note that the NeumannValue is specified as part of the differential equation itself.

We can now plot the solution:

pics = Table[Plot3D[u[t, x, y], Element[{x, y}, omega],
  BoundaryStyle -> Thick, ViewPoint -> {2.4, 2.25, 0.9},
  ColorFunction -> "TemperatureMap", ColorFunctionScaling -> False, 
  PlotRange -> {0, 1.01}], {t, 0, 3, 0.05}];

enter image description here

  • 3
    $\begingroup$ This is a comment to Mark's (@MarkMcClure) answer (I do not have enough rep to comment). The Neumann[0,...] is not needed in this case. Neumann zero boundary conditions drop out of the FEM equations and play no role in the solution. If a boundary has no boundary condition prescribed then they are implicitly set to Neumann zero boundary values. This is why they are called 'natural' boundary conditions. $\endgroup$ – user21 Jul 11 '14 at 12:24

NDSolve requires a rectangular domain, so you have to make a change of coordinates.

$(t,x,y) \in \left( 0, T\right) \times\left\{ a \leq x \leq b, g(x) \leq y \leq h(x) \right\}$

Since you have explicit expressions for your boundaries ($x=a$, $x=b$, $y=g(x)$, and $y=h(x)$), we can use a linear interpolation

x[u_, v_] = (1 - u)*a + u*b
y[u_, v_] = With[{x=x[u,v]},(1 - v)*g[x] + v*h[x]]
$Assumptions = {Element[{g[_], h[_]}, Reals], h[a_] > g[a_]}

In terms of these coordinates your domain becomes

$(t,u,v) \in \left( 0, T\right) \times\left\{ 0 \leq u \leq 1, 0 \leq v \leq 1 \right\}$

In my other answer I described how to transform the PDE. The general steps are

  1. write your metric
  2. find its inverse and determinant
  3. find the basis vectors and basis covectors
  4. transform components of vectors into the new coordinate bases
  5. find expressions for operations (directional derivatives, laplacians, etc.)
  6. transform PDE and IV/BVs

NDSolve should then be able to give you a solution in terms of $u$ and $v$.


You then have to map your result onto your original domain.

u[x_,y_]:= (x-a)/(b-a)
v[x_,y_]:= (y-g[x])/(h[x]-g[x])
F[t_,x_,y_]:= sol[t,u[x,y],v[x,y]]
  • $\begingroup$ I responded yesterday with a similar answer but deleted it when I realized that my approach did not deal correctly with the Neumann type boundary condition at the top. The problem is that your transformation does not preserve orthogonality on the top boundary thus, a Neumann condition in the (x,y) plane doesn't translate to a normal derivative in the (u,v) plane. Furthermore, NDSolve doesn't work when non-normal derivatives are are included as boundary conditions, it issues an NDSolve::bcnorm error. I could temporarily un-delete my answer, if you like. $\endgroup$ – Mark McClure Dec 27 '13 at 18:33
  • $\begingroup$ @MarkMcClure, I was wondering why you deleted your answer. I would un-delete it, since it had good info. If it has a mistake, you can edit a comment, and maybe we can fix it. I agree that the equation representing the boundary condition must be changed along with PDE. The Neumann BC is n[x].grad[F[x,h[x]]==bc[x]. In the answer I linked to, I showed how to transform the directional derivative. $\endgroup$ – Timothy Wofford Dec 27 '13 at 19:12
  • $\begingroup$ @MarkMcClure I see what you mean now. We need to find a better coordinate transformation. One in which not only is the coordinate constant along the boundary, but the coordinate basis vector is normal to the boundary. $\endgroup$ – Timothy Wofford Dec 27 '13 at 23:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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