# Difficulty in solving a differential equation with custom solver for a different initial condition

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

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.

(*constants*)
a = 0.001; hf = 0.01; \[Mu] = 1; \[Beta] = 1;
(*Equations*)
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 = SimplifyPWToUnitStep[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[FrontEndEvaluationNotebook[],ToString[status]]]; clearStatus[] := showStatus[""]; clearStatus[] 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 (*postprocess*) 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! Update: 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. • Is this model taken from the paper "Nonlocal description of evaporating drops"? May 17, 2023 at 14:29 • @AlexTrounev, yes exactly, but the current method is not exactly the same. – lxy May 18, 2023 at 2:12 • This problem can be solved with using Navier-Stocks and FEM as mentioned in my answer on physics.stackexchange.com/questions/509623/… May 18, 2023 at 3:45 • @Alex Troune, it is a promising approach but may require more computation cost than the model. Please elaborate. Thank you in advance! – lxy May 18, 2023 at 5:03 • 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. May 18, 2023 at 9:52

If "the IC only changed in its amplitude", changing sf = 1 to sf = 10 resolves the problem.
• @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?
• @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. May 20, 2023 at 4:33