Mathematica 10.3 Dsolve solving 2D diffusion PDE

Question

I tried to use Mathematica 10.3 DSolve to solve 2D diffusion PDE with Dirichlet boundary conditions. But it seems not working. Could anyone help check the code?

DSolve[
    {
        D[c[x, y, t], t] - D[c[x, y, t], x, x] - D[c[x, y, t], y, y] == 0,
        c[x, y, 0] == 0, 
        c[x, 0, t] == 1, c[0, y, t] == 1
    },
    c[x, y, t], {x, y, t}
]

But for 1D case, it is working perfectly:

DSolve[
    {
        D[c[x, t], t] - D[c[x, t], x, x] == 0,
        c[x, 0] == 0, 
        c[0, t] == 1
    },
    c[x, t], {x, t}
]

I tried a 1D case with time-dependent Dirichlet boundary condition:

DSolve[
    {
        D[c[x, t], t] - DC*D[c[x, t], x, x] == 0,
        c[x, 0] == 0, 
        c[0, t] == 1 - Exp[-β*t]
    },
    c[x, t], {x, t}
]

But the solution I got does not meet the PDE.

Answer 1

The solution

I have found the correct solution to your problem:

u3[x_, t_] = 1 - (Erf[x/(2 Sqrt[t])]) (Erf[y/(2 Sqrt[t])]);

It is a solution of the PDE

D[u3[x, y, t], t] == D[u3[x, y, t], {x, 2}] + D[u3[x, y, t], {y, 2}]

(* Out[38]= True *)

The initial condition is ok:

Limit[u3[x, y, t], t -> 0, Assumptions -> {x > 0, y > 0}]

(* Out[35]= 0 *)

The boundary conditions are ok, as well:

{u3[0, y, t], u3[x, 0, t]}

(* Out[54]= {1, 1} *)

Plotting the solution for two different times

With[{t = 0.1}, 
 Plot3D[u3[x, y, t], {x, 0, 5}, {y, 0, 5}, PlotRange -> {0, 1}]]

With[{t = 1}, 
 Plot3D[u3[x, y, t], {x, 0, 5}, {y, 0, 5}, PlotRange -> {0, 1}]]

How did I find the solution? After determining the general solution (cf. next section) and studying some examples I came across the correct amplitude. Hence the method was not completely constructive.

The derivation

The solution method is the standard separation of variables

Let

c = T[t] X[x] Y[y]

which after substituting into the PDE leads to

T'[t] /T[t] = X''[x]/X[x] + Y''[y]/Y[y]  = - k^2

with -k^2 being the separation constant.

Letting

X''[x]/X[x] = p^2
Y''[y]/Y[y] = q^2 

and requesting c = 0 for x = 0 and y = 0 we find the particular solution

f0 = Exp[- t (p^2 + q^2)] Sin[p x] Sin[q y];

This is a solution for all p and q.

Because of the linearity of the PDE any linear combination of particular solutions is again a solution.

Therefore the general solution can be written as

ua = Integrate[
   a[p, q] Exp[-t (p^2 + q^2)] Sin[p x] Sin[q y], {p}, {q}] // Quiet

with an arbitrary amplitude a[p,q] and some appropiate interval of integration.

I have tried some possible amplitudes and their corresponding solutions

u0[x_, y_, t_] = 
 1 - Integrate[f0  , {p, 0, ∞}, {q, 0, ∞}, 
   Assumptions -> {x > 0, y > 0, t > 0}]

(* Out[501]= 1 - (DawsonF[x/(2 Sqrt[t])] DawsonF[y/(2 Sqrt[t])])/t *)

f1 = Exp[-t (q^2 + p^2)] p q Sin[p x] Sin[q y];

u1[x_, y_, t_] = 
 1 - Integrate[
   f1, {p, -∞, ∞}, {q, -∞, ∞}, 
   Assumptions -> {x > 0, y > 0, t > 0}]

(* Out[505]= 1 - (E^(-((x^2 + y^2)/(4 t))) π x y)/(4 t^3) *)

Taking the the amplitude (2/Pi 1/p)(2/Pi 1/q):

f3 = 
  Exp[-t (q^2 + p^2)] 2/(π p) Sin[p x] 2/(π q) Sin[q y];

u3[x_, y_, t_] = 
 1 - Integrate[f3 , {p, 0, ∞}, {q, 0, ∞}, 
   Assumptions -> {x > 0, y > 0, t > 0}]

(* Out[11]= 1 - Erf[x/(2 Sqrt[t])] Erf[y/(2 Sqrt[t])] *)

Notice: the name u3 was used here because u2 was already occupied, but I detected an error in the "solution" u2 and removed it.