# Conservation of area solving a PDE via finite difference scheme

I have two PDEs that describe the movement of fluid:

1. $h_t + [h^3(1-h)^3((1+\varepsilon h)\sin \theta - \varepsilon h_\theta \cos \theta]_\theta$ = 0
2. $h_t - [h^3(1-h)^3 \varepsilon h_\theta]_\theta$ = 0

You will note that 1) reduces to 2) at $\theta = 0$.

I solve these two PDE's via a finite difference scheme. The first as follows

n = 1001; ϵ = 0.5;
ds = Pi/(n - 1);
ClearAll[dh, t]; endtime = 100;
dh[h_List] := With[{ϵ = 0.5, s = Array[#1 & , n, {-(Pi/2), Pi/2}], d1 = ListCorrelate[{-0.5, 0, 0.5}/ds, #1] & , d2 = ListCorrelate[{1, -2, 1}/ds^2, #1] & },
((-1 + h)^2*h^2*(3*(1 - 2*h)*ϵ*Cos[s]*d1[#1]^2 + (-1 + h)*h*Cos[s]*(1 + h*ϵ - ϵ*d2[#1]) + (-3 + h*(6 + (-5 + 8*h)*ϵ))*d1[#1]*Sin[s]) & )[Join[{0.99}, h, {0.99}]]];

h0 = Table[Piecewise[{{0.01, -0.25510204081632654 < θ < 0.25510204081632654}}, 0.99], {θ, -(Pi/2), Pi/2, ds}];

Clear[h];
sol = With[{h = (h[#1][t] & ) /@ Range[n]}, NDSolveValue[{D[h, t] == dh[h], (h /. t -> 0) == h0}, Head /@ h, {t, 0, endtime}]]; valh = Transpose[(#1["ValuesOnGrid"] & ) /@ sol];
ti = Flatten[sol[[1]]["Grid"]];
s = Array[#1 & , n, {-(Pi/2), Pi/2}]; h = ListInterpolation[valh, {ti, s}, InterpolationOrder -> 1];


This solves reasonably quickly. I check area conservation between $t = 0$ and $t = 100$ with the following:

Table[Integrate[Interpolation[Table[{θ, 0.99 - h[t, θ]}, {θ, -(Pi/2), Pi/2, Pi/3000}], InterpolationOrder -> 1][θ], {θ, -(Pi/2), Pi/2}], {t, {0, 100}}]


An error of $0.96 \%$, which I can reduce to $0.19\%$ by increasing $n$ to 2001. I am happy with this.

However, my problem comes in when solving the second PDE (which I feel in theory should be easier as it is a simplification of the first one). I solve in the same manner (initial condition is flipped due to minus sign and physical reasons).

n = 1001; ϵ = 0.5;
ds = Pi/(n - 1);
ClearAll[dz, t]; endtime = 100;
dz[z_List] := With[{ϵ = 0.5, s = Array[#1 & , n, {-(Pi/2), Pi/2}], d1 = ListCorrelate[{-0.5, 0, 0.5}/ds, #1] & ,
d2 = ListCorrelate[{1, -2, 1}/ds^2, #1] & }, ((-ϵ)*(-1 + z)^2*z^2*((-3 + 6*z)*d1[#1]^2 + (-1 + z)*z*d2[#1]) & )[Join[{0.01}, z, {0.01}]]];

z0 = Table[Piecewise[{{0.99, -0.25510204081632654 < x < 0.25510204081632654}}, 0.01], {x, -(Pi/2), Pi/2, ds}];

Clear[z];
solz = With[{z = (z[#1][t] & ) /@ Range[n]}, NDSolveValue[{D[z, t] == dz[z], (z /. t -> 0) == z0}, Head /@ z, {t, 0, endtime}]];
valz = Transpose[(#1["ValuesOnGrid"] & ) /@ solz];
tiz = Flatten[solz[[1]]["Grid"]];
s = Array[#1 & , n, {-(Pi/2), Pi/2}]; z = ListInterpolation[valz, {tiz, s}, InterpolationOrder -> 1]


Firstly, with $n = 1001$ this takes an age to run on my machine. Secondly, when I check conservation of area for PDE 2.

Table[Integrate[Interpolation[Table[{θ, z[t, θ] - 0.01}, {θ, -(Pi/2), Pi/2, Pi/3000}], InterpolationOrder -> 1][θ], {θ, -(Pi/2), Pi/2}], {t, {0, 100}}]


I get a massive $13.5 \%$ discrepancy between the initial and final error. I could in theory reduce this by increasing $n$, however I tried with $n = 2001$ and the code was running for 8+ hours before I aborted it.

Does anyone have any pointers? My ultimate aim is to have PDE 2 conserving area and solving quickly. Does my FD scheme need some work?

EDIT

By adding the option Method -> {"EquationSimplification" -> "Residual"} to NDSolveValue for equation 2 I can solve for $n = 5001$ points in about 30 minutes. This still gives an intolerable error of $\approx 5 \%$.

• Is there a reason you are using an FDM scheme instead of using NDSolve "out of the box"? Could you also describe what the initial and boundary conditions are for this problem? Spreading of a droplet on some sort of wetting surface? Oct 11, 2016 at 14:12
• @drN I wish to use a FDM as I want to understand each aspect of the code, rather than just putting it into a black box. The initial condition is box function, that - as you correctly state - describes an interface of an initial blob of fluid. The boundary conditions are that far from the blob (ie. ) the interface stays on the surface, so $h(-\pi/2) = h(\pi/2) = 0$. I actually us 0.01 in my code other wise the PDE becomes 0 and nothing moves. Oct 11, 2016 at 14:22
• -1 for showing the wrong equation. There should be a minus sign before the $h_t$ in the second equation. (This mistake ruined my whole night! ) I'll retract the downvote after you correct the mistake anyway… Oct 14, 2016 at 17:26
• @drN Just read the material you linked below, interesting. This and this package might be related, but I can't tell how to use them… Oct 16, 2016 at 4:23
• @ojlm My answer to question 129183 appears to be what you need for your second equation. It produces high accuracy in less than a minute. Thanks to xzczd for pointing out the connection. Dec 11, 2016 at 6:53

Eq. 2 Results

As I noted in a comment above, my answer to question 129183 provides an accurate solution to the second equation in the present question.

where the eleven curves correspond to t == Range[0, 200, 20]. The computation required only 200 grid points, took less than a minute, and conserved area exactly. The method used was essentially the same as that given in the present question but with spatial discretization given by

(ϵ (-(((1 + 1/2 (-z[-1 + i] - z[i]))^3 (-z[-1 + i] + z[i]) (z[-1 + i] + z[i])^3)/(8 d)) +
((1 + 1/2 (-z[i] - z[1 + i]))^3 (-z[i] + z[1 + i]) (z[i] + z[1 + i])^3)/(8 d)))/d


It is easy to see that summing this expression over i leaves only the difference between the first and last terms of the array

((1 + 1/2 (-z[i] - z[1 + i]))^3 (-z[i] + z[1 + i]) (z[i] + z[1 + i])^3)/(8 d)))


which are negligible until well past t == 100. Thus, area conservation is accomplished not by using a large grid but rather by the discretization itself. In contrast, the discretization used in the present question is

-ϵ (-1 + z[i])^2 z[i]^2 (((-3 + 6 z[i]) (-z[-1 + i] + z[1 + i])^2)/(4 d^2) +
((-1 + z[i]) z[i] (z[-1 + i] - 2 z[i] + z[1 + i]))/d^2)


which does not sum conveniently. It is interesting to compare the results above with those from the solutions obtained from the "Edit" at the end of the present question. For n == 5001,

which is similar to the first plot in this answer. However, with the grid size reduced to n == 1001,

which is noticeably different. Results depend on the grid size and only slowly converge to the first plot in this answer. Spatial discretization that mirrors the symmetry of the PDE makes a big difference.

Eq. 1 Results

As requested by ojlm, here are results for the first equation in the question, obtained using a difference scheme similar to that used above for the second equation.

bound = 0.25510204081632654;
upper = 1/100; lower = 1 - upper;
range = {L, R} = {-Pi/2, Pi/2};
endtime = 100;
ϵ = 1/2;
n = 200; d = (R - L)/n;
vars = Table[u[i, t], {i, 2, n}]; u[1, t] = lower; u[n + 1, t] = lower;
eq = Table[dup = (u[i + 1, t] - u[i, t])/d; dum = (u[i, t] - u[i - 1, t])/d;
up = (u[i + 1, t] + u[i, t])/2; um = (u[i, t] + u[i - 1, t])/2;
D[u[i, t], t] == ((up^3 (1 - up)^3) (ϵ dup Cos[L + (i - 1/2) d] -
(1 + ϵ up) Sin[L + (i - 1/2) d]) -
(um^3 (1 - um)^3) (ϵ dum Cos[L + (i - 3/2) d] -
(1 + ϵ um) Sin[L + (i - 3/2) d]))/d, {i, 2, n}];
init = Table[u[i, 0] == Piecewise[{{upper, -bound < L + (i - 1) d < bound}}, lower],
{i, 2, n}];
s = NDSolveValue[{eq, init}, vars, {t, 0, endtime}];
ListLinePlot[Evaluate@Table[Join[{1 - lower},
Table[1 - s[[i - 1]] /. t -> tt, {i, 2, n}], {1 - lower}],
{tt, 0, endtime, endtime/10}], DataRange -> range, PlotRange -> 1]


Conserved area is constant to of order 0.01%. Note that only 200 grid points are required, and increasing the number of grid points to 1000 has no visible impact on the plot.

• Thanks for the reply. So if I understand correctly, you have not differentiated equation 2 before discretising. Suppose now I wished to discretise equation 1 in the same manner, my attempt is as follows (as per your answer 129183) D[u[i, t], t] == (up^3 (1 - up)^3 ((1 + \[Epsilon]*up) Sin[(i - 1) d] - \[Epsilon]* Cos[(i - 1) d] dup) - um^3 (1 - um)^3 ((1 + \[Epsilon]*um) Sin[(i - 1) d] - \[Epsilon]* Cos[(i - 1) d] dum)) /d However this does not work. I presume a problem with the $\sin \theta$ and $\cos \theta$ terms? Dec 14, 2016 at 11:11
• @ojlm Your z0 and h0 are defined differently. Consequently, Eq 2 with initial and boundary conditions is not the same as the limit of Eq 1 with initial and boundary conditions for small θ. Is this intentional? Dec 14, 2016 at 16:12
• @ojlm Try D[u[i, t], t] == ((up^3 (1 - up)^3) (ϵ dup Cos[L + (i - 1/2) d] - (1 + ϵ up) Sin[L + (i - 1/2) d]) - (um^3 (1 - um)^3) (ϵ dum Cos[L + (i - 3/2) d] - (1 + ϵ um) Sin[L + (i - 3/2) d]))/d (with your definition of h0). Dec 14, 2016 at 17:31
• yes, this is intentional, purely due to the physics under consideration during the derivation of equation 1 and equation 2. That was one mistake I had made. Thanks for editing your answer with a full solution for equation 1. It is much much faster than my version. Dec 16, 2016 at 10:45

# Partly NDSolve-based solution

Use a higher even order spatial grid to discretize the PDE to an ODE set seems to be a good approach. The definition of pdetoode can be found here:

bound = 0.25510204081632654;
upper = 99/100; lower = 1 - upper;
range = {L, R} = {-Pi/2, Pi/2};
endtime = 100;
eqn = With[{h = h[t, θ], ϵ = 5/10},
0 == -D[h, t] + D[h^3 (1 - h)^3 ϵ D[h, θ], θ]];
ic = h[0, θ] ==
SimplifyPWToUnitStep@Piecewise[{{upper, -bound < θ < bound}}, lower];
bc = {h[t, L] == lower, h[t, R] == lower};

points = 900;
grid = Array[# &, points, range];
xdifforder = 10;
ptoo = pdetoode[h[t, θ], t, grid, xdifforder];
del = Delete[#, {{1}, {-1}}] &;
var = h /@ grid;
Attributes[h] = NHoldAll;
odeqn = del@ptoo@eqn;
odeic = del@ptoo@ic;
odebc = bc // ptoo;
sollst = NDSolveValue[{odeqn, odeic, odebc} // N, var, {t, 0, endtime},
MaxSteps -> Infinity]; // AbsoluteTiming
(* {19.372995, Null} *)

sol = rebuild[sollst, grid];

(# - #2)/# & @@ Table[NIntegrate[sol[t, θ], {θ, L, R}], {t, {0, endtime}}]
(* 0.00203278 *)

Animate[Plot[sol[t, th], {th, L, R}, PlotRange -> {0, 1}], {t, 0, endtime}]


Notice a even difference order is necessary, if you set e.g. xdifforder = 3, you'll get a erroneous result like:

Result sliced at t = endtime.

# Fully NDSolve-based solution

One may think the solution above is a self implementation of Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> False, "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 900, "MinPoints" -> 900, "DifferenceOrder" -> 10}} option in NDSolve, but it's not true. Check this post for more details. Actually if we directly add this option to NDSolve, we'll obtain a solution quite similar to those we get with odd order in partly NDSolve-based solution!:

Result sliced at t = endtime.

Inspired by the truth, I found that setting an odd "DifferenceOrder" inside NDSolve will resolve the problem:

mol[n_Integer, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}
mol[tf:False|True,sf_:Automatic]:={"MethodOfLines",
"DifferentiateBoundaryConditions"->{tf,"ScaleFactor"->sf}}

sold = NDSolveValue[{eqn, ic, bc}, h, {t, 0, endtime}, {θ, L, R},
Method -> Union[mol@False, mol[900, 9]], MaxSteps -> Infinity]; // AbsoluteTiming
(* {9.774524, Null} *)

(# - #2)/# & @@ Table[NIntegrate[sold[t, θ], {θ, L, R}], {t, {0, endtime}}]
(* 0.00342644 *)

Animate[Plot[sold[t, th], {th, L, R}, PlotRange -> {0, 1}], {t, 0, endtime}]


The animation is similar to the one above so I'd like to omit it here.

• 有你解不出的微分方程吗请问 :)
– yode
Oct 14, 2016 at 18:53
• @yode 有，而且有很多……只是最近弄清楚了NDSolveFiniteDifferenceDerivative的用法所以这几天藉此解出了几个历史遗留问题而已…… Oct 14, 2016 at 18:57
• Very interesting problem: a similarity solution as recommended by Witelski (link 1) (link 2) could help too. Oct 15, 2016 at 15:25
• Thanks for this solution, but it still produces a huge (and unacceptable) error of $\approx 7\%$ Oct 17, 2016 at 12:23
• @ojlm Notice I only used 900 points and it only takes 14 seconds to get the solution. Oct 17, 2016 at 12:29