# Integrating over $x$ in numerically solving a partial integrodifferential equation

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]


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}]];
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}] 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}] If desired, an Interpolatingfunction can be obtained in the same way.

• Can MMA 11 or 12 now solve this form of problem directly?
– MOON
Jun 9 '19 at 9:42
• @MOON I do not think so. Jun 10 '19 at 0:37

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

NDSolveProcessEquations[{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; • @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.. Mar 6 '19 at 17:13
• @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? Jul 23 '20 at 8:05
• @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. Jul 23 '20 at 12:40
• @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... Jul 23 '20 at 13:13