I was having fun modifying a code given to me as an answer to a previous problem here, courtesy of user Alex Trounev (Thank you again), when I encountered a certain error which I had never seen before.

Here is the aforesaid code :

r0 = 0.5;
h = 1;
α = 0.8;

(*region definition*)
reg = Cuboid[{.5, 0., 0.}, {1., 2 Pi, 1.}];

reg3D = ImplicitRegion[
   r0^2 <= x^2 + y^2 <= 1 && 0 <= z <= 1, {x, y, z}];

(*equation + conditions*)
eq1 = D[u[t, r, θ, z], 
    t] - (D[u[t, r, θ, z], r, r] + 
     1/r*D[u[t, r, θ, z], r] - 
     1/(α^2 r^2) D[u[t, r, θ, z], θ, θ] + 
     D[u[t, r, θ, z], z, z]);

ic = u[0, r, θ, z] == 1;

bc = DirichletCondition[u[t, r, θ, z] == Exp[-5 t], r == r0];
nV = NeumannValue[1, r == 1];
pbc = PeriodicBoundaryCondition[u[t, r, θ, z], θ == 0, 
   TranslationTransform[{0, 2 π*α, 0}]];

(*solution computation*)
sol = NDSolveValue[{eq1 == nV, ic, bc, pbc}, 
   u, {t, 0, 2}, {r, θ, z} ∈ reg];

PlotPoints\[Rule]50,PlotLabel\[Rule]Row[{"t = \

When I run the code, after some time, I get greeted with the following error :

NDSolveValue::nlnum: The function value {$Failed} is not a list of numbers with dimensions {39639} at {t,u[t,r,θ,z],(u^(1,0,0,0))[t,r,θ,z]} = {0.0138161,{<<1>>},{-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,<<15>>,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,-4.66626,<<39589>>}}.

When I click on the three dots next to the error, I don't find any information on the error like it's usually the case. I then decide to google some answers. I found some answers here while also trying to comprehend the error by looking at this and finally that answer here.

So if I did understand it correctly, such error arises when you use NDSolve (or NDSolveValue) to get a symbolical solution to your equation, but problems come up when you try to numerically evaluate it for plotting purpose, or when trying to get a symbolical result with a function that requires numerical values ?

In any case, I do not really understand why I get such error as my plot part is currently between (* ... *) so it shouldn't matter. As for the rest of the code, I do not really see an error but I am just a beginner so...

Anyway, can a kind fellow enlighten me please ?

Edit 1 : Yes I forgot to tell you that this is quite the time-consuming computation...sorry.

  • $\begingroup$ I got tired of waiting for it to finish....sorry. $\endgroup$
    – Michael E2
    Jul 5, 2020 at 1:57
  • $\begingroup$ @MichaelE2 I reproduced the error after nearly an hour on my six-processor computer. The computation consumed nearly every available cycle. It produce an InterpolatingFunction over {t, 0, 0.0138}, which appears to go unstable by t = 5 10^-4 $\endgroup$
    – bbgodfrey
    Jul 5, 2020 at 4:03
  • $\begingroup$ There appears to be a sign error in eq1. - 1/(\[Alpha]^2 r^2) D[u[t, r, \[Theta], z], \[Theta], \[Theta]] should be 1/(\[Alpha]^2 r^2) D[u[t, r, \[Theta], z], \[Theta], \[Theta]]. Correct it, and the computation runs correctly. $\endgroup$
    – bbgodfrey
    Jul 5, 2020 at 4:17
  • $\begingroup$ @MichaelE2 It's alright, don't bother, ty nonetheless. $\endgroup$ Jul 5, 2020 at 10:01
  • $\begingroup$ @bbgodfrey Yes, changing the " - " into a "+" fixes the issue...but why does it fix the issue though ? That I'd like to understand. $\endgroup$ Jul 5, 2020 at 10:03

1 Answer 1


For completeness, I summarize my comments here. The computation for eq1, as posted in the question, is violently numerically unstable, due to a sign error. The enormous resulting growth of the solution apparently caused an internal error in NDSolve, and $Failed leaked out. I emphasize, though, that this is a symptom, not the root cause of the failure of the computation. Correcting eq1 to

eq1 = D[u[t, r, θ, z], t] - (D[u[t, r, θ, z], r, r] + 1/r*D[u[t, r, θ, z], r] + 
           1/(α^2 r^2) D[u[t, r, θ, z], θ, θ] + D[u[t, r, θ, z], z, z])

allows the computation to proceed smoothly. To produce the animated plot requested in the question, as opposed to only half of it, replace ArcTan[x, y] by Mod[ArcTan[x, y], 2 Pi] in the final code of the question. (Only every other frame is plotted to reduce the size of the graphic. Alternatives are discussed here.)

enter image description here

  • $\begingroup$ so you do get a half cylinder as well ? That's weird.... $\endgroup$ Jul 5, 2020 at 14:45
  • $\begingroup$ well in any case, thx. $\endgroup$ Jul 5, 2020 at 14:52
  • $\begingroup$ @ConfuzzledStudent I thought you wanted a half cylinder. To get the whole cylinder, use Mod[ArcTan[x, y], 2 Pi]. The problem with your original code is that ArcTan[x, y] returns values in the range {-Pi, Pi}, and the negative values are outside reg and were not computed by NDSolve. I can change my answer accordingly, if you wish. Let me know. $\endgroup$
    – bbgodfrey
    Jul 5, 2020 at 16:58
  • $\begingroup$ @ConfuzzledStudent Alternatively, use reg = Cuboid[{.5, -Pi, 0}, {1., Pi, 1}]. $\endgroup$
    – bbgodfrey
    Jul 5, 2020 at 17:04
  • $\begingroup$ If I was just in the middle of checking the ArcTan documentary, your comment came right on time. If all I need is to replace ArcTan[x, y] with Mod[ArcTan[x, y], 2 Pi] then don't bother editing your answer. The computation does take a lot of time and memory. Thank you for the additional information. $\endgroup$ Jul 5, 2020 at 17:04

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.