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.

  • $\begingroup$ In the expression $\lambda=B(u_0-\left<u\right>)$, is $B$ a constant, or some unknown function? $\endgroup$
    – rajb245
    Commented May 4, 2013 at 19:10
  • $\begingroup$ Also, what is $u_0$ in that expression? $\endgroup$
    – rajb245
    Commented May 4, 2013 at 19:10
  • $\begingroup$ @rajb245 I hope I clarified your questions with my edits. Also see the reference for the paper where this is from -- I might still miss something in my recount of the problem. $\endgroup$
    – Name
    Commented May 4, 2013 at 20:21
  • $\begingroup$ Thanks for the clarifications. Given that the authors used numerical methods to solve the above equation, my guess is that it does not have a closed form solution, so Laplace transform methods will probably not be useful, but I'll take a look into it. $\endgroup$
    – rajb245
    Commented May 5, 2013 at 16:59
  • $\begingroup$ Cheers @rajb245 ... if you happen to know how to recast this problem into something that standard numerical methods can handle I would also be very grateful! Also talked to colleagues of mine about the problem and they suggested (i) imposing the integral condition with a Lagrange multiplier (ii) solving the PDE with the integral as new variable $(u_0 - <u>) = u_i(t)$ so that you get a solution in terms of an integral equation $u(x,t) = h(u_i(t), x, t)$, then solve this integral equation $\endgroup$
    – Name
    Commented May 5, 2013 at 19:13

1 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

pde15s plot

, 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.

  • $\begingroup$ Absolutely fantastic answer. I agree that Runge-Kutta and one of the numerical integration methods you mentioned are most likely the way to go here. Thank you so much for typing this up! From work I've done on a slightly different system, I would expect that there is a range for $u_0$ that gives interesting behavior (due to Turing instability) both with von Neumann b/c and constant initial conditions. At least, this is what I observe for a very similar system (I think only the kinetics differ) described here: math.umn.edu/~ymori/docs/publications/wavepin.pdf $\endgroup$
    – Name
    Commented May 6, 2013 at 20:24
  • $\begingroup$ But generally, you're right that initial conditions with dependence on $x$ give interesting behavior more easily and over a wider range of parameter values. However, let me take a look at the particular example I asked about and see if I can find that regime with Turing instabilities. $\endgroup$
    – Name
    Commented May 6, 2013 at 20:31
  • $\begingroup$ could you help me please math.stackexchange.com/questions/2340696/… $\endgroup$ Commented Jun 29, 2017 at 20:17

You must log in to answer this question.