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!

  • $\begingroup$ Don't you have inconsistencies in ` y[r, z] == 0, y[1, z] == 1, y[r, 1] == 1` ? $\endgroup$ Nov 6, 2012 at 20:02
  • $\begingroup$ @b.gatessucks the y[r,z]==0 is part of the preceding row, or it would be quite trivial :) $\endgroup$
    – ssch
    Nov 6, 2012 at 23:21
  • $\begingroup$ 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! $\endgroup$
    – DumbleKo
    Nov 6, 2012 at 23:36
  • $\begingroup$ @ssch Oops, the mouse tricked me. Thanks for spotting that. $\endgroup$ Nov 7, 2012 at 8:27
  • $\begingroup$ Thank you Mark. This is very helpful!!! I'll let you know if I have more questions. $\endgroup$
    – DumbleKo
    Nov 7, 2012 at 18:16

1 Answer 1


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

enter image description here


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


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


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

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]

enter image description here

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

enter image description here

The "relaxation" of the parabolic solution into the steady state solution can be visualized as an animation.

pics = Table[
    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}];
(* Or
     Export["anim.gif", pics]
   to export to the following animated GIF.

enter image description here

  • $\begingroup$ Er…why does $\frac{1}{r}y^{(1,0)}(r,z)$ become $y^{(2,0)}(r,z)$ when $r\rightarrow 0$? $\endgroup$
    – xzczd
    Mar 8, 2013 at 8:36
  • $\begingroup$ @xzczd By the definition of the derivative! More properly really, it becomes $y^{(1,0)}(0,z)$, as I've indicated in my edit. $\endgroup$ Mar 8, 2013 at 15:21
  • $\begingroup$ 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. $\endgroup$
    – xzczd
    Mar 9, 2013 at 5:18
  • $\begingroup$ @xzczd Exactly! I guess you can also see this in the images by reflecting them about the vertical axis. $\endgroup$ Mar 9, 2013 at 8:49
  • $\begingroup$ I find your method inspiring! Do you think if it works for this case ? mathematica.stackexchange.com/questions/77981/… $\endgroup$
    – mastrok
    Mar 27, 2015 at 4:25

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.