# Solving an integro-differential equation with Mathematica

I try to solve a nonlinear integro-differential equation with this code. Here i used a periodic condition.

L=10; tmax = 2;

NDSolve[{D[u[x, t], t] + u[x, t]*D[u[x, t], x] + D[u[x, t], {x, 2}] +
D[u[x, t], {x, 4}] + 1/(2 L)*NIntegrate[D[u[xp, t],{xp, 3}]*Cot[\[Pi](x - xp)/(2*L)], {xp, -L, x, L}, Method -> {"PrincipalValue"}] == 0,
u[-L, t] == u[L, t], u[x, 0] == 0.1*Cos[\[Pi]/L*x]}, u, {x, -L, L}, {t, 0, tmax}]


which gives me

NDSolve::delpde:Delay partial differential equations are not currently supported by NDSolve"

The warning is understandable because the function u[xp, t] is still unknow when NIntegrate is evaluated. Note that we should use PrincipalValue here in NIntegrate because there is a singularity at $$x=xp$$, which has been specified in the integration range.

• The method I used to solve the earlier problem cited in the question cannot be employed here, because the PDE is nonlinear. It may, however, be possible to solve the equation using the method outlined here, although not without a great deal of effort. – bbgodfrey Feb 24 '19 at 20:08
• My answer to question 175080 may be helpful, although it is a bit simpler. – bbgodfrey Feb 25 '19 at 5:13
• If you solve this problem, please post your solution as an answer. – bbgodfrey Feb 25 '19 at 23:13
• Perhaps, I can try near the end of the week, if you have not already done so. – bbgodfrey Feb 26 '19 at 4:26
• I think it would be advisable to change the order of D[NIntegrate[...]] to NIntegrate[D[...]]. This would mean replacing D[Integrate[u[xp, t]*Cot[\[Pi] (x - xp)/(2*L)], {xp, -L, L}], {x, 3}] with NIntegrate[-((\[Pi]^3 (2 + Cos[(\[Pi] (x - xp))/L]) Csc[(\[Pi] (x - xp))/(2 L)]^4 u[xp, t])/(4 L^3)), {xp, -L, L}], where the integrand in the second form was found with D[u[xp, t]*Cot[\[Pi] (x - xp)/(2*L)], {x, 3}] // FullSimplify. – Roman Feb 26 '19 at 13:31

Based on the hacky way I used in my answer here; I had to split up the NDSolve process, so as not to redefine MapThread too soon:

L = 10; tmax = 2;
sys = {D[u[x, t], t] + u[x, t]*D[u[x, t], x] + D[u[x, t], {x, 2}] +
D[u[x, t], {x, 4}] + 1/(2 L)*int[D[u[x, t], {x, 3}], x, t] == 0,
u[-L, t] == u[L, t], u[x, 0] == 0.1*Cos[\[Pi]/L*x]};
periodize[data_] := Append[data, {N@L, data[[1, 2]]}]; (* for periodic interpolation *)
Block[{int},
(* the integral *)
int[uppp_, x_?NumericQ, t_ /; t == 0] := (cnt++;
NIntegrate[
D[0.1*Cos[\[Pi]/L*xp], {xp, 3}]*Cot[\[Pi] (x - xp)/(2*L)],
{xp, x - L, x, x + L},
Method -> {"InterpolationPointsSubdivision", Method -> "PrincipalValue"},
PrecisionGoal -> 8, MaxRecursion -> 20, AccuracyGoal -> 20]);
int[uppp_?VectorQ, xv_?VectorQ, t_] := Function[x,
cnt++;
NIntegrate[
Interpolation[periodize@Transpose@{xv, uppp}, xp,
PeriodicInterpolation -> True]*Cot[\[Pi] (x - xp)/(2*L)],
{xp, x - L, x, x + L},
Method -> {"InterpolationPointsSubdivision", Method -> "PrincipalValue"},
PrecisionGoal -> 8, MaxRecursion -> 20] (* adjust to suit *)
] /@ xv;
(* monitor while integrating pde *)
Clear[foo];
cnt = 0;
PrintTemporary@Dynamic@{foo, cnt, Clock[Infinity]};
(* broken down NDSolve call *)
InternalInheritedBlock[{MapThread},
{state} = NDSolveProcessEquations[sys, u, {x, -L, L}, {t, 0, tmax},
StepMonitor :> (foo = t)];
MapThread[f_, data_, 1] /; ! FreeQ[f, int] := f @@ data;
NDSolveIterate[state, {0, tmax}];
sol = NDSolveProcessSolutions[state]
]] // AbsoluteTiming


Plot3D[u[x, t] /. sol, {x, -10., 10.}, {t, 0., 2.}]


With the settings PrecisionGoal -> 4, MaxRecursion -> 9 in the NIntegrate, it takes the same amount of time and does more integrations. Breaking down the NDSolve process is explained in the tutorial Components and Data Structures.

• Very impressive (+1). Thanks. – bbgodfrey Mar 6 '19 at 17:21
• @Michael E2, your method is cool and far outside of my knowledge range about mma. Could you explain a little the reason Method -> {"InterpolationPointsSubdivision", Method... } in NIntegrate. I quickly go through tutorial/NIntegrateIntegrationStrategies, and found the general specification NIntegrate[f[x],{x,a,b,c}, Method->{"PrincipalValue",Method->methodspec,"SingularPointIntegrationRadius"->\[Epsilon]}], which appears to be different from yours. Thank you for your answer. – user55777 Mar 7 '19 at 4:04
• By testing the code as it is, i encounter the usual warning: >NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. I don't understand the issue because the singular point has been specified. Please give me suggestion. Thanks a lot. – user55777 Mar 7 '19 at 4:34
• @user55777 Some posts concerning "InterpolationPointsSubdivision" shows use-cases, but it's not well documented. Basically, interpolating functions have (weak) singularities at the interpolation points (nodes), and integration rules are more effective on intervals over which functions do not have singularities. By dividing up the interval of integration at the interpolation points, convergence is more easily obtained. It can be inefficient if there are a lot of interpolation points. – Michael E2 Mar 7 '19 at 19:31
• Technically, it is a preprocessor "strategy" and so can take another Method argument, which may be another "strategy" or "rule." See NIntegrate Introduction, Strategies, Rules, and Preprocessors – Michael E2 Mar 7 '19 at 19:35

After studying during these days, now I could answer the question myself. I admit that both my solution and code are far from good and efficient, even some mistakes or making an unnecessary move. Please give your suggestion if you see anything.

We first create $$2M$$ equidistant grid points $$x_m=(m-M)h$$ with $$m=1,2,...,2M$$. The x-position of the grid points is stored in xtab:

M = 40; L = 10; h = L/M;
xtab = Table[(m - M) h, {m, 1, 2*M}];


Then we should discretize the solution of PDE along $$x$$ into $$2M$$ solutions of a set of coupled ODEs. u[m][t] denotes the solution of function $$u(x,t)$$ at point $$x_m$$. Here, I didn't include the left end-point, since it can be set to be u[0][t]=u[2*M][t] according to the periodicity.

U[t_] = Table[u[m][t], {m, 1, 2*M}];


The spatial derivatives are discretized using 2nd-order central differences, here the periodic condition should be applied. Because I didn't know how to generate these derivatives using ListCorrelate both for boundary points and internal points in one-line command, I manually add the derivatives near the boundary. Please give me some advice if you know how to do this.

1st-derivative wrt x:

internaldUdx = ListCorrelate[{-1, 0, 1}/(2 h), U[t]]; (* for 2<= m <= 19*)
dUdx = Join[{(u[2][t] - u[2*M][t])/(2 h)},
internaldUdx, {(u[1][t] - u[2*M - 1][t])/(2 h)}];


2nd-derivative wrt x:

internaldUdxx = ListCorrelate[{1, -2, 1}/h^2, U[t]]; (* for 2<= m<=19 *)
dUdxx = Join[{(u[2][t] - 2*u[1][t] + u[2*M][t])/h^2},
internaldUdxx, {(u[1][t] - 2 u[2*M][t] + u[2*M - 1][t])/h^2}];


3rd-derivative wrt x

internaldUdxxx = ListCorrelate[{-1, 2, 0, -2, 1}/(2 h^3), U[t]]; (*for 3<= m <= 2*M-2*)
dUdxxx = Join[{(-u[2 M - 1][t] + 2 u[2 M][t] - 2 u[2][t] + u[3][t])/(2 h^3), (-u[2*M][t] + 2 u[1][t] - 2 u[3][t] + u[4][t])/(2 h^3)},
internaldUdxxx, {(-u[2*M - 1 - 2][t] + 2*u[2*M - 1 - 1][t] - 2*u[2*M - 1 + 1][t] + u[1][t])/(2 h^3), (-u[2*M - 2][t] + 2 u[2*M - 1][t] - 2 u[1][t] + u[2][t])/(2 h^3)}];


4th-derivative wrt x:

internaldUdxxxx = ListCorrelate[{1, -4, 6, -4, 1}/h^4, U[t]]; (*for 3 <= m <= 2M-2*)
dUdxxxx = Join[{(u[2*M - 1][t] - 4*u[2*M][t] + 6*u[1][t] - 4*u[1 + 1][t] +
u[1 + 2][t])/h^4, (u[2*M][t] - 4*u[1][t] + 6*u[2][t] - 4*u[2 + 1][t] + u[2 + 2][t])/h^4},
internaldUdxxxx, {(u[2*M - 3][t] - 4*u[2*M - 2][t] + 6*u[2*M - 1][t] - 4*u[2*M][t] + u[1][t])/h^4, (u[2*M - 2][t] - 4*u[2*M - 1][t] + 6*u[2*M][t] - 4 u[1][t] + u[2][t])/h^4}];


To discretize the integral, we may introduce the mid-points: $$x_{m+1/2}=(x_m+x_{m+1})/2$$ for $$m=1,2,....,2M-1$$ with $$x_{1/2}=(-L+x_1)/2$$.

midxtab = Join[{(-L + (1 - M) h)/2}, Table[((m - M) h + (m + 1 - M) h)/2, {m, 1, 2*M - 1}]];
int[midP_] := h/(2 L)*dUdxxxIntP.Cot[\[Pi]*(midxtab[[midP]] - xtab)/(2*L)]


Constructing the system of ODEs and the discrete initial condition:

eqns = Thread[D[U[t], t] == -U[t]*dUdx - dUdxx - dUdxxxx -
Join[Table[1/2*(int[midP] + int[midP + 1]), {midP, 1, 2*M - 1}], {int[2*M] + int[1]}]];



The original PDE now can be solved numerically:

tmax = 10;

lines = NDSolveValue[{eqns, initc}, U[t], {t, 0, tmax},
Method -> {"EquationSimplification" -> "Solve"}] // Flatten;


Then, we can plot by interpolating (appreciating @bbgodfrey's answer to a related question)

surf = Flatten[Table[{(line - M)*h, t, lines[[line]]}, {line, 1, 2*M}, {t, 0,
tmax, 0.2}], 1];

ListPlot3D[surf, PlotRange -> All, AxesLabel -> {"x", "t", "u"}, ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

• You did not use PrincipalValue, although you indicated that you should use. – Alex Trounev Mar 6 '19 at 13:11
• Parameter is not defined dUdxxxIntP. – Alex Trounev Mar 6 '19 at 16:13

We can use iterations. The code is simple but takes time.

L = 10; tmax = 2; del = 10^-6; dx = (L - del)/6 - del;
n = 5;
int[0][x_, t_] := 0
Do[U[i] =
NDSolveValue[{D[u[x, t], t] + u[x, t]*D[u[x, t], x] +
D[u[x, t], {x, 2}] + D[u[x, t], {x, 4}] +
1/(2 L)*int[i - 1][x, t] == 0, u[-L, t] == u[L, t],
u[x, 0] == 0.1*Cos[\[Pi]/L*x]}, u, {x, -L, L}, {t, 0, tmax},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 137}}}];
int[i] = Interpolation[
Flatten[ParallelTable[{{x, t},
NIntegrate[
Derivative[3, 0][U[i]][xp, t]*
Cot[\[Pi] (x - xp)/(2*L)], {xp, -L, x, L},
Method -> "PrincipalValue"]}, {x, -L + del, L - del, dx}, {t,
0, tmax, .2*tmax}], 1]];, {i, 1, n}]

Table[Plot3D[U[i][x, t], {x, -L, L}, {t, 0, tmax}], {i, 1, n}]
Table[Plot3D[
int[i][x, t] - int[i - 1][x, t], {x, -L, L}, {t, 0, tmax},
PlotRange -> All], {i, 1, n}]


• thank you. I found that with the same conditions, your code is much faster than Michael's. Could I consult you three questions: 1. what determine the choice of n = 5? I consider it as an iterator for the calculation of integral, but didn't see a convergence criterion; 2. you actually use FD in x?; 3. Does Table[Plot3D[U[i][x, t], {x, -L, L}, {t, 0, tmax}], {i, 1, n}] only give the intermediate solutions. But how can i plot a convergent u[x,t] with something like Plot3D[u[x, t]/.finalsol, {x, -L, L}, {t, 0, tmax}]? – user55777 Mar 7 '19 at 14:58
• U[i][x, t] this is a solution to the equation at each iteration i. The final solution of the equation U[n][x, t] depends on the selected error. – Alex Trounev Mar 7 '19 at 16:17
• Yes I agree, but it seems that you have not specified an error or convergence criterion? And I think del should not be such a selected error, which is used for avoiding singularity? Am i right? – user55777 Mar 8 '19 at 11:30
• In this example, only convergence is shown. You can add a line of code outside the loop using int[n][x, t] to calculate finalsol. – Alex Trounev Mar 8 '19 at 12:00