This is a follow-up question to this previous post by @FLP, in which an interesting system of equations was solved with the useful pdetoode developed by @xzczd. I have tried to solve this problem with a different initial condition (IC). And I appreciated the help and effort from @xzczd and @FLP. The following code is a combination of the original post and the answer. Due to the symmetry, the differential equation for h[r,t] is to be solved from $r=0$ to a sufficiently large value $r>>1$. Here, the left bound of $r$ is taken as a small non-zero value due to some numerical reason, and the right bound is taken as $50$.

My question:

In the original post, a system about $h[r,t]$ was solved with an IC of the exponential function IC = {h[r, t] == 20 E^-(r/20)^2 + hf} /. t -> 0;. For some weird reason if you change the amplitude of the IC, say, IC = {h[r, t] == 10 E^-(r/10)^2 + hf} /. t -> 0;, it would stuck at a small time (relative to a given Tend) and return an error: step size is effectively zero; singularity or stiff system suspected.

I also tried a parabolic IC defined by If:

IC = {h[r, t] == unitStepExpand@If[r >= lb && r < 1, 3/5*(1 - r^2) + hf, hf]} /. t -> 0

but the error of a singularity remained. Noted that in this case PiecewiseExpand transforms If[...] to an equivalent Piecewise[...]; then SimplifyPWToUnitStep transforms Piecewise[...] to an equivalent combination of UnitStep[...]. This is a necessary transformation since pdetoode assumes the functions in the equations, ICs, and BCs are all Listable.

For the piecewise IC, I guess the issue may be caused by the discontinuous derivative at $r=1$. But I cannot understand why the code fails with the IC only changed in its amplitude: h[r, t] == 10 E^-(r/10)^2 + hf.

An even more strange issue is that if I change the side of the domain with, say rb=100, a singularity occurs again.

The code:

The code for pdetoode neglected here, please find in the abovementioned link. Here, Int[r,t] is an integral function, defined as $$Int(r,t)=\int_0^\infty K(r, r^\prime)\frac{\partial}{\partial r^\prime }\left(\frac{h_f}{h} \right)^3 \mathrm{d}r^\prime, $$ where $K(r,r^\prime)=\frac{2}{\pi}r[\mathbf{K}(r^\prime/r)-\mathbf{E}(r^\prime/r)]$ for $r^\prime<r$

or $K(r,r^\prime)=\frac{2}{\pi}r^\prime[\mathbf{K}(r/r^\prime)-\mathbf{E}(r/r^\prime)]$ for $r^\prime>r$, in which $\mathbf{K}$ and $\mathbf{E}$ are the complete elliptic integrals.

In order to avoid the singularity of $K(r,r^\prime)$ at $r=r^\prime$, $r$ takes the values at the grid point lb+i*step and $r^\prime$ takes the values at the staggered grid lb+(i+0.5)*step.

a = 0.001; hf = 0.01; \[Mu] = 1; \[Beta] = 1;
With[{h = h[r, t], P = P[r, t], Q = Q[r, t], J = J[r, t]},
EqnP = P == D[h, r, r] + D[h, r]/r + a^2/h^3;
EqnQ = Q == 1/3*h^3*r*D[P, r];
EqnC = D[h, t] + (\[Mu] D[Q, r])/r == -J;
EqnJ = J == (\[Beta] D[Int[r, t], r])/r;]

(*Computation Domain*)
lb = 1/100; rb = 50;

(*unitStepExpand is a function that transforms If[\[Ellipsis]] to mathematically equivalent combination of UnitStep[\[Ellipsis]], useful in the computation of the following integral*)
unitStepExpand = Simplify`PWToUnitStep[PiecewiseExpand[#]] &;

(*initial and boundary conditions*)
IC = {h[r, t] == 20 E^-(r/20)^2 + hf} /. t -> 0;
(*IC = {h[r, t] == 10 E^-(r/10)^2 + hf} /. t -> 0;*)
(*IC = {h[r, t] == unitStepExpand@If[r >= lb && r < 1, 3/5*(1 - r^2) + hf, hf]} /. t -> 0;*)
BC = {{D[h[r, t], r] == 0, D[h[r, t], r, r, r] == 0, D[h[x, t], x] == 0, D[h[x, t], x, x, x] == 0} /. {r -> lb, x -> rb}};

(*Discretization setting*)
difforder = 2;
points = 200;
grid = Array[# &, points, {lb, rb}];

hgrid = Map[h[#][t] &, grid];
hdgrid = fdd[Derivative[1], grid, hgrid];

(*the Integral Function*)
Kernel[r_, i_] := (2*r)/\[Pi]*unitStepExpand@If[r > i, r (EllipticK[i/r] - EllipticE[i/r]), i (EllipticK[r/i] - EllipticE[r/i])];

Integrand = D[(hf/h[r, t])^3, r] /.Thread[{D[h[r, t], r], h[r, t]} -> {hdgrid, hgrid}];
step = (rb - lb)/points;
Int[r_][t_] := Total@Table[step*Kernel[r, lb + (i + 0.5)*step]*0.5 (Integrand[[i]] + Integrand[[i + 1]]), {i, 1, points - 1}];

(*Discretising the rest of the system*)
ptoofunc = pdetoode[{h[r, t], P[r, t], Q[r, t], J[r, t], Int[r, t]}, t, grid, difforder];
removeredundant = #[[3 ;; -3]] &;

odeadd = Map[ptoofunc, {EqnP, EqnQ, EqnJ}, {2}];
ode = Block[{P, Q, J}, Set @@@ odeadd; ptoofunc@EqnC] // removeredundant;
odIC = ptoofunc@IC;

With[{sf = 1}, odBC = diffbc[t, sf]@BC // ptoofunc];

(*Solving the System*)
var = Map[h, grid];

showStatus[status_] := LinkWrite[$ParentLink,SetNotebookStatusLine[FrontEnd`EvaluationNotebook[],ToString[status]]];
clearStatus[] := showStatus[""];
monitor[t_] := EvaluationMonitor :> showStatus["t = " <> ToString[CForm[t]]]

Tend = 200;

sollst = NDSolveValue[{ode, odIC, odBC}, var, {t, 0, Tend}, monitor[t], Method -> {EquationSimplification -> MassMatrix}]; // AbsoluteTiming

solfunc = rebuild[sollst, grid];

Manipulate[Plot[solfunc[t, r], {r, lb, rb}, PlotRange ->All], {t, 0, Tend}]

Some notes:

  1. This code works well with the original IC, but when you change the amplitude of the IC, say, from $20$ to $10$, a singularity will be encountered. Increasing the point number and order of FD has little improvement.

  2. I am wondering if a nonuniform mesh (or adaptive mesh) could be implemented in the framework of pdetoode such that a refined resolution is bunched towards the region where the change of slope of h[r,t] is large based on physical reason?

  3. As mentioned by @FLP, the original code has some additional problems as well. I'm not sure if the treatment of IC will be one of the potential issues. So I post this question here to draw attention.

Please give me some suggestions. Thank you in advance!


  1. With the suggestion of xzczd, the code can run with IC = {h[r, t] == 10 E^-(r/20)^2 + hf} /. t -> 0;, which only changes in its amplitude, up to T=502 approximately, and then singularity occurs at the radius position of the blob, where the slope changes most rapidly.

  2. By changing the width of the initial blob, I found smaller the initial width earlier the singularity occurs. For example, using IC = {h[r, t] == unitStepExpand@If[r >= lb && r < 1, 3/5*(1 - r^2) + hf, hf]} /. t -> 0, the singularity occurs at t=0.000024 approximately. Increasing points improves things a little.

  • $\begingroup$ Is this model taken from the paper "Nonlocal description of evaporating drops"? $\endgroup$ May 17, 2023 at 14:29
  • $\begingroup$ @AlexTrounev, yes exactly, but the current method is not exactly the same. $\endgroup$
    – lxy
    May 18, 2023 at 2:12
  • $\begingroup$ This problem can be solved with using Navier-Stocks and FEM as mentioned in my answer on physics.stackexchange.com/questions/509623/… $\endgroup$ May 18, 2023 at 3:45
  • $\begingroup$ @Alex Troune, it is a promising approach but may require more computation cost than the model. Please elaborate. Thank you in advance! $\endgroup$
    – lxy
    May 18, 2023 at 5:03
  • $\begingroup$ Actually it could be much faster to solve linear equation $\nabla^2 c=0$ with appropriate boundary conditions on the moving border and at infinity then nonlinear integrodifferential equation. $\endgroup$ May 18, 2023 at 9:52

1 Answer 1


If "the IC only changed in its amplitude", changing sf = 1 to sf = 10 resolves the problem.

  • $\begingroup$ Is it answer or comment? $\endgroup$ May 20, 2023 at 3:10
  • $\begingroup$ @xzczd, good to see you. I tried sf=10 with i.c. IC = {h[r, t] == 10*E^-(r/10)^2 + hf} /. t -> 0 or IC = {h[r, t] == unitStepExpand@If[r >= lb && r < 1, 3/5*(1 - r^2) + hf, hf]} /. t -> 0. The singularity issue remains in each case. With a larger sf, are you trying to resolve a problem of inconsistent i.c. and b.c.s? $\endgroup$
    – lxy
    May 20, 2023 at 3:52
  • $\begingroup$ @xzczd, what did you mean by period in the answer? $\endgroup$
    – lxy
    May 20, 2023 at 3:53
  • $\begingroup$ @Alex An answer of course, if "the IC only changed in its amplitude". $\endgroup$
    – xzczd
    May 20, 2023 at 4:28
  • $\begingroup$ @lxy Again, you've made careless mistake when writing the question, 10*E^-(r/10)^2 + hf is not an "IC only changed in its amplitude", period. $\endgroup$
    – xzczd
    May 20, 2023 at 4:33

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.