I have the following code for a 2d heat c equation:

 soln = NDSolve[{        D[y[r, f, t], r, r] + 
                   (1/r)*D[y[r, f, t], r] +
                 (1/r^2)*D[y[r, f, t], f, f] == 
                                (1/1.27*10^-4)*D[y[r, f, t], t], 
          y[r, f, 0] == 0, 
          y[10^-2, f, t] == 0,
          (D[y[r, f, t], r] /. 
               r -> 2*10^-6) == -(1/(3.14*10^-8))*(1 + Sin[2*2*3.14*t/2]), 
          (D[y[r, f, t], f] /. f -> 0) == -(10/238)*y[r, 0, t], 
          (D[y[r, f, t], f] /. f -> 3.14/2) == 0}, 
                    y, {r, 2*10^-6, 10^-2}, {f, 0, 3.14/2}, {t, 0, 1}]

and I'm getting the two following error messages:

NDSolve::ibcinc: Warning: Boundary and initial conditions are inconsistent

NDSolve::mxst:Maximim number of 10000 steps reached at the point t==1.629*10^-10

NDSolve::eerr:Scaled local spatial error estimate of 4203.0461 at t=1.840*10^-8 in the direction of independent variable r is much greater than prescribed error tolerance. Grid spacing with 51 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or you may want to specify a smaller grid spacing using the MaxStepSize or MinPoints method options

I saw an answer for the first one so i'm pretty ok with that but the second one is giving me trouble limitting the useful time domain.I cannot figure out why I'm getting this. Any help will be greatly appreciated! Claudiu

  • 3
    $\begingroup$ It's possible that the techniques outlined in the answer to this question might help. $\endgroup$ Mar 3, 2014 at 14:13
  • 1
    $\begingroup$ If you set your time domain very small {t,0,10^-9} you will get "a" solution, which is blowing up at r=.1,f=0. Take a look at that and think why your bcs are driving a singularity at the corner. Describing the physical problem in your question might be useful if you need more help. $\endgroup$
    – george2079
    Mar 3, 2014 at 15:34
  • 1
    $\begingroup$ @MarkMcClure notice you mention finite elements in the linked answer - have you seen this: "Method" -> {"PDEDiscretization" -> "FiniteElement"} ? It doesn't work and isn't documented that i can find. $\endgroup$
    – george2079
    Mar 3, 2014 at 19:22
  • 1
    $\begingroup$ @george2079 The new reference site is set up to hold lots of new information of finite element methods. Most of these pages don't yet exist, but presumably will in the near future: reference.wolfram.com/search/?q=fem $\endgroup$ Mar 3, 2014 at 21:55
  • 1
    $\begingroup$ Thanks looking forward to seeing those docs. @user12353 my feeling is you have formulated this correctly and NDSolve is simply not up to the task. $\endgroup$
    – george2079
    Mar 4, 2014 at 16:19

2 Answers 2


Partial answer: the inconsistent boundary condition is that you are specifying df/dr non zero for all t, which is inconsistent with the zero (hence df/dr=0) initial condition.

Ive modified the third condition as:

    (D[y[r, f, t], r] /. 
           r -> 2*10^-6) == -(1/(Pi 10^-8))
                (1 +  Sin[2 Pi t]) Piecewise[{{t 10, t < 1/10  }, {1, True}}]

"ramping up" that BC from zero eliminates the inconsistency warning. However, the solution is still diverging in the corners:

 Plot3D[(y /. soln[[1]])[r, f, 1*^-8] , {r, 2*^-6, 10*^-2}, {f, 0, 
       3.14/2}, PlotRange -> {-100, 1000}]

enter image description here


I'd like to point out OP is probably simulating a point heat source with infinite strength at the center of the circle with the b.c. (D[y[r, f, t], r] /. r -> 2*10^-6) == -(1/(3.14*10^-8))*(1 + Sin[2*2*3.14*t/2]) and this seems to be the most troublesome part of the problem. (To learn more about "point source", check this post. ) I'm not sure about why OP chose this approximation, and guess a better approximation is necessary if one wants to solve the PDE in a not-that-slow way. Anyway, the following is my solution.

First, some clean-up for the equation and parameters, just to make the code more readable:

rrange = {rL, rR} = {2 10^-6, 10^-2};
frange = {fL, fR} = {0, Pi/2};
tbegin = 0;
tend = 1;
eqn = D[y[r, f, t], r, 
     r] + (1/r) D[y[r, f, t], r] + (1/r^2) D[y[r, f, t], f, f] == (100/127 10^-4) D[
     y[r, f, t], t];

ic = y[r, f, tbegin] == 0;
bc = {y[r, f, t] == 0 /. r -> rR, 
   D[y[r, f, t], r] == -((1 + Sin[(2 2 π t)/2])/(π/10^8)) /. r -> rL,
   D[y[r, f, t], f] == -(10/238)*y[r, f, t] /. f -> fL, 
   D[y[r, f, t], f] == 0 /. f -> fR};

Then in principle the PDE can be solved with

showStatus[status_] := 
   SetNotebookStatusLine[FrontEnd`EvaluationNotebook[], ToString[status]]];
clearStatus[] := showStatus[""];

mol[n_Integer, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

soln = NDSolve[{eqn, ic, bc}, y, {r, rL, rR}, {f, fL, fR}, {t, tbegin, tend}, 
    MaxSteps -> Infinity, Method -> Union[mol[True, 100], mol[10, 2]]]; // AbsoluteTiming

Definition of showStatus is from here.

MaxSteps -> Infinity is set to deal with the mxst warning. "DifferentiateBoundaryConditions"->{True,"ScaleFactor"->100} is set to avoid the inconsistent b.c. being severely ignored. (For more information check this post.) "MaxPoints" -> 10, "MinPoints" -> 10, "DifferenceOrder" -> 2 is set to speed up the calculation. According to the monitor, the time step passes 1.840*10^-8 without difficulty, but sadly the calculation is still slow and memory-consuming and my laptop can't bear with it. If you have a PC with more memory, you may have a try, but in the rest part of this answer, I'll decrease the computation burden further by discretizing the PDE to a set of ODE with a coarser grid. (NDSolve won't allow this coarse grid if the PDE is directly sent into NDSolve.)

I'll use the pdetoode function here for discretization.

rpoints = 8;
fpoints = 8;
rgrid =Array[# &, rpoints, rrange];
fgrid = Array[# &, fpoints, frange];
var = Outer[y, rgrid, fgrid];
xdifforder = 2;
torder = 1;
ptoo = pdetoode[y[r, f, t], t, {rgrid, fgrid}, xdifforder];

del = Delete[#, {{1}, {-1}}] &;

odeqn = del /@ del@ptoo@eqn;
odeic = ptoo@ic;
odebc = With[{sf = 100}, diffbc[{t, torder}, sf]@MapAt[del, ptoo@bc, {{1}, {2}}]];

SetAttributes[y, NHoldAll]

sollst = NDSolveValue[{odeqn, odeic, odebc} // N, var, {t, tbegin, tend}, 
    Method -> {"EquationSimplification" -> "Solve"}, MaxSteps -> Infinity, 
    EvaluationMonitor :> showStatus["t = " <> ToString[CForm[t]]]]; // AbsoluteTiming
(* {95.391782, Null} *)

sol = rebuild[sollst, {rgrid, fgrid}, 3];

Animate[RevolutionPlot3D[sol[r, f, t], {r, rL, rR}, {f, fL, fR}, 
  PerformanceGoal -> "Quality", ColorFunctionScaling -> False, Mesh -> None, 
  ColorFunction -> (ColorData["TemperatureMap"][#3/(1.5 10^5)] &), 
  PlotRange -> {Automatic, Automatic, {0, 1.5 10^5}}], {t, 0, tend}]

With such a coarse grid, the result still looks not bad:

enter image description here

  • $\begingroup$ A great answer, I'd like to know what MMA can't do with this pdetoode tool. $\endgroup$
    – user69323
    Jun 3, 2020 at 4:32

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.