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!


1 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]

Mathematica graphics

The code for numeric integration is modified from this answer.

  • $\begingroup$ 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. $\endgroup$
    – Enter
    Feb 13, 2015 at 15:42
  • $\begingroup$ @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. $\endgroup$
    – xzczd
    Feb 13, 2015 at 16:27
  • $\begingroup$ 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! $\endgroup$
    – Enter
    Jun 10, 2017 at 11:45
  • $\begingroup$ @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. ) $\endgroup$
    – xzczd
    Jun 10, 2017 at 12:56
  • $\begingroup$ thanks for your advice. It is true that the b.c.s are the key issues. $\endgroup$
    – Enter
    Jun 10, 2017 at 13:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.