Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.)

share|improve this question
2  
Could you make a specific, concrete example ? –  b.gatessucks 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

4 Answers 4

up vote 12 down vote accepted

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 enclose 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 :
- ~1 if (x^4 + y^6) < 1
- >>1 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]

enter image description here

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}]

enter image description here

share|improve this answer
    
I'd like this better if you used Piecewise for the conductivity (not sure if it works, but worth a try) –  george2079 Mar 16 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]]

enter image description here

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", %]

enter image description here

enter image description here

share|improve this answer
    
nice. Stay tuned rumor has it FEM capability is coming soon to NDSove...reference.wolfram.com/search/?q=fem –  george2079 Mar 16 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 at 8:17
    
Great work. :) Thank you! –  DaoTRINH Apr 11 at 13:06

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}]

enter image description here

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]

enter image description here

Note that mesh is an ElementMesh object. It has some a "Wireframe" property that can be accessed directly for visualization.

mesh["Wireframe"]

enter image description here

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}]

enter image description here

share|improve this answer

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["NDSolve`FEM`"]
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.

share|improve this answer
2  
Thanks - mystery man! –  Mark McClure Jul 10 at 19:27
    
Most, welcome... –  user21 Jul 10 at 19:28

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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