Numerically solving an inhomogeneous partial differential equation

Question

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!

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

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.
*)