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

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!

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[NIntegrateGaussRuleData[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.

• Hi, @xzczd, thanks. However, I don't agree with your D on my PDE. As can be seen, $\lambda(t)$ is a parameter depending on $t$, so when you take derivative with respect to x, it will be gone totally. In fact, $\lambda(t)$ is a critical parameter acting as a Lagrange multiplier. Mathematically, after taking derivative, the nature of PDE will be changed usually. For example,DSolve[y'[x] + y[x] == a f[t], y[x], x] and DSolve[D[y'[x] + y[x] == a f[t], x], y[x], x] will give different answers. Thus, I think they are different from each other. Feb 13, 2015 at 15:42
• @lxy No, they're the same, when being a definite solution problems. The solutions of your 2 examples are different, but it's only true for general solution. When you set proper b.c.s to get a unique solution, they'll be the same, try DSolve[{y'[x] + y[x] == a f[t], y[0] == 0}, y[x], x] and DSolve[{D[y'[x] + y[x] == a f[t], x], y[0] == 0, y'[0] == a f[t]}, y[x], x]. (y'[0] == a f[t] can be easily deduced from the first equation set.) And this is exactly your case. Of course if Derivative[2, 0][h][0, t] == Derivative[2, 0][h][L, t] is mistakenly given, then it's another story. Feb 13, 2015 at 16:27
• Hi, @xzczd, a simple extension of your nice trick. How should I D the equation to eliminate the integral in the differential equation if the problem becomes the following 2D version w.r.t.$h(x,y,t)$: $h_t=(h_{xx}+h_{yy})-V_h-\lambda(t)$ with $\lambda(t)=-l^{-2}\int_0^l\int_0^l V_h dxdy$, where $V_h$ is again a corresponding function of $h(x,y,t)$. I am wondering which one is correct: D[..., x, y] or D[...,{{x,y}}]? I know the latter gives a gradient vector? Thanks! Jun 10, 2017 at 11:45
• @Enter I think a single D[…, x] or D[…, y]` is enough, because $\lambda(t)$ is independent of $x$ and $y$, so after differentiate the equation once, the $\lambda(t)$ term no longer exists. (Of course the b.c. needs to be determined carefully then. ) Jun 10, 2017 at 12:56
• thanks for your advice. It is true that the b.c.s are the key issues. Jun 10, 2017 at 13:51