# How could I solve this Reaction-Diffusion PDE using mathematica?

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[{\"(\", RowBox[{\"1\", \",\", \"0\", \",\", \"0\"}], \")\"}], Derivative], MultilineFunction->None][t, x, y] == u[t, x, y] + \ \*SuperscriptBox[\"u\", TagBox[ RowBox[{\"(\", RowBox[{\"0\", \",\", \"0\", \",\", \"2\"}], \")\"}], Derivative], MultilineFunction->None][t, x, y] + \*SuperscriptBox[\"u\", TagBox[ RowBox[{\"(\", RowBox[{\"0\", \",\", \"2\", \",\", \"0\"}], \")\"}], Derivative], MultilineFunction->None][t, x, \ y]$$,True,True,True,True,u[0,x,y]==100}"


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

Thank you for supporting!

• Take a look at NDSolve and NeumannValue. – user484 Nov 4 '15 at 19:00
• @Rahul I'm Using mathematica 9, this way I don't have the NeumannValue feature avaliable to me... – Herr Schrödinger Nov 4 '15 at 19:04
• NeumannValue is convenient but not necessary to solve this system of equations. To obtain more thorough responses, post your equations in Mathematica format. – bbgodfrey Nov 4 '15 at 21:55
• Note that the solution to this equation will grow exponentially at rate r once the initial distribution spreads out over the domain. – Chris K Nov 4 '15 at 22:10
• 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. – obsolesced Jul 18 '16 at 7:41

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?

• 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. – Herr Schrödinger Nov 9 '15 at 1:28
• @WaynerKlën : How is this different from the line u[0,x,y] == 100 (with whatever appropriate constant instead of "100")? – Eric Towers Nov 9 '15 at 1:32
• 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. – Herr Schrödinger Nov 9 '15 at 1:45
• @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. – Eric Towers Nov 9 '15 at 2:58
• Thanks for explanation, this way, you've answered my question fully! Thank you very much for your supporting! – Herr Schrödinger Nov 9 '15 at 3:00