# How do I solve a PDE with a strange boundary condition?

How do I solve the PDE with boundary value like this

$$u(t,x,y)=0, \textrm{when } F(x,y)=0$$

using DSolve?

As a specific example, I want to solve heat equation

$$\frac{\partial u}{\partial t}=\alpha\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)$$

with a strange boundary condition $u(t,x,y)=0$ where $x^6+y^4=1$. (Or I want to solve it numerically with a strange closed curve.)

• Could you make a specific, concrete example ? – b.gates.you.know.what Feb 2 '13 at 13:39
• @b.gatessucks OK,I will edit it. – Golbez Feb 2 '13 at 13:42
• Have you tried NDSolve for a numerical solution? – user21 Feb 2 '13 at 20:29

EDIT from the author of this answer, 20 March 2016

This old answer was preceding Mathematica Version 10 which affords the finite Element Method. So this answer is now obsolete. See other answers below. The only interesting information that is remaining is that the Method of lines associated with the TensorProductGrid discretization can eventually handle PDEs in a inhomogeneous medium.

end of the EDIT

Here is a solution for your specific example.
I didn' t succeed in finding a more general solution.

The idea is to use NDSolve.
As NDSolve accept only rectangular boundary conditions, I take a rectangular boundary that encloses your domain, I make $u(t, x, y) = 0$ on this new boundary, and I make the thermal conductivity between the old and the new boundary very high.

The new problem is therefore of the kind :
Heat equation in a non-homogeneous media with rectangular boundary

After many attempts, I succeed in implementing a non-homogeous conductivity equal to :
- $\approx1$ if $(x^4 + y^6) < 1$
- $\gg1$ otherwise

Any improvement would be very appreciated

sol = NDSolve[
With[{thermalConductivity = (1 + (x^4 + y^6)^4)},
{(* heat equation *)
D[u[x, y, t], t] ==
D[thermalConductivity * D[u[x, y, t], x], x]
+ D[thermalConductivity * D[u[x, y, t], y], y],
(* boundary conditions *)
u[-1.5, y, t] == 0,
u[1.5, y, t] == 0,
u[x, -1.5, t] == 0,
u[x, 1.5, t] == 0,
(* initial condition *)
u[x, y, 0] == (2 Exp[-8 (x^2 + y^2)] )
}],
{u},
{t, 0, 0.13}, {x, -1.5, 1.5}, {y, -1.5, 1.5},
AccuracyGoal -> 2,
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {
"TensorProductGrid",
"MinPoints" -> {30, 30}
},
"TemporalVariable" -> t}]


the solution at $t=0.12$ :

Plot3D[Evaluate[u[x, y, 0.12] /. sol], {x, -1.5, 1.5}, {y, -1.5, 1.5},
PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}, {-.1, .5}},
MeshFunctions -> {#3 &},
Mesh -> {{0.02, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5}}, PlotPoints -> 30] The foot of the bump is the shape of $x^4+y^6$, ie the boundary $u(x,y,t)=0$ :

ContourPlot[x^4 + y^6 == 1, {x, -1.5, 1.5}, {y, -1.5, 1.5}] • I'd like this better if you used Piecewise for the conductivity (not sure if it works, but worth a try) – george2079 Mar 16 '14 at 14:53

After some self-learning of finite element method (FEM), I'm able to write a solver for OP's question without NDSolve now. The following code is highly inspired by ruebenko's excellent answer here. To understand the code, you must have some basic understanding for FEM, which anyone with a conscience won't say it's trivial 囧. The best introduction for FEM I found so far is in Chapter 10, Chapter 11, Appendix B of Olver's Introduction to Partial Differential Equations, which is no longer available in his site, but actually Chapter 11 of it still exists as Chapter 14 of his lecture notes on Numerical Analysis.

The finite element equations for Dirichlet problem of heat conduction equation ($u_t=\alpha \text{$\Delta $u}+f$ in $\Omega$, $u=g$ on $\partial \Omega$) is as follows:

$$C \frac{d\phi }{dt}+K \phi =P$$ where

$$K_{i j}=\sum _e K_{i j}^e$$ $$C_{i j}=\sum _e C_{i j}^e$$ $$P_i=\sum _e P_{i}^e$$ $$K_{i j}^e=\int _{\Omega ^e}\alpha \left(\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}+\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z}\right)d\Omega$$ $$C_{i j}^e=\int _{\Omega ^e}N_iN_jd\Omega$$ $$P_i^e=\int _{\Omega ^e}fN_id\Omega$$

The first trouble we encountered when implementing FEM in Mathematica is how to discretize the whole area and index the elements. Ruebenko did mention a post that gives several pieces of code for meshing, but when writing this answer, what's in my mind is mainly how to make use of existed tools in Mathematica as much as possible: does 2D meshing algorithm already integrated in Mathematica? Of course! It's used every time we're ploting:

Clear["*"]
grid = RegionPlot[x^6 + y^4 <= 1, {x, -1.5, 1.5}, {y, -1.5, 1.5}(*, PlotPoints -> 100*)] /.
{a_, b_, c_, d_} -> Sequence[{a, b, c}, {c, d, a}]; // AbsoluteTiming

coord = First@Cases[grid, GraphicsComplex[a_, __] -> a];
index = Flatten[Cases[grid, Polygon[a_] -> a, Infinity], 1];
bound = First@Cases[grid, Line[a__] -> a, Infinity];

grid /. Polygon[a__] -> Sequence[EdgeForm[Black], Polygon[a]] If you feel confused with this piece of code, just check the InputForm of grid. (You can check it by selecting the cell grid exists and pressing Ctrl+Shift+I. )

Then we need to get the stiffness of elements ($K_{i j}^e$). For this part Ruebenko has showed an approach that's smart but somewhat hard to understand and extend to the calculation of $C_{i j}^e$ so I decide to create my own way. The following function is almost a interpretation of the introduction in Olver's lecture note (Though not as good as Ruebenko's, it's also fast :D):

w[x_, y_] = a x + b y + c;
wl[{x_, y_}, {{xi_, yi_}, {xj_, yj_}, {xl_, yl_}}] =
w[x, y] /. Solve[w[xl, yl] == 1 && w[xi, yi] == 0 && w[xj, yj] == 0, {a, b, c}][];
grad[{{xi_, yi_}, {xj_, yj_}, {xl_, yl_}}] =
{D[#, x], D[#, y]} &@ wl[{x, y}, {{xi, yi}, {xj, yj}, {xl, yl}}];

eindex = {i0, i1, i2};
ecoord = {{x0, y0}, {x1, y1}, {x2, y2}};
argu1 = {#, _Integer} & /@ eindex;
argu2 = {#, _Real} & /@ Flatten@ ecoord;
help1 = Tuples[eindex, 2];
allpairs = Map[Join[Complement[ecoord, {#1}], {#1}] &, Tuples[ ecoord, 2], {-2}];
help2 = With[{α = 1},
getpos = With[{help = help1}, Compile[#, help] &@argu1];
getstiffval = With[{help = help2}, Compile[#, help] &@argu2];


Here I stole a function from ruebenko's answer (notice I've made a tiny change):

matrixAssembly[rule_, dim_] :=
Block[{matrix},
SystemSetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
matrix = SparseArray[rule, dim];
SystemSetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}];
matrix]


OK, let's assemble the elements and get $K$:

coordflattened = Flatten /@ (coord[[#]] & /@ index);
AbsoluteTiming[stiffpos = getpos @@@ index;]
AbsoluteTiming[stiffval = getstiffval @@@ coordflattened;]
AbsoluteTiming[totalstiff = matrixAssembly[Thread[stiffpos -> stiffval], Max@index];]


Next part is to get the capacity matrix $C$. The general steps are similar with the building of $K$, but here the calculating of $\int _{\Omega ^e}N_iN_jd\Omega$ turned out to be a little troublesome: at beginning I wanted to make use of the concise syntax of Integrate with Boole for irregular area, but sadly found it doesn't work for this integration, then I turned to NIntegrate, only to find it's a little too slow and the formula for the determination of whether a point is inside a triangle or not provided by MathWorld is wrong. Finally I had to analyse the integration on a scratch paper and modified the integration into regular forms and turned back to Integrate, and the integration didn't finish in an acceptable period of time until I tried different assumptions for times:

int = wl[{x, y}, #] wl[{x, y}, #2] & @@@ allpairs;
AbsoluteTiming[
help3 = With[{ymin = ((y2 - y0) (x - x0))/(x2 - x0) + y0,
ymax1 = ((y1 - y0) (x - x0))/(x1 - x0) + y0,
ymax2 = ((y2 - y1) (x - x1))/(x2 - x1) + y1},
Simplify[(Abs[
Integrate[#1, {x, x0, x1}, {y, ymin, ymax1}, Assumptions -> x0 < x1 < x2] +
Integrate[#1, {x, x1, x2}, {y, ymin, ymax2},Assumptions -> x0 < x1 < x2]] &) /@
int, (y0 | y1 | y2) ∈ Reals]]]


The integration takes about 8 minutes in my computer, if you don't want to wait, just use the following result:

{1/12 Abs[((x0 - x1) (x0^2 + x0 x1 + x1^2 - 3 (x0 + x1) x2 +
3 x2^2) (x2 (y0 - y1) + x0 (y1 - y2) + x1 (-y0 + y2)))/(x0 - x2)^3],
1/24 Abs[x2 (-y0 + y1) + x1 (y0 - y2) + x0 (-y1 + y2)],
1/24 Abs[x2 (-y0 + y1) + x1 (y0 - y2) + x0 (-y1 + y2)],
1/24 Abs[x2 (-y0 + y1) + x1 (y0 - y2) + x0 (-y1 + y2)],
Piecewise[{{(1/12) Abs[x2 (-y0 + y1) + x1 (y0 - y2) + x0 (-y1 + y2)],
(y1 >= 0 && y1 (x0 y1 - x2 y1 + Sqrt[((-x1) y0 + x2 y0 - x0 y2 + x1 y2)^2]) < 0) ||
(y1 < 0 && y1 ((-x0) y1 + x2 y1 + Sqrt[((-x1) y0 + x2 y0 - x0 y2 + x1 y2)^2]) > 0)}}],
1/24 Abs[x2 (-y0 + y1) + x1 (y0 - y2) + x0 (-y1 + y2)],
1/24 Abs[x2 (-y0 + y1) + x1 (y0 - y2) + x0 (-y1 + y2)],
1/24 Abs[x2 (-y0 + y1) + x1 (y0 - y2) + x0 (-y1 + y2)],
1/12 Abs[((x1 - x2) (3 x0^2 + x1^2 + x1 x2 + x2^2 -
3 x0 (x1 + x2)) (x2 (-y0 + y1) + x1 (y0 - y2) + x0 (-y1 + y2)))/(x0 - x2)^3]}


Just a friendly reminder: don't imagine to merge the assumption in Simplify into Integrate.

OK, now let's get $C$:

getcapacityval = With[{help = help3}, Compile[#, help] &@argu2];
indexordered = MapThread[#[[#2]] &, {index, Ordering /@ (coord[[#]] & /@ index)}];
coordordered = Sort /@ (coord[[#]] & /@ index);
capacitypos = getpos @@@ indexordered;
capacityval = getcapacityval @@@ (Flatten /@ coordordered);
totalcapacity = matrixAssembly[Thread[capacitypos -> capacityval], Max@index];


With similar steps, we piecefully get $P$:

load[{x_, y_}] =0.;
i -> load /@ co Det[{{1, 1, 1}, Sequence @@ Transpose[co]}]/6]


The last step is to set the initial condition and boundary condition:

ic[{x_, y_}] = Exp[-2 (x^2 + y^2)];
bc = ConstantArray[0., Length@bound];


Again, I stole a function from ruebenko's answer (notice the tiny modification):

dirichletBoundary[matrix2_, load2_, fPos_List, fValue_List] :=
matrix[[All, fPos]] = matrix[[fPos, All]] = 0.;
matrix += SparseArray[Transpose[{fPos, fPos}] -> Table[1., {Length[fPos]}],
Dimensions[matrix], 0];


Here I choose backward difference for the discretization of $\frac{d\phi }{dt}$ to avoid possible stability problem:

sol = With[{dt = 0.005},
NestList[
LinearSolve @@ (dirichletBoundary[totalstiff, totalload, bound, bc] +
{totalcapacity/dt, #.totalcapacity /dt}) &, ic /@ coord, 25]]; // AbsoluteTiming


This step takes 2.462 seconds in my 2GHz laptop with PlotPoints -> 100 in RegionPlot while with same i.c. and time domain andre's solution takes 16.316 seconds but I won't tell you the code above wastes me almost 20 hours.

Finally, a demonstration of the fruits:

<< Developer
(* Here I stole the third piece of code from ruebenko. *)
Graphics[GraphicsComplex[coord, Polygon[index],
VertexColors -> ToPackedArray@(List @@@ (ColorData["TemperatureMap"][#] & /@ #))]] & /@ sol
Export["my20.gif", %]
ListPlot3D[Flatten /@ Transpose[{coord, #}], PlotRange -> {0, 1},
Mesh -> None, ColorFunctionScaling -> False, ColorFunction -> "TemperatureMap"] & /@ sol
Export["hours.gif", %]  • nice. Stay tuned rumor has it FEM capability is coming soon to NDSove...reference.wolfram.com/search/?q=fem – george2079 Mar 16 '14 at 14:58
• @george2079 Wow, this time it seems to be true! If so, it'll be the best news for me about v10 :D – xzczd Mar 17 '14 at 8:17
• Great work. :) Thank you! – Nam Nguyen Apr 11 '14 at 13:06
• @xzczdYou worked this problem from a robust fundamental theory put forth by Boris Galerkin!! Bad news for the folks at mathematica ;).. for you'll put them out of business ;)!!! With your permission, I would like to demonstrate your answer to my undergraduate class (that I teach) on FEM. – dearN Jan 16 '16 at 13:07
• @LCFactorization Theoretically, it can, but I won't say it's easy. If you just want to solve a PDE with FEM, I'll recommend you to use NDSolve in v10, my solution won't surpass it, even when dealing with nonlinear cases. (I guess you've already read e.g. this post.) – xzczd Oct 13 '16 at 10:56

Using V10's new FEM functionality, this problem can be solved as follows

<< NDSolveFEM;
omega = ImplicitRegion[x^6 + y^4 <= 1, {x, y}];
mesh = ToElementMesh[omega,
"MaxCellMeasure" -> {"Area" -> 0.005, "Length" -> 0.1}];
gamma = DirichletCondition[u[t, x, y] == 0, x^6 + y^4 == 1];
u = NDSolveValue[
{Derivative[1, 0, 0][u][t, x, y] ==
Derivative[0, 0, 2][u][t, x, y] + Derivative[0, 2, 0][u][t, x, y],
gamma, u[0, x, y] == 1 - (x^6 + y^4)}, u,
Element[{x, y}, mesh], {t, 0, 1}] We can then visualize the solution as follows:

pics = Table[
Quiet[Plot3D[u[t, x, y], {x, -1, 1}, {y, -1, 1},
MeshFunctions -> ({#3 &}), PlotPoints -> 101,
RegionFunction -> Function[{x, y}, x^6 + y^4 <= 1],
BoundaryStyle -> Thick, ColorFunction -> "TemperatureMap",
ColorFunctionScaling -> False, PlotRange -> {0, 1.01}]],
{t, 0, 1, 0.01}];
ListAnimate[pics] Note that mesh is an ElementMesh object. It has some a "Wireframe" property that can be accessed directly for visualization.

mesh["Wireframe"] I don't know that there is a simple, built in way to view the plot using the mesh directly, but it can be done.

elements = #[[1 ;; 3]] & /@ mesh["MeshElements"][[1, 1]];
pts = mesh["Coordinates"];
pts3D[t_] := {#[], #[], u[t, #[], #[]]} & /@ pts
vertexColor[t_Real, k_Integer] :=
ColorData["TemperatureMap"][u[t, ##] & @@ pts[[k]]];
Graphics3D[GraphicsComplex[pts3D[0.05],
{Polygon[elements,
VertexColors -> Map[vertexColor[0.05, #] &, elements, {2}]]}
], PlotRange -> {0, 1}] I'd like to add to Mark's answer but do not have enough rep to do so in a comment.

One can write:

Needs["NDSolveFEM"]
omega = ImplicitRegion[x^6 + y^4 <= 1, {x, y}];
m = ToElementMesh[omega];
uif = NDSolveValue[{Derivative[1, 0, 0][u][t, x, y] ==
Derivative[0, 0, 2][u][t, x, y] +
Derivative[0, 2, 0][u][t, x, y],
DirichletCondition[u[t, x, y] == 0, True],
u[0, x, y] == 1 - (x^6 + y^4)}, u, Element[{x, y}, m], {t, 0, 1}];


Note that in the case where the DirichletCondition is valid on all of the boundary one can use True as the predicate.

If one has already computed a mesh then that can be used within the plotting functions:

pics = Table[
Plot3D[uif[t, x, y], {x, y} \[Element] m, BoundaryStyle -> Thick,
ColorFunction -> "TemperatureMap", ColorFunctionScaling -> False,
PlotRange -> {0, 1.01}], {t, 0, 1, 0.1}];
ListAnimate[pics]


I did not use as many frames.

• Most, welcome... – user21 Jul 10 '14 at 19:28
• How do you handle such cases if the boundary consitions are Neumann type zero boundary conditions on the same region? – LCFactorization Feb 11 '16 at 13:55
• Just remove the DirichletCondition. Neumann zero is the default if nothing is specified. – user21 Feb 12 '16 at 0:55