# NDSolveValue gives up without warning or solution on PDEtoODE system

I have the following pair of PDEs that I want to solve in the half-space x>=0:  I followed this post to decompose the 4th-order height equation into two 2nd-order ones to respect the hyperbolic conservation law and it solves beautifully. However, I would like to solve directly for the pressure p = D[h,{x,2}] (instead of placing it in Block) as I require 2nd-order spatial derivatives of it (ie 4th spatial derivative of h) for further computations. With my below code, NDSolveValue just exits without any warning. sol simply remains as the NDSolveValue call.

Providing an ic to p doesn't resolve the issue either.

1. Could it be because there are no time derivatives of p in eqp?

Of course, as an alternative, I can easily increase the InterpolationOrder of the solutions, but it would be interesting to solve directly for the curvature too.

Thank you so much!

Code: (pdetoode can be found here)

ϵ = 0.1;
b = 0.1;
Ca = 10^-3;
Pe = 10^6;
α = 0.1;
σ[Γ_] := (α + 1)/(1 + (((α + 1)/α)^(1/3) - 1)*Γ)^3 - α
Hfunc[x_] := 1/2*(1 + Tanh[100*x])
hinit = (1 + b - x^2)*Hfunc[1 - x] + b*Hfunc[x - 1];
Γinit = Hfunc[1 - x];

uniform = Subdivide[30, 2000];
rb = uniform[[-1]];

With[{h = h[x, t], Γ = Γ[x, t],
p = p[x, t], qh = qh[x, t], q2 = q2[x, t]},
eqp = p == D[h, {x, 2}];
eqqh = qh == h^2/2*D[σ[Γ], x] + Ca*h^3/3*p;
eqh = D[h, t] == -D[qh, x];
eqΓ = D[Γ, t] == -D[Γ*h*D[σ[Γ], x] + Ca*(Γ*h^2)/2*p, x] + 1/Pe*D[Γ, {x, 2}];
ic = {h == hinit, Γ == Γinit} /. t -> 0;
bc = Flatten[{{D[h, x] == 0, D[h, {x, 3}] == 0, D[Γ, x] == 0} /. x -> 0, {D[h, x] == 0, D[h, {x, 3}] == 0, D[Γ, x] == 0} /. x -> rb}]; (*bc if I do not solve directly for pressure*)
bcp = Flatten[{{D[h, x] == 0, D[p, x] == 0, D[Γ, x] == 0} /. x -> 0, {D[h, x] == 0, D[p, x] == 0, D[Γ, x] == 0} /. x -> rb}];
]

domain = {0, rb}
grid = uniform;
difforder = 2;
remove1 = #[[2 ;; -2]] &;

tfunc = pdetoode[{h, Γ, p, qh}[x, t], t, grid,
difforder, False];
odemid = Block[{p}, Map[tfunc, {eqqh}, {2}]];
odep = tfunc@eqp // remove1;
odeh = Block[{qh}, Set @@@ odemid; tfunc@eqh] // remove1;
odeΓ = Block[{qh}, Set @@@ odemid; tfunc@eqΓ] // remove1;
ode = {odep, odeh, odeΓ};

odeic = tfunc@ic;
With[{sf = 1}, odebc = diffbc[t, sf]@bcp // tfunc];
var = Outer[#[#2] &, {p, h, Γ}, grid];

Monitor[sol = NDSolveValue[{ode, odeic, odebc}, var, {t, 0, 1}, Method -> {"EquationSimplification" -> "Solve"}, EvaluationMonitor :> (time = t)], time];

{psol, hsol, csol} = rebuild[#, grid, 2] & /@ sol;

(*visualisation*)
Manipulate[
Plot[{hsol[r, t], csol[r, t]}, {x, ##}, PlotPoints -> 100] & @@ domain, {t, 0, 1, 0.001}]

Expected output for h (blue) and Γ (red)

At the request of @Alex Trounev, here is working code to reproduce the above plot of h and Γ, the difference is only in the discretisation procedure:

remove2 = #[[3 ;; -3]] &;
remove1 = #[[2 ;; -2]] &;
tfunc = pdetoode[{h, Γ, p, qh}[x, t], t, grid, difforder, False];
odemid = Map[tfunc, {eqp, eqqh}, {2}];
odeh = Block[{p,qh}, Set @@@ odemid; tfunc@eqh] // remove2;
odeΓ = Block[{p,qh}, Set @@@ odemid; tfunc@eqΓ] // remove1;
ode = {odeh, odeΓ};

odeic = tfunc@ic;
With[{sf = 1}, odebc = diffbc[t, sf]@bc // tfunc];
var = Outer[#[#2] &, {h, Γ}, grid];

Monitor[sol = NDSolveValue[{ode, odeic, odebc}, var, {t, 0, 1}, Method -> {"EquationSimplification" -> "MassMatrix"}, EvaluationMonitor :> (time = t)], time];

{hsol, csol} = rebuild[#, grid, 2] & /@ sol;

Update: Testing out the method proposed by @xzcd, I obtained the following solution for the pressure. My question though, is why is the value so different compared to the one obtained from differentiating hsol directly? (psol is the new solution, psoldif=D[hsol[x,t],{x,2}] is the one from differentiation) Am I missing something?

Manipulate[Plot[psol[x, t], {x, 0, 10},PlotPoints -> 100], {t, 0, 5, 0.01}]

Manipulate[Plot[psoldif /. t -> a, {x, 0, 10}, PlotPoints -> 100], {a, 0, 5, 0.0001}]

In response to my previous statement that it could be due to the error message I obtained during discretisation of the IC of p:General::munfl: 4.9226*10^-155^2 is too small to represent as a normalized machine number; precision may be lost. I doubt that it is the source of the issue because it occurs only far away from the initial "drop" deposited where the pressure is very very small.

Update 2 I was careless and did not plot the full PlotRange, thanks @xzczd for pointing that out. However, I noticed that the pressure solution (light blue) is greater than that rebuilt using fdd (red): why is this so?

• Could you upload working code you used to plot picture above? Commented Jan 11, 2023 at 4:26
• @AlexTrounev sorry for the delay, I have uploaded the code! Commented Jan 13, 2023 at 10:17
• (-1) The update is careless. Commented Jan 13, 2023 at 11:04
• @xzczd sorry what do you mean? Commented Jan 13, 2023 at 11:08
• 1. Definition of remove2 is missing. 2. The code text of the 2 Manipulate[…]s are not given. 3. The General::munfl is posted without any explanation. There may be more. Commented Jan 13, 2023 at 11:16

Could it be because there are no time derivatives of p in eqp?

Yes. One of the weakness of DAE solver of NDSolve is, "EquationSimplification" -> "Solve" never tries to differentiate the algebraic equation, it simply Solve it. When certain derivative term is missing, it returns the input silently. A simpler example illustrating the issue:

NDSolve[{x'[t] == 0, y[t] == 0, x[0] == 0}, {x, y}, {t, 0, 1},
Method -> "EquationSimplification" -> "Solve"]
(* Input returned *)

By modifying definition of ic to

ic = {h == hinit, Γ == Γinit, p == D[hinit, x, x]} /. t -> 0;

and definition of odep to

odep = tfunc@D[eqp, t] // remove1;

and set

Method -> {"EquationSimplification" -> "MassMatrix"}
(* Not necessary in principle, but "Solve" is too slow for large system *)

for NDSolve, the problem will be resolved. (Though not necessary, uniform can be reduced to e.g. uniform = Subdivide[2, 500]; to save some time. )

However, it doesn't make much sense to "solve directly for the pressure p". Perhaps what's in your mind is something like:

How to find a numerical antiderivative with NIntegrate methods?

↑ In this case the method makes a difference (to some degree) because we've turned to the adapative ODE solver to integrate, while in our case, no matter how you handle the p == D[h, {x, 2}], the spatial derivative will be calculated by fdd (NDSolveFiniteDifferenceDerivative) as long as you're using pdetoode. If you just want to avoid differentiating the InterpolatingFunction of h, calculate p with

hlst = DeveloperToPackedArray@#["ValuesOnGrid"] & /@ sol[[1]];
gridlst = {grid, Flatten[sol][[1]]["Coordinates"][[1]]};
psol = ListInterpolation[fdd[{2, 0}, gridlst, hlst, DifferenceOrder -> difforder],
gridlst];

is enough.

• Thank you for the answer (and for all the edits, I apologise for them)! Sorry for the delay, I tested it out but have a new question regarding the difference in the solutions output, could you take a look at that? Commented Jan 13, 2023 at 10:19