I'm modeling a problem with PDEs, So I gotta solve numerically this Reaction-Diffusion Partial Differential Equation

$$ \frac{\partial u(t,x,y)}{\partial t}=D\Big( \frac{\partial^{2}u(t,x,y) }{\partial x^{2}} +\frac{\partial^{2}u(t,x,y) }{\partial y^{2}}\Big)+r\ u(t,x,y) \,\, , $$ where $D$ is the diffusion coeficient and the $r$ is a reaction coefficient (both constants). My boundary conditions are:

$$ \frac{\partial u(t,0,y)}{\partial x}=\frac{\partial u(t,x,0)}{\partial y}=\frac{\partial u(t,a,y)}{\partial x}=\frac{\partial u(t,x,b)}{\partial y}=0 \,\, , $$

and my inicial condition is:

$$ u(0,x,y)=u_{0}\delta (x) \,\, , $$

Where $\delta (x)$ is the Dirac Delta Function. I'm not able to solve this equation on mathematica yet, How could I solve this equation numerically? (I'm not very sure about the boundary conditions, if there's something worong with them, feel free to tell me.)

The code I've tried to make is this:

It defines the EDP I'm dealing.

eqd = D[u[t, x, y], t] == 
D[u[t, x, y], x, x] + D[u[t, x, y], y, y] + u[t, x, y]

It defines the domain of my problem.

\[CapitalOmega] = ImplicitRegion[0 <= x <= 100 && 0 <= y <= 100, {x, y}];

So I've tried to solve this using NDSolve, where I've changed the Dirac Delta by the constant function 100 as inicial condition, and the boundary conditions I tried to express as functions D[u[t,0,y],x]==0 and so on, this way:

NDSolve[{eqd, D[u[t, x, y], x]/.x->0 == 0, D[u[t, x, y], x]/.x->100 == 0, 
D[u[t, x, y], y]/.y->0 == 0, D[u[t, x, y], y]/.y->100 == 0, 
u[t, x, y] == 100}, u, {x, 0, 100}, {y, 0, 100}, {t, 0, 10}]

What has produced the following message:

"Equation or list of equations expected instead of True in the first \
argument {\!\(\*SuperscriptBox[\"u\", TagBox[
RowBox[{\"1\", \",\", \"0\", \",\", \"0\"}], \")\"}],
MultilineFunction->None][t, x, y] == u[t, x, y] + \
\*SuperscriptBox[\"u\", TagBox[
RowBox[{\"0\", \",\", \"0\", \",\", \"2\"}], \")\"}],
MultilineFunction->None][t, x, y] + \*SuperscriptBox[\"u\", TagBox[
RowBox[{\"0\", \",\", \"2\", \",\", \"0\"}], \")\"}],
MultilineFunction->None][t, x, \

and didn't give me an answer for my problem...

Thank you for supporting!

  • $\begingroup$ Take a look at NDSolve and NeumannValue. $\endgroup$
    – user484
    Nov 4, 2015 at 19:00
  • $\begingroup$ @Rahul I'm Using mathematica 9, this way I don't have the NeumannValue feature avaliable to me... $\endgroup$ Nov 4, 2015 at 19:04
  • 1
    $\begingroup$ NeumannValue is convenient but not necessary to solve this system of equations. To obtain more thorough responses, post your equations in Mathematica format. $\endgroup$
    – bbgodfrey
    Nov 4, 2015 at 21:55
  • $\begingroup$ Note that the solution to this equation will grow exponentially at rate r once the initial distribution spreads out over the domain. $\endgroup$
    – Chris K
    Nov 4, 2015 at 22:10
  • 1
    $\begingroup$ Just wanted to point out the OP's code doesn't work because expressions like /.x->0 == 0 evaluate to /.x->True etc, and the initial condition was wrongly specified as u[t, x, y] == 100 instead of u[0, x, y] == 100. Fixing the latter, and using parentheses to fix replacements, e.g. (D[u[t, x, y], x]/.x->0) == 0 gives a solution equivalent to the one in the accepted answer. $\endgroup$
    – obsolesced
    Jul 18, 2016 at 7:41

1 Answer 1


I note that your diffusion constant has gone away (been set to $1$). I get working solutions with

soln = NDSolve[{
  D[u[t, x, y], t] == 
    D[u[t, x, y], x, x] + D[u[t, x, y], y, y] + u[t, x, y], 
  Derivative[0, 1, 0][u][t, 0, y] == 0,
  Derivative[0, 1, 0][u][t, 100, y] == 0,
  Derivative[0, 0, 1][u][t, x, 0] == 0,
  Derivative[0, 0, 1][u][t, x, 100] == 0,
  u[0, x, y] == 100
}, u, {x, 0, 100}, {y, 0, 100}, {t, 0, 10}]

Plot3D[u[t, x, 1] /. soln, {t, 0, 10}, {x, 0, 100}]

Is there any particular analysis you're interested in?

  • $\begingroup$ Thanks for your Answer! :D Actually There's one, I would like to solve this problem with the inicial condition $u(0,x,y)==u_{0}\delta (x)$, what I can't express in mathematica language. $\endgroup$ Nov 9, 2015 at 1:28
  • $\begingroup$ @WaynerKlën : How is this different from the line u[0,x,y] == 100 (with whatever appropriate constant instead of "100")? $\endgroup$ Nov 9, 2015 at 1:32
  • $\begingroup$ The function $\delta(x)$ is a Dirac delta function, what expresses $u_{0}$ animals in a defined place from domain, consequently it changes the distribution of animals with the time along the domain, I guess. $\endgroup$ Nov 9, 2015 at 1:45
  • 1
    $\begingroup$ @WaynerKlën : The strict purpose of using the delta function in a boundary (or initial) condition is to correct for units. We can only ever integrate densities, so to turn a number into a density we must multiply it by a function that has vanishing support (place where it is nonzero) but finite integral. The Dirac delta function does this. However, once we abstractualize away from units, we no longer need that device -- we directly specify that the value of $u$, a density, is its initial value and need not use the device of delta to correct units. $\endgroup$ Nov 9, 2015 at 2:58
  • $\begingroup$ Thanks for explanation, this way, you've answered my question fully! Thank you very much for your supporting! $\endgroup$ Nov 9, 2015 at 3:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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