How do I solve a PDE with a strange boundary condition?
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}][[1]];
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},
1/2 Det[{{1, 1, 1}}~Join~Transpose[ecoord]] (α grad[#].grad[#2] &) @@@ allpairs];
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},
System`SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
matrix = SparseArray[rule, dim];
System`SetSystemOptions["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.;
getload[i_] := Block[{co = coord[[i]]},
i -> load /@ co Det[{{1, 1, 1}, Sequence @@ Transpose[co]}]/6]
AbsoluteTiming[totalload = matrixAssembly[getload /@ index, Max@index];]
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] :=
Module[{matrix = matrix2, load = load2},
load -= matrix[[All, fPos]].fValue;
load[[fPos]] = fValue;
matrix[[All, fPos]] = matrix[[fPos, All]] = 0.;
matrix += SparseArray[Transpose[{fPos, fPos}] -> Table[1., {Length[fPos]}],
Dimensions[matrix], 0];
{matrix, load}]
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", %]
Using V10's new FEM functionality, this problem can be solved as follows
<< NDSolve`FEM`;
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_] := {#[[1]], #[[2]], u[t, #[[1]], #[[2]]]} & /@ 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}]
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}]
