I try to solve an integro-differential equation using 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, 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
>NIntegrate::inumr: "The integrand ... has evaluated to non-numerical values for all sampling points"

The warning is understandable because the function `u[xp, t]` is still unknow when `NIntegrate` is evaluated.

My problem may closely relate to [this][1], which is about a linear equation and more like a wave equation in the sense of $u_{tt}$. But my pde is nonlinear and has 1st-derivative in $t$. So I post my problem here for a more general solution. Thank you for any suggestion.

**Update:**

1. As indicated by bbgodfrey in a comment: it is necessary to convert the convolution to the same range as used in `NDSolve`. For a spatial periodic solution, I can change the convolution with integral limits of infinite to a `NIntegrate` over a finite interval.

2. Note it should use `PrincipalValue -> True` here in `NIntegrate`, because there is a singularity at $x=xp$;

Here are my thoughts based on `MethodOfLines` introduced [here][2]

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

    M = 50; dx = 2*L/(2*M); xtab = Table[(m - M) dx, {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, for example:

    dUdxx=ListCorrelate[{1, -2, 1}/dx^2, U[t]]

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$, which are saved in `midxtab`. Correspondingly, we define `u[m+1/2][t]=(u[m][t]+u[m+1][t])/2`.

    midxtab = Join[{(-L+(1-M) dx)/2}, Table[((m - M) dx + (m + 1 - M) dx)/2, {m, 1, 2*M-1}]];

and the corresponding discreted solution at the midpoints

    midU[t_] = Table[Subscript[u, i + 1/2][t], {i, 0, 2*M - 1}];

The value of the midpoint-solution `midU[t_]` can be evaluated to be
`(u[m][t]+u[m+1][t])/2`.    

Now, we discretized the integral at these mid-points, here `dUdxxx[t]` represents 3rd-derivative wrt `x`:

    usum = dUdxxx[t].Cot[\[Pi](midxtab - xtab)/(2*L)]*h;

After `Thread` the system of ODEs (? problem here), and construct the discreted initial condition

    initc = Thread[U[0] == 0.1*Cos[\[Pi]/L*xtab]];

The original PDE could be solved numerically

    lines = NDSolveValue[{eqns, initc}, U[t], {t, 0, tmax}] // Flatten;


  [1]: https://mathematica.stackexchange.com/questions/107538/solving-pde-involving-hilbert-transform-numerically
  [2]: https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html