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]
x
. $\endgroup$x
? $\endgroup$