# 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? – xzczd Dec 5 '15 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)) – Dr. Wolfgang Hintze Dec 5 '15 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? – Yilun Dec 6 '15 at 0:28
• @xzczd It is a square with semi-infinite condition at x==0 && y==0 – Yilun Dec 6 '15 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. – Dr. Wolfgang Hintze Dec 6 '15 at 18:08

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.