How can I invoke the solution of NDSolve to determine a parameter in my equation just inside NDSlove?

Question

I am trying to solve a differential equation by NDSlove for $h(x,t)$. It reads $$h_t=h_{xx}-V_h-\lambda(t)$$ where $V_h$ is a given function of $h(x,t)$ denoted by vdh[x_,t_] in my code, $\lambda(t)$ is a time-dependent parameter which is determined by a definite integration.

$$\lambda(t)=-\frac{1}{L}\int_0^L V_h dx$$ In addition, $h(x,t)$ subjects to simple periodic boundary condition and initial condition, which is defined on the periodic interval $[0,L]$, see my code.

First, Defining constants, function vdh[x_,t_] and $\lambda(t)$

ClearAll[x, t]
amp = 5787/1000;
tmax = 1000;
L = 20;
vdh[x_?NumericQ, t_?NumericQ] := 1/(1 + h[x, t]) - h[x, t] + Log[h[x, t]/(1 + h[x, t])];
\[Lambda][t_Real] := -(1/L)*NIntegrate[vdh[x, t], {x, 0, L}];

Then, the main part

 myfunel = First[h /. NDSolve[{
 D[h[x, t], t] == D[h[x, t], {x, 2}] - vdh[x, t] - \[Lambda][t],
 h[0, t] == h[L, t],
 Derivative[1, 0][h][0, t] == Derivative[1, 0][h][L, t],
 Derivative[2, 0][h][0, t] == Derivative[2, 0][h][L, t],
 h[x, 0] == 1 + 1/Sqrt[2*\[Pi]]*Exp[-((x - 10)^2/2)],
 WhenEvent[h[L/2, t] >= amp, "StopIntegration"]},
 h, {x, 0, L}, {t, 0, tmax},
 Method -> {"MethodOfLines", 
 "SpatialDiscretization" -> {"TensorProductGrid", 
 "MinPoints" -> 101, "MaxPoints" -> 101, 
 "DifferenceOrder" -> 4}}, AccuracyGoal -> 20, 
 WorkingPrecision -> MachinePrecision, StepMonitor :> Print[t]]]

Which dose not work. NDSolve prompts a mass of errors, such as:

NIntegrate::inumr: "The integrand 1-h[x,y,t]-h[x,y,t]/(1+h[x,y,t])+Log[h[x,y,t]/(1+h[x,y,t])] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,20},{0,20}}."

Note 1: In every time step, I want to determine $\lambda(t) $ numerically using \[Lambda][t_Real], where lies the problem. It is well known that when NDSolve uses (N)Integrate internally, we will suffer from such kind of difficult very often. However, the equation I want to solve contain parameter $\lambda(t)$ found by integrating a function of the intermediate solution of my NDSlove. Fortunately, the NIntegrate does not involve t variable, which can be regarded as the average value of ${dV \over dh}$. I believe this point causes these problems. I'm not sure whether I have some other mistakes.

Note 2: I also need to see the evolution of the following function:

f[t_] := NIntegrate[h[x, t]^2 + 1/2*(D[h[x, t], x])^2, {x, 0, L}]
Plot[Evaluate[f[t]/.myfunel], {t, 0, tstop}, PlotRange -> All]

where tstop is the stop moment of the main part. I have used WhenEvent to stop my integration, but how to get the stop time in my main part rather than check StepMonitor :> Print[t] after its stop?

Note 3: I also tried LaplaceTransform, but it is obvious that this equation does nto has a closed form solution, like this problem

Any ideas? Also any suggestion for code improvements is more than welcome!

Answer 1

This problem can be solved by disretizing the PDE to a set of ODEs. This approach was too cumbersome for me in 2015, but much easier nowadays because I've created pdetoode.

vdh = {x, t} \[Function] 1/(1 + h[x, t]) - h[x, t] + Log[h[x, t]/(1 + h[x, t])];
λ = t \[Function] -(1/L) int@t;
eqn = D[h[x, t], t] == D[h[x, t], {x, 2}] - vdh[x, t] - λ[t];
ic = h[x, 0] == 1 + 1/Sqrt[2*π]*Exp[-((x - 10)^2/2)];

amp = 5787/1000; tmax = 1000; L = 20;

points = 25; difforder = 4; domain = {0, L};

{nodes, weights} = Most[NIntegrate`GaussRuleData[points, MachinePrecision]];

midgrid = Rescale[nodes, {0, 1}, domain];

intrule = int@t -> -Subtract @@ domain weights.Map[Function[x, vdh[x, t]], midgrid] /. 
   h[x_, t] :> h[x][t];

grid = Flatten[{0, midgrid, L}];

(*Definition of pdetoode isn't included in this post,please find it in the link above.*)


ptoofunc = pdetoode[h[x, t], t, grid, difforder, True];

ode = ptoofunc[eqn] /. intrule;
odeic = ptoofunc@ic;
hmid = h@grid[[Length@grid/2 // Round]];

sollst = NDSolveValue[{ode, odeic, 
    WhenEvent[hmid[t] >= amp // Evaluate, "StopIntegration"]}, h /@ grid, {t, 0, tmax}];
{{tmin, trealmax}} = sollst[[1]]["Domain"];
sol = rebuild[sollst, grid, -1];
Plot3D[sol[x, t], {x, 0, L}, {t, tmin, trealmax}, PlotRange -> All]

The code for numeric integration is modified from this answer.