Xmax = 5;

Tmax = 5;

eq1 = D[u[x, t], t] == D[u[x, t], x, x] + (x - NIntegrate[x u[x, t], {x, 0, Xmax}]) u[x, t]

iv5 = {u[x, 0] == 2/(Sqrt[Pi]*Exp[x^2])};

bcs = {u[0, t] == 2/Sqrt[Pi], u[Xmax, t] == 0};

s10 = NDSolve[ {eq1, iv5, bcs} , {u[x, t] } , {x, 0, Xmax} , {t, 0, Tmax} ];

y = Table[u[x, t] /. s10, {x, 0, Xmax}, {t, 0, Tmax}]

Plot3D[u[x, t] /. s10, {x, 0, Xmax}, {t, 0, Tmax}, PlotRange -> All]

2 Answers 2


NDSolve is not capable of solving this sort of problem as a PDE. Thus, it is necessary to perform the computation by discretizing the PDE in x. This procedure is discussed in Introduction to Method of Lines. A while ago, I solved a somewhat similar problem, 78493, that involved an integral over u in one of the boundary conditions. Here, the integral of x u enters into the PDE itself. The code nonetheless resembles that in the earlier problem.

xmax = 5; tmax = 5;
n = 100; h = xmax/n;
U[t_] = Table[u[i][t], {i, 1, n + 1}];
xtab = Table[(i - 1) h, {i, 1, n + 1}];

creates the list of dependent variables and their corresponding positions in x. Then,

usum = xtab.U[t] h;
stab = Join[{0}, Thread[(xtab - usum) U[t]][[2 ;; n]], {0}];

generates the result of the discretized integral of x u and constructs the source term. (x - NIntegrate[x u[x, t], {x, 0, xmax}]) u[x, t]. (Observe that the source term is not applied to the boundary equations.) Next,

eqns = Thread[D[U[t], t] == stab + Join[{0}, ListCorrelate[{1, -2, 1}/h^2, U[t]], {0}]];
initc = Thread[U[0] == 2/(Sqrt[Pi]*Exp[xtab^2])];
lines = NDSolveValue[{eqns, initc}, U[t], {t, 0, tmax}] // Flatten;

constructs the coupled ODEs that represent the PDE, the initial conditions, and solves them numerically. The result is,

ParametricPlot3D[Evaluate@Thread[{xtab, t, lines}], {t, 0, tmax}, 
    PlotRange -> All, AxesLabel -> {"x", "t", "u"}, BoxRatios -> {2, 2, 1}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

enter image description here

3D Plot

In response to the comment below, a smooth 3D surface plot can be obtained by

Flatten[Table[{(m - 1) h, t, lines[[m]]}, {m, n + 1}, {t, 0, tmax, .1}], 1];
ListPlot3D[%, AxesLabel -> {"x", "t", "u"}, BoxRatios -> {2, 2, 1}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

enter image description here

If desired, an Interpolatingfunction can be obtained in the same way.

  • $\begingroup$ Can MMA 11 or 12 now solve this form of problem directly? $\endgroup$
    – MOON
    Commented Jun 9, 2019 at 9:42
  • $\begingroup$ @MOON I do not think so. $\endgroup$
    – bbgodfrey
    Commented Jun 10, 2019 at 0:37

Here's a hacky way to do it, based on inspecting this:

NDSolve`ProcessEquations[{eq1, iv5, bcs}, {u[x, t]}, {x, 0, Xmax}, {t, 0, Tmax}]

There's a couple of places where you see MapThread[rhsFN, data, 1], that maps the right-hand side of the first-orderized differential equation onto the state data. Since in this case, the RHS is vectorized, we can override MapThread and apply the RHS directly with a integration slipped in for a dummy function int[]. Maybe not the safest way to do this, but I thought it was cool enough to share.

Xmax = 5;
Tmax = 5;
eq1 = D[u[x, t], t] == D[u[x, t], x, x] + (x - int[u[x, t], x, t]) u[x, t]; (* N.B. *)
iv5 = {u[x, 0] == 2/(Sqrt[Pi]*Exp[x^2])};
bcs = {u[0, t] == 2/Sqrt[Pi], u[Xmax, t] == 0};

Block[{int, xx},
  int[u_, x_, t_ /; t == 0] = (* IC - fools ProcessEquations, thinks int[] a good num.fn. *)
    NIntegrate[2/(Sqrt[Pi]*Exp[x^2]), {x, 0, Xmax}];
  int[u_?VectorQ, x_?VectorQ, t_?NumericQ] := 
    Integrate[Interpolation[Transpose@{x, x*u}][xx], xx] /. xx -> Xmax;
   MapThread[f_, data_, 1] /; ! FreeQ[f, int] := f @@ data;
   s10 = NDSolve[{eq1, iv5, bcs}, {u[x, t]}, {x, 0, Xmax}, {t, 0, Tmax}];

Plot3D[u[x, t] /. s10, {x, 0, 5}, {t, 0, 5}]

enter image description here

  • $\begingroup$ Certainly, I am impressed (+1). Can your method be applied to 192123? I can solve it by the method I used in the question above, but I have not found the time to write it up yet. $\endgroup$
    – bbgodfrey
    Commented Mar 6, 2019 at 5:15
  • $\begingroup$ @bbgodfrey I gave iit a shot. Since the integral isn't vectorized/listable, I had to adjust some things. Also MapThread got used in ProcessEquations and my overload of it interfered with that stage.. $\endgroup$
    – Michael E2
    Commented Mar 6, 2019 at 17:13
  • $\begingroup$ @MichaelE2 For a slight extension, if we have 2 integrals, say, int1[u[x, t], x, t] and int2[D[u[x,t],{x,2}], x, t], is it correct to use MapThread[f_, data_, 1] /; ! FreeQ[f, int1] := f @@ data; MapThread[f_, data_, 1] /; ! FreeQ[f, int2] := f @@ data; to slip both of them in? $\endgroup$
    – user55777
    Commented Jul 23, 2020 at 8:05
  • $\begingroup$ @user55777 You shouldn't have to change anything except int to int1 (or int2 or int1 | int2). Substitute MapThread[f_, data_, 1] /; ! FreeQ[f, int] := Return[f, NDSolve] in my example above (and omit the Plot). Then s10 will be set to what f (= rhsFN) is. In your case, rhsFN should contain both int1 and int2. Now, one limitation is that the rhsFN must evaluate correctly when fed a vector of x values instead of being mapped over each individual one. $\endgroup$
    – Michael E2
    Commented Jul 23, 2020 at 12:40
  • $\begingroup$ @MichaelE2 Among others, I don't under the use of Return. In the documentation, it only takes one argument (i.e., Return[expr]), which is used to return the value of expr from a function... $\endgroup$
    – user55777
    Commented Jul 23, 2020 at 13:13

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.