Mathematica 10.3 Dsolve solving 2D diffusion PDE

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.

• Are you solving the equation in a quarter infinite region or just a square? Commented Dec 5, 2015 at 12:54
• @Yilun I'm not sure if the boundary and Initial conditions are consistent. Here's a solution which seems to come close to the one requested except for small t: c = 1 - ((-1 + E^(x^2/t))*(-1 + E^(y^2/t))*Pi)/(E^((x^2 + y^2)/t)*(4*t)) Commented Dec 5, 2015 at 18:20
• @Dr.WolfgangHintze Thanks for the solution. I tried the 1D case (see added above) and it is working fine. Might be some limit on Mathematica? Commented Dec 6, 2015 at 0:28
• @xzczd It is a square with semi-infinite condition at x==0 && y==0 Commented Dec 6, 2015 at 0:29
• @xzczd Thank you for your question and your suggestion. I have extended my comment into an answer. Today, in EDIT #1, I gave the correct solution (wiich is different from the one in the comment). BTW, in my version 10.2 MMA can't solve the 1D case as was done in the OP. Commented Dec 6, 2015 at 18:08

1 Answer

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.