How to solve partial integro-differential equation?

Question

Suppose the following partial integro-differential equation for a function $u(x,t)$ with $t\geq0$, $x \in [0,L]$:

$\partial_t u = \partial_{xx} u + f(u,\lambda)$

$\lambda = B\left(u_0 - \int_{x=0}^L {u(x,t)dx} \right) = B \left(u_0 - <u> \right)$ (the latter for brevity).

$B,u_0$ are constants. $B$ is the global interaction sensitivity, $u_0$ can be thought of as a measure of how much mass can go into $u$.

$f$ has nice properties and can be assumed to be polynomial in $u$ and $\lambda$. E.g. take $f$ as

$f(u, \lambda) = -u^3 + u + \lambda$

Assume nice initial and boundary conditions, e.g.

$u(x,0) = u_i = \text{const}$

$\partial_t u(0,t) = \partial_t u(L,t) = 0$

(this problem is from: http://link.aip.org/link/doi/10.1063/1.466399)

How can I go about solving this? Is there a way to turn this equation into something more easily solvable?

One approach I've tried was using the Laplace Transform with respect to $x$. Here's my attempt doing this but I am far from avid with this technique so please correct me:

$\mathcal{L}[\partial_t u](s) = \int_{0}^\infty{e^{-xs}\partial_t u(x,t) dx} = \partial_t \overline{u}(t | s) = \frac{d \overline{u}(t | s)}{dt}$,

where I read $\overline{u}(t | s)$ as a function of $t$, given $s$ ($s$ is arbitrary but fixed as I understand).

For the spatial derivative I get (using a table)

$\mathcal{L}[\partial_{xx} u](s) = s^2 \overline{u}(t|s) - s u(0,t) - \partial_x u (0,t)$,

and for $\lambda$ I am not sure what to do.

Answer 1

The equation for $u(x,t)$ is $$ \partial_t u = \partial_{xx}u - u^3 + u + B\left(u_0-\frac{1}{L}\int_{0}^{L} u\;dx\right) $$

This is a non-linear integro-differential equation (if not for the $u^3$, it would be linear). Laplace transform techniques don't get you very far with this equation, because the Laplace transform of $u^3$ is an integral in the Laplace domain. So even after transforming, you have an integro-differential equation; it doesn't simplify to an algebraic equation or ODE, as would happen with linear equations.

It won't be simple to develop your own, but numerical solutions are the way to go here. The authors of the paper from which this comes hastily mention that they use spatial discretization to reduce the problem to a system of ODEs in time, which they then solve with Runge-Kutta methods. The equation above is directly amenable to such an approach, but you have to formulate it. There are also a myriad of other techniques for numerical solutions to such equations. A good package I've been using recently is Chebfun for MATLAB, which can solve this using pseudospectral methods.

MATLAB and Chebfun

This is the easy way. The pde15s routine from the Chebfun package seems to solve this problem correctly. For a bump shaped initial condition with $L=4$, $B=10$, $u_0=1/2$, I get plots like

, which has the diffusion kind of properties I would expect from this equation. My numerical experiments tell me that in order to see interesting behavior, you have to have the equation on a ring with periodic boundary conditions (like the authors of that paper), or you need non constant initial conditions. For the constant initial conditions and Neumann boundary conditions like you propose, the solution doesn't even depend on $x$, i.e. the lack of initial $x$ variations can't be "smoothed out" any further by this kind of diffusion equation, so the solution is just a constant in $x$ along the rod that changes in time.

Finite Differences and Runge-Kutta

Another option is to do what the authors of that paper did. The idea is to sample $u(x,t)$ on a set of $N$ equidistant-spaced points (distance $\Delta x$), say $x_n\in[0,L]$, $n\in[1,2\ldots N]$. These are called nodes in a mesh. Let $u_n(t)=u(x_n,t)$ be the value of $u$ at the point $x_n$. You approximate your spatial derivatives and integrals with finite approximations, which in this case boils down to

$$ \partial_{xx} u(x_n,t)\approx \frac{u_{n+1}(t)-2u_{n}(t)+u_{n-1}(t)}{\left(\Delta x\right)^2} $$ $$ \int_{0}^{L} u\;dx\approx \sum_{n=1}^N u_n(t) \Delta x $$ The integral could be approximated with a higher order rule like the trapezoid or Simpson's rule. There is probably some optimum order for this problem (something like matching the approximation order used for the derivative), you'll have to find the details on how that works. Plugging these into your equation and enforcing the equation on the nodes gives a system of $N$ equations

$$ u_n' = \frac{u_{n+1}-2u_{n}+u_{n-1}}{\left(\Delta x\right)^2} - u_n^3 + u_n + B\left(u_0-\sum_{n=1}^N u_n \Delta x\right),\;\;\forall n\in[1,2\ldots N] $$

I've dropped the explicit time dependence and let the prime denote a time derivative. These are $N$ ODEs in the $N$ unknowns, the $u_n(t)$ functions. You'll also need some more boundary conditions to specify what $u_{0}(t)$ and $u_{N+1}(t)$ are. Now you can use the Runge-Kutta method to solve this system of ODEs at discrete time steps. This is what the authors of the paper did.

You also have the option to directly discretize in space and time, which I think is the more standard approach. This is essentially a heat/diffusion type equation, so the discretization methods used there all apply. Check out the heat-equation example of finite differences on Wikipedia.