# NDSolve refuses to initialise when solving an integro-differential equation with custom pdetoode solver scheme

Problem Background

Recently I'm attempting to replicate the result of the following research paper on the Nonlocal description of Evaporating Drops. The equation of motion of a evaporating, spreading drop is described as

Where the kernel $$K(r,r')$$ is defined as such,

My current approach

To deal with the coupled system of integro-differential equations, I've decided to discretise the system with the useful pdetoode into a massive system of coupled ODEs. By expressing the integral as a sum across spatial coordinates, I can simplify the problem to a coupled ODE problem.

As mentioned by this post, the high order of the differentiation leads to problem with conservation, necessitating numerical calculation of the respective intermediates. The equations are therefore defined as such:

(*Defining Variables*)

a = 0.001;(*pseudo hanmaker's constant*)
hf = 0.01;(*precursor film thickness*)
μ = 1;(*viscosity*)
β = 1;(*Evaporation Constant*)

(*Defining 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] + (μ D[Q, r])/r == -J;
EqnJ = J == (β D[Int[r, t], r])/r;]

(*Defining Initial and Boundary Conditions*)

IC = {h[r, t] == 20 E^-(r/20)^2 + 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}}


The integral, $$\int_{0}^{\infty} K(r,r')\frac{\partial}{\partial r'} (\frac{h}{hf})^3 \,dr'$$ is completed as such:

1. $$\frac{\partial}{\partial r'} (\frac{h}{hf})^3$$ is discretised into a list of differences of $$h(r,t)$$, and the integrand is integrated via the midpoint rule;
2. For $$r'\approx r$$, there is an integrable singularity for $$K(r,r')$$. We approximate $$K(r,r')$$ as $$\frac{r}{\pi}\ln|r'-r|$$ and integrate it symbolically to find the integral for this region.
(*Discretising the Integral Function*)

lb = 1/100; rb = 50;
difforder = 2;
points = 100;

unitStepExpand = SimplifyPWToUnitStep@PiecewiseExpand@# &;
step = (rb - lb)/points;
grid = Array[# &, points, {lb, rb}];

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

Kernel[r_, i_] := (
2 r)/π unitStepExpand@
If[r > i, r (EllipticK[i/r] - EllipticE[i/r]),
i (EllipticK[r/i] - EllipticE[r/i])];

IntegratedSingularity[r_, i_] :=
r/π Sqrt[(r - i)^2] (Log[Sqrt[(r - i)^2]] - 1)

Integrand =
D[(hf/h[r, t])^3, r] /. Thread[{D[h[r, t], r], h[r, t]} -> {hdgrid, hgrid}];

Int[r_][t_] :=
Total@Table[
unitStepExpand@
If[Sqrt[(r - lb - (i + 0.5) step)^2] <=
3 step, (IntegratedSingularity[r, lb + (i + 1) step] -
IntegratedSingularity[r, lb + i step])*0.5 (Integrand[[i]] +
Integrand[[i + 1]]),
step*Kernel[r,
lb + (i + 0.5) step]*0.5 (Integrand[[i]] +
Integrand[[i + 1]])], {i, 1, points - 1}];


The rest of the system is then discretised as such:

(*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;

(*For some weird reason if you apply ptoofunc to a boundary condition at a finite value ie /
D[h[r,t],r]/.r->50==0,it would discretise all values of D[h[r,t],r] prior to 50 (ie /
49,48,47,46...) as 0 as well. I have to manually choose the correct boundary conditions from /
the given list*)

(*Edit: Resolved with fractional lb instead of MachinePrecision*)

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

(*Solving the System*)
var = h /@ grid;
Monitor[sollst =
NDSolveValue[{ode, odIC, odBC}, var, {t, 0, 1},
EvaluationMonitor :> (time = t)], time]


Current Situation

At the NDSolveValue step, however, the code refuses to initialise and the "time" under the EvaluatioMonitor does not change into any numbers. Leaving NDSolveValue to run on its own eventually causes the computer to lag severely and eventually crash. In fact, NDSolveValue does not stop evaluating even if it encountered an error, and refuses to abort after running for a certain amount of time.

Personally, I'm not sure whether this is an issue with computational power, or that there are certain parts of my approach that is erroneous but I am personally unable to see. I would therefore like to ask the following questions:

1. Is my approach to this problem computationally sound?
2. Is there any severe inefficiency in my approach that I can simplify for the problem?
3. Are there any other potential issues with my code?

Thank everyone in advance for any help you could offer!

Edits

At the advice of @xzczd, Ive looked through the integration function again and realised that the integral would take on indeterminate values whenever $$r=r'$$. This is because

$$\int_{a}^{b}\frac{r}{\pi}\ln|r'-r|\,dr'=\frac{r}{\pi}(r-b)(\ln|r-b|-1)-\frac{r}{\pi}(r-a)(\ln|r-a|-1)$$

Notice that whenever $$a=r$$ or $$b=r$$, $$\ln|a-r|=-\infty$$ resulting in the indefinite value as given.

In fact in our case, we don't even have to consider the case of the integrable singularity at all. Since we are taking the trapezoidal rule, $$r$$ can only take on the values lb+step*points while $$r'$$ can only take on the values lb+(step+0.5)*points, so the singularity does not occur in our entire integration. Using purely the trapezoidal rule to integrate $$K(r,r')$$, and comparing the results to NIntegrate:

The fit is honestly not bad.

The integral function discretisation is now redefined as such:

(*Discretising the Integral Function*)

lb = 0.01; rb = 50;
difforder = 2;
points = 100;

unitStepExpand = SimplifyPWToUnitStep@PiecewiseExpand@# &;
step = (rb - lb)/points;
grid = Array[# &, points, {lb, rb}];

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

Kernel[r_, i_] := (
2 r)/π 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}];

Int[r_][t_] :=
Total@Table[
step*Kernel[r,
lb + (i + 0.5) step]*0.5 (Integrand[[i]] +
Integrand[[i + 1]]), {i, 1, points - 1}];


However, the computational problem persists. When the full solver is ran, the computer begins lagging severely and eventually crashes completely by switching off.

A cursory check via Task Manager indicates that the memory usage of NDSolve reaches nearly 98% to 99% of all available memory prior to the crash, potentially explaining the lagging behaviour of the computer.

In that case, is there any method for reducing the memory usage of NDSolve when computing the expression? Perhaps is it possible to conduct numerical preprocessing of the equations before plugging them into NDSolve, or are there any other possible methods?

• 1. It's not a good idea to create sock puppet in this site. 2. "For some weird reason if you apply ptoofunc to a boundary condition at a finite value ie / D[h[r,t],r]/.r->50==0,it would discretise all values of D[h[r,t],r] prior to 50 (ie / 49,48,47,46...) as 0 as well." It's because you've used MachinePrecision number to define lb, which causes error accumulation. Defining lb as e.g. 1/100 or 0.0116 will fix the problem. Dec 17, 2022 at 3:31
• 1. Apologies! Honestly didn't know that it is bad to post under multiple accounts, will note that in the future! (However, since this question is already posted I'll just leave this here to avoid confusion) 2. Thanks for the advice - I'll try that out!
– FLP
Dec 17, 2022 at 4:02
• You can merge the accounts: mathematica.stackexchange.com/help/merging-accounts 3. The Indeterminate generated by Int is causing problem in NDSolveFiniteDifferenceDerivative. (Looks like a bug to me. ) You need to properly modify Int to avoid generating Indeterminate as output of Int to circumvent the issue. Dec 17, 2022 at 4:10
• 4. There's a simple mistake in definition of EqnC. 5. Since you've used diffbc, you should not remove anything from the discretized i.c.. Dec 17, 2022 at 4:32
• 4.,5. Thanks and edited!
– FLP
Dec 17, 2022 at 4:36

## 1 Answer

Aside from those simple mistakes already corrected in the question, there're 2 issues worth elaborating a bit, so let me write an answer.

## 1

For some weird reason if you apply ptoofunc to a boundary condition at a finite value ie D[h[r, t] ,r] /. r->50 == 0,it would discretise all values of D[h[r,t],r] prior to 50 (ie 49, 48, 47, 46...) as 0 as well. I have to manually choose the correct boundary conditions from the given list.

Edit: Resolved with fractional lb instead of MachinePrecision.

This is because pdetoode internally uses PossibleZeroQ to detect whether an equation lies on the boundary of the computational domain, but:

With[{lb = 0.01, rb = 50, points = 100},
Array[# &, points, {lb, rb}][[-1]] - rb]
(* -7.10543*10^-15 *)
% // PossibleZeroQ
(* False *)


As we can see, the numeric error caused by MachinePrecision number is so large that PossibleZeroQ doesn't think the 2 numbers equal. Actually this is a common problem when using MachinePrecision number for indexing, you may check the following posts for more info:

Strange behavior when defining a symbol/function pointwise

Strange behavior of MemberQ, Position

Position function not always returning an answer even with no apparent problems

Can Someone Please Explain Internal$SameQTolerance? So, by changing definition of lb to lb = 1/100 or lb = 0.0116, and adding Method -> {EquationSimplification -> Solve} or Method -> {EquationSimplification -> Residual} or SolveDelayed -> True or Method -> {EquationSimplification -> MassMatrix} to NDSolveValue, the equation system will be solved without difficulty: showStatus[status_]:=LinkWrite[$ParentLink,
SetNotebookStatusLine[
FrontEndEvaluationNotebook[],ToString[status]]];
clearStatus[]:=showStatus[""];
clearStatus[]
jianshi[t_]:=EvaluationMonitor:>showStatus["t = "<>ToString[CForm[t]]]

tend = 200;

sollst =
NDSolveValue[{ode, odIC, odBC}, var, {t, 0, tend}, jianshi[t],
Method -> {EquationSimplification -> MassMatrix}]; // AbsoluteTiming
(* {4.63116, Null} *)

solfunc = rebuild[sollst, grid];

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


## 2

The other issue, though already circumvented by the code in Edits, deserves to be mentioned in an answer so it can be easily searched because it's caused by a surprsing and long-standing bug of NDSolveFiniteDifferenceDerivative.

As mentioned above, the original Int sometimes give Indeterminate as the output, but:

grid = Range@6/10;
values = 0.1 ReplacePart[h /@ grid, {1} -> Indeterminate];
NDSolveFiniteDifferenceDerivative[1, grid, values]


As we can see, when Indeterminate involves in, NDSolveFiniteDifferenceDerivative arbitrarily changes the precision of values, including the those inside h. What's more surprising:

grid = Range@6/10;
values = 0.1 ReplacePart[h /@ grid, {1} -> Indeterminate];
code := NDSolveFiniteDifferenceDerivative[0, grid, values]
code
(* Start a new cell here. *)
code


Cell division influences the output! These behaviors definitely look like bugs. Already reported to WRI.

• already read this link you suggested, I tried a piecewice IC, IC= {h[r, t] == If[r >= lb && r <= 1, 3/5 (1 - r^2) + hf, hf]} /. t -> 0, with an additional modification odIC = ptoofunc@unitStepExpand@IC; NDSolveValue gives an error <The number of constraints (198) (ics) is not equal to the total differential order of the system plus the number of discrete variables (100)>. (to continue)
– lxy
Mar 20 at 14:42
• @jsxs Check carefully about how I use unitStepExpand. If you still feel confused, observe the output of PiecewiseExpand@IC and think about what trouble will arise then. Mar 20 at 15:11
• Oops, I should have extracted the i.c. in the form of an If function with IC[[1, 2]]. I have modified the error, but when executing odIC = ptoofunc@unitStepExpand@IC[[1, 2]], i got another error >"Part 1 of {} does not exist". I have also checked unitStepExpand@IC[[1, 2]], which looks correct. Thank you for the help!
– lxy
Mar 21 at 1:26
• @jsxs …It should be IC= {h[r, t] == unitStepExpand@If[r >= lb && r <= 1, 3/5 (1 - r^2) + hf, hf]} /. t -> 0`. Mar 21 at 2:56
• @jsxs by the way if you are still interested in this simulation I think the original code has some additional problems as well that was detected later on... Didn't really include it into the edits (the post is already long enough) but may post another answer in the near future to include all necessary edits.
– FLP
Apr 6 at 7:13