I'm in the midst of solving the given problem

Given the setup I tried writing up the following:

β = 0.05;
sol = NDSolveValue[{
        D[Subscript[U, z][η, τ], τ] == -Sin[τ] + β (1/η D[η*Subscript[U, z][η, τ], η, η]),
        Subscript[U, z][η, 0] == 0, Subscript[U, z][1, τ] == 0,},
        Subscript[U, z], {τ, 0, 10}]

How can I model the second boundary condition accordingly, and what should I do next to model the solution on a 3D plot?

  • 3
    $\begingroup$ Have a look at Menu/Help/WolframDocumentation/NeumannValue and there go to Applications/Time-Dependent Problems. There you will find examples which are no exactly equal to yours (that is, not radial), but very close and with comparable type of Neumann and initial conditions. $\endgroup$ Nov 19, 2019 at 8:48
  • $\begingroup$ The book does not reference that function so I will not use that $\endgroup$
    – TexMexDex
    Nov 19, 2019 at 21:22

1 Answer 1


A couple things. Since $\eta$ is in the denominator of your pde, your solution will have a problem when $\eta$ is 0. Luckily since zero is one of the end points we can pick a number close to zero without changing the answer much. A numerical solution is an approximation after all. Also your NDSolveValue needs the limits of $\eta$ as well as $\tau$. Fixing that and choosing a very small number for zero, we have:

β = 0.05;
zero = 10^-12;

sol = NDSolveValue[{D[U[η, τ], τ] == -Sin[τ] + β*(1/η)*D[η*D[U[η, τ], η], η], 
    U[η, 0] == 0, U[1, τ] == 0, Derivative[1, 0][U][zero, τ] == 0}, U, {η, zero, 1}, 
   {τ, 0, 10}]

tp = Table[
   Plot[Evaluate[sol[η, τ]], {η, zero, 1}, 
    PlotRange -> {-2, 1}], {τ, 0, 10, .1}];

enter image description here

In this particular case we can include $\eta = 0$ in the solution because the derivative wrt to $\eta$ is zero there. Expanding the pde we have

Derivative[0, 1][U][η, τ] == (β*Derivative[1, 0][U][η, τ])/η + 
   β*Derivative[2, 0][U][η, τ] - Sin[τ]

the term with $\eta$ in the denominator is $0/0$ at $\eta = 0$, so we can use L'Hospital's rule. Rewriting the pde with Piecewise functions to separate the zero case from the general case and solving we get:


pde = Derivative[0, 1][U][η, τ] == 
   Piecewise[{{β*Derivative[2, 0][U][η, τ], η == 0}, 
      {(β*Derivative[1, 0][U][η, τ])/η, True}}] + β*Derivative[2, 0][U][η, τ] - 

β = .05

sol = NDSolveValue[{pde, U[η, 0] == 0, U[1, τ] == 0, 
   Derivative[1, 0][U][0, τ] == 0}, 
  U, {η, 0, 1}, {τ, 0, 10}];

The resulting plot looks the same as before, but we can now start the solution at $\eta = 0$.

  • $\begingroup$ What is the purpose of the code: β*Derivative[2, 0][U][η, τ] $\endgroup$
    – TexMexDex
    Nov 22, 2019 at 0:21
  • $\begingroup$ It comes from L'Hospital's rule. When you have f[x]/g[x] = 0/0 at x = 0. You can often find the limit as x->0 by finding f'[x]/g'[x]. So when you have in this case an expression like u'[x]/x, find the limit at 0 by taking the numerator and denominator derivatives and get u''[x]/1. $\endgroup$
    – Bill Watts
    Nov 22, 2019 at 1:24

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.