4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Are you solving the equation in a quarter infinite region or just a square? $\endgroup$ – xzczd Dec 5 '15 at 12:54
  • 1
    $\begingroup$ @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)) $\endgroup$ – Dr. Wolfgang Hintze Dec 5 '15 at 18:20
  • $\begingroup$ @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? $\endgroup$ – Yilun Dec 6 '15 at 0:28
  • $\begingroup$ @xzczd It is a square with semi-infinite condition at x==0 && y==0 $\endgroup$ – Yilun Dec 6 '15 at 0:29
  • 1
    $\begingroup$ @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. $\endgroup$ – Dr. Wolfgang Hintze Dec 6 '15 at 18:08
6
$\begingroup$

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

enter image description here

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

enter image description here

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.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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