3

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 :

(*parameters*)
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];


(frames=Table[DensityPlot3D[sol[t,Sqrt[x^2+y^2],ArcTan[x,y],z],{x,y,

z}∈reg3D,ColorFunction[Rule]"Rainbow",OpacityFunction[Rule]

None,Boxed[Rule]False,Axes[Rule]False,PlotRange[Rule]{0,1.5},

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

",t}],ColorFunctionScaling[Rule]False],{t,.05,1,.05}]
ListAnimate[frames])

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.

AsukaMinato
  • 9,758
  • 1
  • 14
  • 40
ConfuzzledStudent
  • 105
  • 9
  • I got tired of waiting for it to finish....sorry. – Michael E2 Jul 05 '20 at 01:57
  • @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 – bbgodfrey Jul 05 '20 at 04:03
  • 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. – bbgodfrey Jul 05 '20 at 04:17
  • @MichaelE2 It's alright, don't bother, ty nonetheless. – ConfuzzledStudent Jul 05 '20 at 10:01
  • @bbgodfrey Yes, changing the " - " into a "+" fixes the issue...but why does it fix the issue though ? That I'd like to understand. – ConfuzzledStudent Jul 05 '20 at 10:03
  • Regarding the sign error. Diffusion coefficients should be intrinsically positive. What is the physical meaning of a negative diffusion coefficient? – Tim Laska Jul 05 '20 at 11:42
  • @TimLaska might just be an error on my end, I need to thoroughly check my calculations on my notepad. – ConfuzzledStudent Jul 05 '20 at 12:34
  • @bbgodfrey since you managed to run the whole thing, does it show a half cylinder to you as well ? – ConfuzzledStudent Jul 05 '20 at 12:34
  • The sign error in eq1 caused the computation to be violently unstable., which apparently led to an internal failure in NDSolve. Fixing the error allowed the computation to run smoothly, producing the desired animation. – bbgodfrey Jul 05 '20 at 13:24
  • @MichaelE2 Although the ultimate cause of the error cited in the question was violent growth of the solution due to an error in the PDE, it seems to me that $Failed should have been trapped inside of NDSolve instead of being allowed to leak out. I would consider this a (minor) bug. What do you think? – bbgodfrey Jul 05 '20 at 13:27
  • @bbgodfrey I guess it's a bug. I haven't looked into it, but $Failed normally is caught – Michael E2 Jul 05 '20 at 20:20

1 Answers1

6

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

bbgodfrey
  • 61,439
  • 17
  • 89
  • 156
  • so you do get a half cylinder as well ? That's weird.... – ConfuzzledStudent Jul 05 '20 at 14:45
  • well in any case, thx. – ConfuzzledStudent Jul 05 '20 at 14:52
  • @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. – bbgodfrey Jul 05 '20 at 16:58
  • @ConfuzzledStudent Alternatively, use reg = Cuboid[{.5, -Pi, 0}, {1., Pi, 1}]. – bbgodfrey Jul 05 '20 at 17:04
  • 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. – ConfuzzledStudent Jul 05 '20 at 17:04
  • I'll try both methods and will get back to you after the computation ends (in 30 mins ) – ConfuzzledStudent Jul 05 '20 at 17:07
  • It works great, thank you very much. – ConfuzzledStudent Jul 05 '20 at 18:15