# Problems with solving PDEs

I am using NDSolve to solve the two equations:

func[x_, t_] := ((v[x, t] - 9/100)*(1 + u[x, t]) - 1/2*u[x, t])/u[x, t];

equ = D[u[x, t], t] - 1/20*func[x, t] + D[u[x, t]^3*(D[u[x, t], x] + D[u[x, t], {x, 3}]), x] == 0;
eqv = D[v[x, t], t] + u[x, t]^2/2*(3*D[u[x, t] + D[u[x, t], {x, 2}], x])*D[v[x, t], x] - 1/10*D[v[x, t], {x, 2}] == (1 - v[x, t])*(func[x, t]^2 - (1/10*Log[(1 - v0)/(1 - v[x, t])])^2);


with the parameters and trial code (it takes about 30s on my i5 laptop):

L = 10 \[Pi]; c = 1/40; v0 = 1/5; tmax = 200;
eqns = {equ, eqv, u[0, t] == u[L, t], v[0, t] == v[L, t], u[x, 0] == 1 + c*Cos[x], v[x, 0] == 1/10};
{solnu, solnv} = NDSolveValue[eqns, {u, v}, {x, 0, L}, {t, 0, tmax}]


By plotting the solution, we found that after $$t>100$$, the solutions begin to disappear partially (see red lines for $$t=101$$):

Plot[{solnu[x, 0], solnu[x, 50], solnu[x, 100], solnu[x, 101]}, {x, 0, L}, PlotRange -> {{0, L}, All}, ImageSize -> 400, PlotPoints -> 80,
Frame -> True, Axes -> False, AspectRatio -> 0.4, PlotStyle -> {Black, Green, Blue, Red}] Plot[{solnv[x, 0], solnv[x, 50], solnv[x, 100], solnv[x, 101]}, {x, 0, L}, PlotRange -> {{0, L}, All}, ImageSize -> 400, Frame -> True,
Axes -> False, AspectRatio -> 0.4, PlotStyle -> {Black, Green, Blue, Red}] My thought:

By evaluating the functions at a certain position, say, $$x=4$$, I can understand that the disappearance of the solutions is because the solutions become complex numbers as $$t$$ increases.

Table[solnu[4, tn], {tn, 0, tmax, 20}]

(*{0.983661, 0.894927, 0.873114, 0.859587, 0.853879, 0.828441,
0.729746 - 0.0799211 I, 0.291994 + 0.0356712 I, 0.170041 - 0.089365 I, 0.166466 - 0.0930112 I, 0.166531 - 0.0929651 I}*)

Table[solnv[4, tn], {tn, 0, tmax, 20}]

(*{0.1, 0.313017, 0.316713, 0.316615, 0.315822, 0.432728,
0.999884 - 0.000152834 I, 0.781819 - 0.00317719 I, 0.999983 - 5.23657*10^-6 I, 1. - 1.08127*10^-7 I, 1. - 6.74084*10^-8 I}*)


In my problem, I need the solutions to be real over the entire domain. I suspected that the solutions with the imaginary parts result from the last term, Log[(1 - v0)/(1 - v[x, t])], of the 2nd equation, because the Log function itself has the domain of definition to be a real function. In addition, the function v[x, t] should not be larger than $$1$$, which is a physical constraint. See also the lower fig, the curves of solnv begin to vanish as the peak approaches to $$1$$).

My first question is:

Is it possible to add a physical constraint of 0<v[x, t]<=vlim[x,t] in the code to make the solution reasonable? Here, vlim[x, t] is determined by the solutions:

vlim[x_, t_] := Exp[4800*(1/350 - solnu[x, t]/((solnv[x, t] - 1/10)*solnu[x, t]))];


Or, is it acceptable if I only take the solutions before they vanish (i.e. $$t\le 100$$). In other words, is it equivalent to the constraint if one only takes the solutions from the current code, which still satisfy the constraint (e.g. by checking v[x, t])?

If you think there is something wrong with my idea, please also leave your comments and suggestion.

My second problem (solved by @Alx):

By defining a function based on the solutions, I need to numerically integrate the function at a given time, say, $$t=50$$. But NIntegrate runs for a very long time ... (aborted by myself). Please help me with defining such a function and then NIntegrate it correctly. Thank you in advance!

funcSol[x_, t_] := ((solnv[x, t] - 1/10)*(1 + solnu[x, t]) - solnu[x, t])/solnu[x, t];

NIntegrate[funcSol[x, 50], {x, 0, L}]

• For the second problem use funcSol[x_?NumericQ, t_?NumericQ]:= ... so that NIntegrate will get numerical input as it expects. – Alx Sep 9 '19 at 15:24
• Is there reason to believe that the exact solution actually is real over the entire domain? – bbgodfrey Sep 10 '19 at 2:58
• @bbgodfrey yes, both solutions u and v should be real because they are for real physical quantities. The problem is how to add the constraint in my first question. Thank you. – user55777 Sep 10 '19 at 3:01
• Physically, how close to 1 can v become? – bbgodfrey Sep 10 '19 at 3:08
• @user55777 From a physical point of view you should use Log[Abs[(1 - v0)/(1 - v[x, t])]] – Alex Trounev Sep 10 '19 at 4:05

## 3 Answers

As noted in the question, the computation fails when v[x, t] >= 1. This is easy to fix by replacing (1 - v[x, t]) by Max[1 - v[x, t], 10^-8] in the right side of eqv. (Results below do not change, if 10^-8 is replaced by 10^-6.) For completeness, here is a plot of the problem term. It is not large in value.

Plot[(1 - z)*(-(1/10*Log[(1 - v0)/(1 - z)])^2), {z, 0, 1},
AxesLabel -> {v, Null}, ImageSize -> Large, LabelStyle -> {Black, Bold, 15}] However, the solution still behaves badly for t greater than about 100. Examining the details of these solutions suggests that the default value of MaxStepFraction is too large to resolve sharp variations in the solution. So, I steadily decreasedMaxStepFraction until the computations proceeded stably at least as far as t == tmax. This required MaxStepFraction -> 1/1000, which led to solution plots,

Plot3D[solnu[x, t], {x, 0, L}, {t, 0, tmax}, AxesLabel -> {x, t, u},
PlotRange -> All, ImageSize -> Large, LabelStyle -> {Black, Bold, 15}, PlotPoints -> 200]
Plot3D[solnv[x, t], {x, 0, L}, {t, 0, tmax}, AxesLabel -> {x, t, v},
ImageSize -> Large, LabelStyle -> {Black, Bold, 15}, PlotPoints -> 200]  A phase transition appears to occur at about t == 140. A computation with MaxStepFraction -> 1/2000 produced similar plots, but with the transition at about t == 120.

Addendum: 1D solution

One might expect that setting c = 0 would yield a solution independent of x.

L = 10 \[Pi]; c = 0; v0 = 1/5; tmax = 200;
eqns = {equ, eqv, u[0, t] == u[L, t], v[0, t] == v[L, t], u[x, 0] == 1 + c*Cos[x],
v[x, 0] == 1/10, WhenEvent[Max@Table[v[L (i - 1/2)/7, t], {i, 7}] > 2/5,
tend = t; "StopIntegration"]};
{solnu, solnv} = NDSolveValue[eqns, {u, v}, {x, 0, L}, {t, 0, tmax},
MaxStepFraction -> 1/4000]


where WhenEvent is used to stop the integration when v[x, t] starts growing toward 1, as now requested in the question. Doing so saves a lot of run time too. The resulting plots, similar to those above but cut off at tend = 126.977, are  Evidently, these oscillations are growing from truncation errors and do not yield a 1D solution. A true 1D solution is obtained by explicitly eliminating x-dependence from the equations.

func[x_, t_] = (((v[x, t] - 9/100)*(1 + u[x, t]) - 1/2*u[x, t])/u[x, t])
/. {u[x, t] -> u[t], v[x, t] -> v[t]};
equ = Unevaluated[D[u[x, t], t] - 1/20*func[x, t] + D[u[x, t]^3*(D[u[x, t], x]
+ D[u[x, t], {x, 3}]), x] == 0] /. {u[x, t] -> u[t], v[x, t] -> v[t]};
eqv = Unevaluated[D[v[x, t], t] + u[x, t]^2/2*(3*D[u[x, t] + D[u[x, t], {x, 2}], x])*
D[v[x, t], x] - 1/10*D[v[x, t], {x, 2}] == (1 - v[x, t])*(func[x, t]^2 -
(1/10*Log[(1 - v0)/(1 - v[x, t])])^2)] /. {u[x, t] -> u[t], v[x, t] -> v[t]};

v0 = 1/5; tmax = 5 10^4;
eqns = {equ, eqv, u == 1, v == 1/10};
{solnu, solnv} = NDSolveValue[eqns, {u, v}, {t, 0, tmax}, MaxStepSize -> 10^-2]

Plot[solnu[t], {t, 0, tmax}, AxesLabel -> {t, u}, ImageSize -> Large,
LabelStyle -> {Black, Bold, 15}, PlotPoints -> 200, PlotRange -> All]
Plot[solnv[t], {t, 0, tmax}, AxesLabel -> {t, v}, ImageSize -> Large,
LabelStyle -> {Black, Bold, 15}, PlotPoints -> 200, PlotRange -> All]
{solnu[tmax], solnv[tmax]}  (* {0.287336, 0.201601} *)


which is very close to the symbolic asymptotic solution, {(9 - 100 v0)/(-59 + 100 v0), v0}. With v0 = 0.2, these are

(* {0.282051, 0.2} *)


We see from the last two plots that the asymptotic solution is not achieved for tmax less than about 10^4, much larger than tmax achievable for the 2D computations. The robustness and relatively long wavelength of the growing 2D oscillations cause me to wonder whether there is an error in the spatial derivative terms of the PDEs.

• Perhaps you can adjust temperal grid and spatial grid separately using e.g. MaxStepFraction -> {1/30, 1/1000}? (More discussion can be found here. ) – xzczd Sep 10 '19 at 7:05
• @xzczd I am aware of this capability, but the variation in x and in t both are extreme. – bbgodfrey Sep 10 '19 at 12:15
• @user55777 The solution rapidly reaches a steady state but apparently becomes numerically unstable at much later time. – bbgodfrey Sep 10 '19 at 12:52
• @bbgodfrey when observing u we can see its peaks increase in time near tmax, why did you say the solution reaches a steady state rapidly? Could you update by adding some descriptions? – user55777 Sep 10 '19 at 13:18
• @user55777 It reaches a steady state at around t == 50 but goes unstable at around t == 100. – bbgodfrey Sep 10 '19 at 13:32

For such a system of nonlinear equations it is necessary to use a special method of solution and a thin grid with a sufficiently large number of cells. Then the concentration function will be automatically limited and the solution will look completely different than on a coarse grid of 25 cells, which is used in the automatic solution method.

func[x_, t_] := ((v[x, t] - 9/100)*(1 + u[x, t]) - 1/2*u[x, t])/
u[x, t];

equ = D[u[x, t], t] - 1/20*func[x, t] +
D[u[x, t]^3*(D[u[x, t], x] + D[u[x, t], {x, 3}]), x] == 0;
eqv = D[v[x, t], t] +
u[x, t]^2/2*(3*D[u[x, t] + D[u[x, t], {x, 2}], x])*
D[v[x, t], x] -
1/10*D[v[x, t], {x, 2}] == (1 -
v[x, t])*(func[x,
t]^2 - (1/10*Log[Abs[(1 - v0)/(1 - v[x, t])]])^2);

L = 10 \[Pi]; c = 1/40; v0 = 1/5; tmax = 200;
eqns = {equ, eqv, u[0, t] == u[L, t], v[0, t] == v[L, t],
u[x, 0] == 1 + c*Cos[x], v[x, 0] == 1/10};
{solnu, solnv} =
NDSolveValue[eqns, {u, v}, {x, 0, L}, {t, 0, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 137, "MaxPoints" -> 137,
"DifferenceOrder" -> "Pseudospectral"}}, MaxSteps -> 10^6]

Plot[{solnu[x, 0], solnu[x, 50], solnu[x, 100], solnu[x, 101]}, {x, 0,
L}, PlotRange -> {{0, L}, All}, ImageSize -> 400, PlotPoints -> 80,
Frame -> True, Axes -> False, AspectRatio -> 0.4,
PlotStyle -> {Black, Green, Blue, Red}]

Plot[{solnv[x, 0], solnv[x, 50], solnv[x, 100], solnv[x, 101]}, {x, 0,
L}, PlotRange -> {{0, L}, All}, ImageSize -> 400, Frame -> True,
Axes -> False, AspectRatio -> 0.4,
PlotStyle -> {Black, Green, Blue, Red}] {Plot3D[solnu[x, t], {x, 0, L}, {t, 0, tmax}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic],
Plot3D[solnv[x, t], {x, 0, L}, {t, 0, tmax}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic,
PlotRange -> All]} This method can be used up to t=270 without loss of stability, but without violation of condition v<1 only up to t=230. The result is similar to what @bbgodfrey got. • Using Pseudospectral does delay the onset of rapidly growing oscillations a bit but it does not eliminate them. In fact, to reach the steady state that the OP apparently wants would require tmax of order 10,000! The robustness and relatively long wavelength of these growing oscillations cause me to wonder whether there is an error in the spatial derivative terms of the PDEs. – bbgodfrey Sep 10 '19 at 16:51
• I agree that in this case the instability is suppressed at t<200, but grows rapidly at t>250. Consequently we still need to limit the concentration. Obviously, the model is unsuitable for modeling real effects. – Alex Trounev Sep 11 '19 at 5:09

This problem seems to be 4th one I find in this site, on which the default spatial discretization of NDSolve doesn't work well even if nothing like shock wave involves in. (BTW here are 1st one, 2nd one and 3rd one. )

The following method gives stable solution up to tmax = 200. The key point is the discretization of D[u[x, t]^3 (D[u[x, t], x] + D[u[x, t], {x, 3}]), x] in equ. We need to discretize it before the automatic symbolic differentiation happens. In other words, terms like $$\frac{\partial (f g)}{\partial x}$$ should be discretized as

$$\frac{\partial (f g)}{\partial x}\Biggl|_{x=x_i} \approx \frac{f(x_i+h) g(x_i+h)-f(x_i-h)g(x_i-h)}{2h}$$

Notice this scheme is different from:

$$\frac{\partial (f g)}{\partial x}\Biggl|_{x=x_i}=\left(f\frac{\partial g}{\partial x}+g\frac{\partial f}{\partial x}\right)\Biggl|_{x=x_i}\approx f(x_i) \frac{g(x_i+h)-g(x_i-h)}{2h}+g(x_i)\frac{f(x_i+h)-f(x_i-h)}{2h}$$

I'll use pdetoode for discretization, notice the periodic b.c.s are set by 5th argument of pdetoode:

func[x_, t_] := ((v[x, t] - 9/100) (1 + u[x, t]) - 1/2 u[x, t])/u[x, t];

equ = D[u[x, t], t] - 1/20 func[x, t] + D[mid[x, t], x] == 0;
midexpr = u[x, t]^3 (D[u[x, t], x] + D[u[x, t], {x, 3}]);

eqv = D[v[x, t], t] + u[x, t]^2/2 (3 D[u[x, t] + D[u[x, t], {x, 2}], x]) D[v[x, t], x] -
1/10 D[v[x, t], {x, 2}] == (1 -
v[x, t]) (func[x, t]^2 - (1/10 Log[(1 - v0)/(1 - v[x, t])])^2);

L = 10 π; c = 1/40; v0 = 1/5; tmax = 200;

ic = {u[x, 0] == 1 + c*Cos[x], v[x, 0] == 1/10};

points = 300; difforder = 2; domain = {0, L}; grid = Array[# &, points, domain];

(* Definition of pdetoode isn't included in this post,
please find it in the link above. *)
ptoofunc = pdetoode[{u, v, mid}[x, t], t, grid, difforder, True];

rule = Thread[ptoofunc@mid[x, t] -> ptoofunc@midexpr];
ode = ptoofunc@{equ, eqv} /. rule;
odeic = ptoofunc@ic;
sollst = NDSolveValue[{ode, odeic}, Outer[#@#2 &, {u, v}, grid], {t, 0, tmax}(*,
Method->{"EquationSimplification"->"Solve"}*)]; // AbsoluteTiming
(* {2.0083, Null} *)
{solu, solv} = rebuild[#, grid, -1] & /@ sollst

With[{solnu = solu, solnv = solv},
Plot3D[solnu[x, t], {x, 0, L}, {t, 0, tmax}, AxesLabel -> {x, t, u}, PlotRange -> All,
PlotPoints -> 100]] With[{solnu = solu, solnv = solv},
Plot3D[solnv[x, t], {x, 0, L}, {t, 0, tmax}, AxesLabel -> {x, t, v}, PlotPoints -> 100,
PlotRange -> All]] You may have noticed it's not necessary to modify the Log[…] term in eqv with this method.

Unfortunately, for larger tmax e.g. tmax = 300, the method still becomes unstable, perhaps a even better difference scheme helps, but I haven't found any so far.

• it seems that you never use the boundary conditions? – user55777 Sep 11 '19 at 11:31
• @user55777 Periodic b.c. is set by 5th argument of pdetoode. – xzczd Sep 11 '19 at 11:32
• thank you for the reply. I need time to understand the method. Btw, when using a spline ic: iniu[x_] = 1 + c*BSplineFunction[RandomReal[{-1, 1}, 20], SplineClosed -> True][x/L]; I got the error "Initial condition is not a number or a rectangular array of numbers", even I set its attribute with SetAttributes[iniu0, Listable]; and used iniu0[x_] = iniu[x]; in ic = {u[x, 0] == iniu0[x], v[x, 0]...};? Please give some suggestion. – user55777 Sep 11 '19 at 12:10
• @user55777 Just check what's inside your ic, and re-visit how I make use of Listable. – xzczd Sep 11 '19 at 12:16
• it is ic = {u[x, 0] == iniu[x], v[x, 0] == 1/10}, where iniu[x] is a spline function. Sorry, i cannot really understand your question. But I understood that we should add Listable to the list of attributes of the symbol iniu, as you did. – user55777 Sep 11 '19 at 12:40