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]
  • 4
    $\begingroup$ Welcome to Mathematica.SE. Right now, your question lacks any helpful description. The good thing is, you can simply edit your question and include as much description as you like. Please take the time to explain what you are doing where you hit a problem. This will help others to give you advice on how to solve your problem. If you are new here, start by taking time to read the FAQ and the guidelines you find in the "Help Center" at the end of the FAQ. $\endgroup$ – halirutan Jun 11 '18 at 13:35

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.

| improve this answer | |
  • $\begingroup$ Can MMA 11 or 12 now solve this form of problem directly? $\endgroup$ – MOON Jun 9 '19 at 9:42
  • $\begingroup$ @MOON I do not think so. $\endgroup$ – bbgodfrey Jun 10 '19 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

| improve this answer | |
  • $\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 Mar 6 '19 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 Mar 6 '19 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 Jul 23 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 Jul 23 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 Jul 23 at 13:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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