8
$\begingroup$

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!

Improve this question
$\endgroup$

1 Answer 1

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
8
  • $\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$ Commented 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$ Commented 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$ Commented 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$ Commented Jun 10, 2017 at 12:56
  • $\begingroup$ thanks for your advice. It is true that the b.c.s are the key issues. $\endgroup$ Commented Jun 10, 2017 at 13:51

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.