0

I am currently trying to solve a transient cylindrical diffusion problem with space variable (r,z). Using NDSolve, I get a certain Interpolation function but I can't plot it no matter what I try.

Here is my code :


ClearAll["Global`*"]
r0 = 0.5;
h = 1;

eq1 = D[u[t, r, z], 
    t] - (D[u[t, r, z], r, r] + 1/r*D[u[t, r, z], r] + 
     D[u[t, r, z], z, z]);

ic = {u[0, r, z] == 0};

bc = {u[t, r0, z] == 0, 
   u[t, 1, z] == 0, (D[u[t, r, z], r] /. r -> r0) == 
    0, (D[u[t, r, z], r] /. r -> 1) == 1, u[t, r, 0] == u[t, r, h]};


sol = NDSolve[{eq1 == 0, ic, bc}, 
  u[t, r, z], {t, 0, 10}, {r, r0, 1}, {z, 0, h}, 
  MaxSteps -> Infinity , MaxStepFraction -> 1/10]

f[t_, r_, z_] := u[t, r, z] /. First[sol];

Manipulate[Plot3D[f[t, r, z], {t, 0, 10}, {r, r0, 1}], {z, 0, 1}]

I do end up with the error NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent. but I think it arises from the approximations of NDSolve, so I dismissed it.

I tried looking for a solution online. I noticed that most blank plots like mine end up being the result of what they call "scoping" errors. From what I gathered, for instance from this, a solution can be found if you correctly define the function that you wish to plot, so that the variables don't end up being local variables in Manipulate.

Therefore, I have tried to define a function from the output of NDSolve using this question as an example.

But even with my definition of the function f, it appears that nothing is being plotted. What am I doing wrong ? Can someone please enlighten me ?

ConfuzzledStudent
  • 105
  • 9

1 Answers1

1

This fixes the plotting issue, but you will have to think about boundary and initial conditions and make them consistent.

ClearAll["Global`*"]
r0 = 0.5;
h = 1;

eq1 = D[u[t, r, z], 
    t] - (D[u[t, r, z], r, r] + 1/r*D[u[t, r, z], r] + 
     D[u[t, r, z], z, z]);

ic = {u[0, r, z] == 0};

bc = {u[t, r0, z] == 0, 
   u[t, 1, z] == 0, (D[u[t, r, z], r] /. r -> r0) == 
    0, (D[u[t, r, z], r] /. r -> 1) == 1, u[t, r, 0] == u[t, r, h]};

sol = NDSolveValue[{eq1 == 0, ic, bc}, 
   u, {t, 0, 10}, {r, r0, 1}, {z, 0, h}, MaxSteps -> Infinity, 
   MaxStepFraction -> 1/10];
Manipulate[Plot3D[sol[t, r, z], {t, 0, 10}, {r, r0, 1}], {z, 0, 1}]
user21
  • 39,710
  • 8
  • 110
  • 167
  • Wow, the answer seems so simple. So if i want to use manipulate and plot the result it would be in my best interest to use NDSolveValue instead of NDSolve, and directly use sol in Plot3D ? In any case, thank you for the answer. – ConfuzzledStudent Jun 15 '20 at 14:37
  • 1
    You get the same result with NDSolve, if you use Set (=) instead of SetDelayed (:=) in definition of f f[t_, r_, z_] = u[t, r, z] /. First[sol]; . – Akku14 Jun 15 '20 at 14:43
  • @Akku14 it does indeed also fix the problem. I was not aware of the problems that could arise from using SetDelayed instead of Set. Thanks a lot. – ConfuzzledStudent Jun 15 '20 at 14:55