# Fourth order, nonlinear PDE depicting dynamics of liquid lenses

PROBLEM STATEMENT

Recently, I was trying to verify the solution of a set of thin film spreading equations given by this paper.

Where $$p1=-(h_{1xx}+\sigma h_{2xx})$$

$$p2=-\sigma h_{2xx}-\Pi[h2-h1]$$

$$\Pi[h]=\frac{k}{h_f}[(\frac{h_f}{h})^4-(\frac{h_f}{h})^3]$$

However, given the high order and non-linear nature of the equations, Mathematica keeps returning me enormous spatial errors even under extremely fine meshes of 1000 points and a difforder of 5 in the finite difference scheme.

NDSolveValue::eerr: Warning: scaled local spatial error estimate of 68548.44786079538 at t = 0.03579162478117444 in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 1000 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.


The enormous error, in turn, leads to non-physical behaviours in the solution field, which hinders me from solving it under reasonable time-scales.

While increasing both the mesh size and the difference order modestly reduce the error, both step increase the solve time significantly. (it now takes more than 2 hours to solve for 0.1s of spreading time for 5000 points and a difference order of 11) Since I am planning to solve an even more complex version of this equation, is there any way to reduce the error of the system through alternative ways?

CODE

My current attempt is given in the code below.

(*Setting Known Functions*)

k = 1; hf = 0.01;
Π[h_] := k/hf ((hf/h)^4 - (hf/h)^3)

G = 1; μ = 1; λ = 5; hini = 5;

With[{h1 = h1[x, t], h2 = h2[x, t]},
P1[x_, t_] := -D[h1, x, x] - G D[h2, x, x];
P2[x_, t_] := -G D[h2, x, x] - Π[h2 - h1]];

(*Equations*)

With[{h1 = h1[x, t], h2 = h2[x, t], P1 = P1[x, t], P2 = P2[x, t]},
Q = h1 (h1 - h2) D[P2, x] - D[P1, x]/2 h1^2 -
D[P2, x]/μ h1 (h1/2 - h2);
eqn1 = {D[h1, t] ==
D[(h1^3/3 D[P1, x] - h1^2/2 D[P2, x] (h1 - h2)), x]};
eqn2 = {D[(h2 - h1), t] ==
D[(D[P2, x]/μ (h2^3/3 + h1^3/6 - (h1^2 h2)/2) - Q (h2 - h1)),
x]};

bc = {{D[h1, x] == 0, D[h1, x, x, x] == 0, D[h2, x] == 0,
D[h2, x, x, x] == 0} /.
x -> -10, {D[h1, x] == 0, D[h1, x, x, x] == 0, D[h2, x] == 0,
D[h2, x, x, x] == 0} /. x -> 10};
ic = {h1 == hini,
h2 == hf + hini +
1/3 (λ/π)^(1/2) E^(-λ x^2)} /. t -> 0];

(*Solving the equations*)
timer = 0.100;
Scheme[n : _Integer | {_Integer ..}, o_] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}};

time = 0; points = 5000; difforder = 11;

Monitor[{solh1, solh2} =
NDSolveValue[{eqn1, eqn2, ic, bc}, {h1, h2}, {t, 0, timer}, {x, -10,
10}, Method -> Scheme[points, difforder],
StepMonitor :> (time = t)], time]

• Your problem is x. Jul 7 at 15:41
• @JohnManko Er… what do you mean by x? Jul 8 at 3:43

Yet another problem related to the hyperbolic conservation law. Based on the experience obtained in e.g. here, let's avoid symbolic expansion of D. I'll use pdetoode for the task.

k = 1; hf = 0.01;
Π[h_] := k/hf ((hf/h)^4 - (hf/h)^3)
G = 1; μ = 1; λ = 5; hini = 5;

With[{h1 = h1[x, t], h2 = h2[x, t]}, P1[x_, t_] := -D[h1, x, x] - G D[h2, x, x];
P2[x_, t_] := -G D[h2, x, x] - Π[h2 - h1]];

With[{h1 = h1[x, t], h2 = h2[x, t], P1 = P1[x, t], P2 = P2[x, t], mid1 = mid1[x, t],
mid2 = mid2[x, t], midP2 = midP2[x, t], midP1 = midP1[x, t]},
Q = h1 (h1 - h2) D[midP2, x] - D[midP1, x]/2 h1^2 - D[midP2, x]/μ h1 (h1/2 - h2);
eqnadd1 = mid1 == (h1^3/3 D[midP1, x] - h1^2/2 D[midP2, x] (h1 - h2));
eqnadd2 = mid2 == (D[midP2, x]/μ (h2^3/3 + h1^3/6 - (h1^2 h2)/2) - Q (h2 - h1));
eqn1 = {D[h1, t] == D[mid1, x]};
eqn2 = {D[(h2), t] == D[mid2 + mid1, x]};
ic = {h1 == hini, h2 == hf + hini + 1/3 (λ/π)^(1/2) E^(-λ x^2)} /. t -> 0];


I've removed the original bc. I'll use periodic b.c.s instead in subsequent coding, because

1. Considering the physics background, periodic b.c. should be a good enough approximation for boundary at infinity in this case.

2. Periodic b.c. is convenient to set in pdetoode via last argument.

points = 100; difforder = 2;
domain = {-10, 10}; grid = Array[# &, points, domain];
(* Definition of pdetoode isn't included in this post,
ptoofunc =
pdetoode[{h1, h2, mid1, mid2, midP1, midP2}[x, t], t, grid, difforder, True];
ode = Block[{mid1, mid2, midP1, midP2}, Set @@@ odeadd; ptoofunc@{eqn1, eqn2}];
odeic = ptoofunc@ic;
var = Outer[#[#2] &, {h1, h2}, grid];
timer = 0.1;
sollst = NDSolveValue[{ode, odeic}, var, {t, 0, timer}]; // AbsoluteTiming
(* {1.04179, Null} *)
{h1func, h2func} = rebuild[#, grid, 2] & /@ sollst;

Manipulate[
Plot[{h1func[x, t], h2func[x, t]}, {x, ##}, PlotRange -> {4.5, 5.5}] & @@
domain, {t, 0, timer}]


• Hello! I have recently been playing around with the code and is wondering if it is possible to allow some of the boundary conditions to be defined while letting others remain periodic? Thanks!
– FLP
Jul 10 at 19:33
• @flp Are you sure this is physically correct? Anyway, it's possible, but you cannot make use of the last argument of pdetoode in this case. You need to write down the periodic b.c. in traditional form e.g. u[a, t]==u[b, t], Derivative[1,0][u][a,t]==Derivative[1,0][u][b,t],… then use pdetoode to discretize it. Notice pdetoode cannot directly handle traditional periodic b.c. so you need to discretize it in a manner similar to odeadd above. Jul 11 at 2:16