# Solve a PDE on a domain $\Omega$ with given boundary conditions

I'm starting to study the behavior of some PDEs and I would like to run simulations in mathematica to help me visualize solutions.

For example, a prime example that I would like to study is $$\left\{ \begin{array}{cc} -\Delta\phi + u\cdot\nabla\phi=1 & \text{in }\Omega\\ \phi=0& \text{on }\partial\Omega \end{array}\right.$$ for a fixed (smooth enough) $u$ which is divergence free and flows tangentially to the boundary.

For one dimensional problems, inputting boundary conditions is easy since there are typically only two bounary points at most. But I would like to study this PDE (in 2-d) with a variety of different $\Omega$'s. I can't find anything in the mathematica documentation about how to input arbitrary boundary conditions.

To make things simpler, I will restrict myself to reasonably nice $\Omega$'s which have nice boundaries. Suppose that $\partial\Omega$ is given implicitly by $F(x,y)=0$, e.g. $x^2+y^2-1=0$. What would be the easiest/best/most efficient way to input the boundary condition $\phi(x,y)=0$ when $F(x,y)=0$?

Edit: Feel free to use any reasonable discretization technique, or to assume that you have a function like isInRegion[x,y] if you need it.

• This is one of the questions I'm considering :). It's a pity that NDSolve can not yet handle irregular rigion, and this seems to be the only question that tried to deal with this issue in this site so long. Oct 28 '13 at 2:25
• Besides the сhanging of the coordinate system you can use irregular grids. See this answer and link in it. Jan 19 '14 at 21:31
• May be this is also relevant mathematica.stackexchange.com/a/54492/1089 Aug 24 '14 at 14:18

Define your region $\Omega$ and its boundary $\partial\Omega$ by $$\left\{ \begin{array}{cc} f(x,y)>0 & \text{in }\Omega\\ f(x,y)=0& \text{on }\partial\Omega \end{array}\right.$$ And define another function $g(x,y)$ so that the ordered pair $(f(x,y),g(x,y))$ is invertible, i.e. you can write $x=x(f,g)$ and $y=y(f,g)$.

In the $x,y$ coordinate system, the level curves of $f$ form something like concentric circles, though irregularly shaped. So $f$ is something like a radial coordinate. If you choose $g$ to be the streamlines of the gradient of $f$, then the curves of $g$ are orthogonal to the curves of $f$, like the radial rays of constant angle in polar coordinates.

In the $f,g$ coordinate systems, the boundary conditions are rectangular, $0<f<F$ and $0<g<2\pi$ with periodic boundary conditions for the $g$ dimension: $\phi(f,0)=\phi(f,2\pi)$ and your given boundary condition for the $f$ dimension: $\phi(0,g)=0$ along with the physical condition $\phi(F,g)<\infty$.

You can look up expressions for the laplacian and directional derivative in general curvilinear coordinates, or derive it yourself. Your PDE will now look really complicated, but the boundaries are rectangular, so NDSolve should be able to handle them.

EDIT: The only really hard part is coming up with the coordinate transformations. You may find acceptable expressions using Assumptions with FullSimplify on the results of Solve. Here is an example

(* our new coordinates *)
c = {\[ScriptF], ℊ};


The boundary defines one family of curves. $F(x,y)=f$

(* Begin example *)
\[GothicF][x_, y_] := 1 - (x^2 + y^2)


The other family of coordinate curves is arbitrary, but needs to be compatible. $G(x,y)=g$

\[GothicG][x_, y_] := ArcTan[x, y]


Try to let the computer do the hard work of inverting the two equations above. $x=X(f,g)$ and $y=Y(f,g)$

{x[\[ScriptF]_, ℊ_], y[\[ScriptF]_, ℊ_]} = {x, y} /.
Assuming[(0 < \[ScriptF] < 1) &&
(0 <= ℊ < 2*Pi),
FullSimplify@First@Solve[{
\[GothicF][x, y] == \[ScriptF],
\[GothicG][x, y] == ℊ}, {x, y}]]
(* End example *)


Once you've decided on your coordinates and the transformation, the rest is tedious and algorithmic. First we find the coefficients of the metric tensor and its inverse.

$$ds^2=g_{\text{ij}}dx^idx^j$$

dsSquared = Expand[(#1 . #1 & )[Dt[{x @@ c, y @@ c}]]];
g = FullSimplify@Table[
(1/2 + 1/2 KroneckerDelta[i, j])Coefficient[dsSquared, Dt[c[[i]]]*Dt[c[[j]]]],
{i, Length[c]}, {j, Length[c]}];
gI = Inverse[g];


Then we are ready to calculate the Laplace-Beltrami operator

$$\Delta \phi =\frac{1}{\sqrt{|g|}}\frac{\partial }{\partial x^i}\left(g^{\text{jk}}\sqrt{|g|}\frac{\partial \phi }{\partial x^k}\right)$$

Laplacian[ϕ_, c_, g_] :=
FullSimplify@Expand@Sum[
(1/Sqrt[Abs[Det[g]]])*D[Sqrt[Abs[Det[g]]]*gI[[j, k]]*D[ϕ @@ c, c[[k]]], c[[j]]],
{j, Length[c]}, {k, Length[c]}]


To handle the directional derivative, we first need to find the components of the vector in our new coordinate system. $$\overset{\rightharpoonup }{u}=u^x\overset{\rightharpoonup }{e}_x+u^x\overset{\rightharpoonup }{e}_y=u^f\overset{\rightharpoonup }{e}_f+u^g\overset{\rightharpoonup }{e}_g$$ The basis vectors $\overset{\rightharpoonup }{e}_f,\overset{\rightharpoonup }{e}_g$

eSub = {D[x @@ c, c[]], D[y @@ c, c[]]};
eSub = {D[x @@ c, c[]], D[y @@ c, c[]]};


The basis covectors $\overset{\rightharpoonup }{e}^f,\overset{\rightharpoonup }{e}^g$

eSup = (eSub*g[[2, 2]] - eSub*g[[2, 1]])/Det[g];
eSup = (g[[1, 1]]*eSub - g[[1, 2]]*eSub)/Det[g];


The components of the vector are then $$u^f=\overset{\rightharpoonup }{e}^f\cdot \overset{\rightharpoonup }{u}=\overset{\rightharpoonup }{e}^f\cdot \overset{\rightharpoonup }{e}_xu^x+\overset{\rightharpoonup }{e}^f\cdot \overset{\rightharpoonup }{e}_yu^y$$

{uf, ug} = {{eSup . {1, 0}, eSup . {0, 1}},
{eSup . {1, 0}, eSup . {0, 1}}} . {ux, uy};


And the directional derivative is $$\overset{\rightharpoonup }{u}\cdot \nabla \phi =u^f\frac{\partial }{\partial f}\phi +u^g\frac{\partial }{\partial g}\phi$$

DirectionalDerivative[ϕ_, c_, u_] := Sum[
u[[i]]*D[ϕ @@ c, c[[i]]],
{i, Length[c]}]


Finally, we construct your differential equation

eqn = -Laplacian[ϕ, c, g] + DirectionalDerivative[ϕ, c, {uf, ug}] == 1

• I'm not quite satisfied with this answer. It pushes all the hard work onto me finding these change of coordinates and what the new PDE looks like. I would like to push all the work onto the computer. Nov 1 '13 at 12:46
• @nullUser, It doesn't push the hard work onto you. Just go ahead and push the hard work onto the computer. (See my edit.) Since Mathematica doesn't accept coordinate-free notation, you will have to set up your PDE with respect to a coordinate system anyway. It shouldn't be any more onerous to use a general curvilinear coordinate system than any of the standards (cartesian, polar, elliptic, parabolic, hyperbolic, etc). Nov 1 '13 at 16:33
• very good answer ! Aug 24 '14 at 13:12