# Numerically solving an inhomogeneous partial differential equation

I'm trying to solve a cylindrical partial differential equation with boundary conditions. But I got an error message saying NDSolve::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable. Is there any ways I can solve it?

\begin{align*}\frac1{r}\frac{\partial}{\partial r}\left(r\frac{\partial Y}{\partial r}\right)+\frac{\partial}{\partial z}\frac{\partial Y}{\partial z}-Y&=0\\ r=0,\;Y=1,\;r=1,\;\frac{\partial Y}{\partial r}&=0\\ z=0,\;Y=1,\;z=1,\;\frac{\partial Y}{\partial z}&=0 \end{align*}

And my Mathematica input is :

eqn1 = {1/r*D[y[r, z], r] + D[y[r, z], {r, 2}] + D[y[r, z], {z, 2}] -
y[r, z] == 0, y[1, z] == 1,
y[r, 1] == 1, (D[y[r, z], {r, 1}] /. r -> 0) ==
0, (D[y[r, z], {z, 1}] /. z -> 0) == 0};

NDSolve[eqn1, y, {r, 10^-20, 1}, {z, 10^-20, 1}];


Thank you!

• Don't you have inconsistencies in  y[r, z] == 0, y[1, z] == 1, y[r, 1] == 1 ? – b.gates.you.know.what Nov 6 '12 at 20:02
• @b.gatessucks the y[r,z]==0 is part of the preceding row, or it would be quite trivial :) – ssch Nov 6 '12 at 23:21
• I'm not sure why y[r,z] would be 0... The only way I can think of to solve this question is using Bessel's equation. But it still didn't work and might need to use numerical method... Does anyone know how to solve this with numerical method? Thank you! – DumbleKo Nov 6 '12 at 23:36
• @ssch Oops, the mouse tricked me. Thanks for spotting that. – b.gates.you.know.what Nov 7 '12 at 8:27
• Thank you Mark. This is very helpful!!! I'll let you know if I have more questions. – DumbleKo Nov 7 '12 at 18:16

Edit of July 10, 2014

As of V10, this equation can now be solved with a single, simple call to NDSolve:

y = NDSolveValue[{
r D[y[r, z], z, z] + D[y[r, z], r] + r D[y[r, z], r, r] == r y[r, z],
y[1, z] == 1, y[r, 1] == 1
}, y, {r, 0, 1}, {z, 0, 1}];
ContourPlot[y[r, z], {r, 0, 1}, {z, 0, 1},
ColorFunction -> "TemperatureMap",
ColorFunctionScaling -> False,
Contours -> Range[0, 1, 0.02]] Pre-V10

Below is the technique I suggested pre-V10.

You are trying to solve an elliptic, purely spatial problem. As described in the documentation on numerical PDEs, Mathematica uses the numerical method of lines, which requires a temporal variable with initial (not boundary) conditions. That is the proximate cause of your problem.

Honestly, in this case, you might consider trying another system that implements the finite element method - the preferred numerical method for dealing with elliptic equations. I'm sure that FEM will eventually be in Mathematica, but it is not, as of version 9.

The standard way to solve this type of problem in Mathematica, though, is the method of relaxation. The idea is to realize your solution as the steady state limit of a parabolic equation. There is set up involved, however. First, let's examine your equation

$$y^{(0,2)}(r,z)+\frac{1}{r}y^{(1,0)}(r,z)+y^{(2,0)}(r,z)-y(r,z)=0.$$

Clearly, the singularity at $r=0$ is likely to cause a problem. If we take the limit as $r\rightarrow 0$, however, this second term becomes $y^{(2,0)}(0,z)$. To see this, note that, since $y$ describes a rotationally symmetric function, it extends to negative values of $r$ naturally by $y(-r,z)=y(r,z)$. This extension defines $y$ to be an even function and it's easy to show that the derivative of a differentiable even function is zero at the origin. Thus, $y^{(1,0)}(0,z)=0$ and $$\lim_{r\rightarrow 0} \frac{1}{r}y^{(1,0)}(r,z) = \lim_{r \rightarrow 0} \frac{y^{(1,0)}(r,z) - y^{(1,0)}(0,z)}{r-0} = y^{(2,0)}(0,z).$$

Thus, at $r=0$, we might write the equation as

$$y^{(0,2)}(r,z)+2y^{(2,0)}(r,z)-y(r,z)=0.$$

In Mathematica speak, this could be accomplished using discontinuous helper functions, f1 and f2:

Clear[y];
f1[r_ /; r > 0] := 1;
f1[r_ /; r == 0] := 2;
f2[r_ /; r > 0] := 1/r;
f2[r_ /; r == 0] := 0;
-y[t, r, z] + Derivative[0, 0, 2][y][t, r, z] +
f2[r]*Derivative[0, 1, 0][y][t, r, z] +
f1[r]*Derivative[0, 2, 0][y][t, r, z] ==
Derivative[1, 0, 0][y][t, r, z] Note that I've included the temporal variable t and the first order derivative with respect to t. Experience with parabolic equations, like the heat equation, suggests that the transient part of the solution will decay exponentially to your steady state. Also, note how nicely expressions involving Derivative typeset. This is perfectly good Input, as well, and that typeset version is my preferred method of input.

Now, for the initial condition. Ideally, it should be consistent with your boundary conditions. The simplest such function I could figure is

$$y(0,r,z)=1 - (1 - r^2)(1 - z^2)$$

Using all this, we can write your equations as

eqns = {
-y[t,r,z] + Derivative[0,0,2][y][t,r,z] +
f2[r]*Derivative[0,1,0][y][t,r,z] +
f1[r]*Derivative[0,2,0][y][t,r,z] ==
Derivative[1,0,0][y][t,r,z],
y[t,1,z] == 1, y[t,r,1] == 1,
Derivative[0,1,0][y][t,0,z] == 0,
Derivative[0,0,1][y][t,r,0] == 0,
y[0,r,z] == 1 - (1 - r^2)*(1 - z^2)
};


Now, we solve it.

y[t_,r_,z_] = y[t,r,z] /. First[
NDSolve[eqns, y[t,r,z],{t,0,1.2},{r,0,1},{z,0,1},
Method->{"MethodOfLines", Method->"StiffnessSwitching",
"DifferentiateBoundaryConditions"->{True, "ScaleFactor"->1}}]];


As your equation is expressed in cylindrical coordinates but independent of $\theta$, I guess that, for relatively large fixed $t_0$, $y(t_0,r,z)$ represents a rectangular slice of the cylinder with one edge right down the middle of the cylinder. We can visualize it like so.

ContourPlot[y[1.2, r, z], {r, 0, 1}, {z, 0, 1},
ColorFunction -> "TemperatureMap",
ColorFunctionScaling -> False,
Contours -> Range[0, 1, 0.02]] The "relaxation" of the parabolic solution into the steady state solution can be visualized as an animation.

pics = Table[
Labeled[
ContourPlot[y[t, r, z], {r, 0, 1}, {z, 0, 1},
ColorFunction -> "TemperatureMap",
ColorFunctionScaling -> False,
Contours -> Range[0, 1, 0.02],
ContourStyle -> Opacity[0.5]],
Row[{"t = ", Pane[t, {50, 12}]}]],
{t, 0, 1.2, 0.02}];
ListAnimate[pics]
(* Or
Export["anim.gif", pics]
to export to the following animated GIF.
*) • Er…why does $\frac{1}{r}y^{(1,0)}(r,z)$ become $y^{(2,0)}(r,z)$ when $r\rightarrow 0$? – xzczd Mar 8 '13 at 8:36
• @xzczd By the definition of the derivative! More properly really, it becomes $y^{(1,0)}(0,z)$, as I've indicated in my edit. – Mark McClure Mar 8 '13 at 15:21
• Oh, I see, it's symmetric and smooth at $r=0$, right? Now I can understand it both by considering $y$ as a steady state of something like heat distribution and considering the Cartesian form of the equation. – xzczd Mar 9 '13 at 5:18
• @xzczd Exactly! I guess you can also see this in the images by reflecting them about the vertical axis. – Mark McClure Mar 9 '13 at 8:49
• I find your method inspiring! Do you think if it works for this case ? mathematica.stackexchange.com/questions/77981/… – mastrok Mar 27 '15 at 4:25