# Scipy OdeInt solver with Neumann boundary conditions

I'm using scipy.odeint to solve Fisher-Kolmogorov equation:

$$u_t = u_{xx}+u(1-u)$$

The code can be found here.

From Ablowitz and Zeppetella we know that the analytical solution reads:

$$u(x,t)=\frac{1}{1+e^{{(-\frac{5}{6} t+ \frac{\sqrt{6}x}{6}})^2}}$$

Now, if we compare the analytical solution to the numerical solution obtained with ODE integrator we get the following:

Red: Analytical, Blue: OdeInt

By default scipy.odeint uses Dirichlet boundary conditions.

I would like to know the following:

1) What kind of boundary conditions should I use for odeInt to compare to the analytical solution?

2) How should I include these boundary conditions using scipy.odeint?

• scipy.odeint only handles the time-stepping for the discretized in space system $u_t = A(t,u)$ -- you need to incorporate the boundary condition in the definition of $A$ which you pass to odeint. – Christian Clason Jul 23 '15 at 9:50
• Sorry but I don't see where are the dirichlet bc included. odeint gets a function which computes the rhs and no values for bc are passed. – ilciavo Aug 3 '15 at 2:28
• They are included in the definition of the right-hand side function (basically, the value this function takes in the components corresponding to the boundary points). You should really read up on the basics of finite-difference discretization (e.g., in LeVeque's book Geoff linked to). – Christian Clason Aug 3 '15 at 2:42

From Ablowitz and Zeppetella we know that the analytical solution reads:

$$u(x,t)=\frac{1}{1+e^{{(-\frac{5}{6} t+ \frac{\sqrt{6}x}{6}})^2}}$$

Usually, analytical solutions require the imposition of boundary conditions. Are there any used in the derivation of this expression? If so, you must use the same boundary conditions when implementing the numerical solution intended to replicate this analytical solution.

By default scipy.odeint uses Dirichlet boundary conditions.

I think there might be a misunderstanding here. After semi-discretizing a PDE in space, boundary conditions are enforced by discretization, not by the ODE solver. The ODE solver requires initial conditions, not boundary conditions.

I would like to know the following:

1) What kind of boundary conditions should I use for odeint to compare to the analytical solution?

You should use the same boundary conditions that are used to derive the analytical solution. (See above.)

2) How should I include these boundary conditions using scipy.odeint?

You should include these boundary conditions via your discretization. For examples of how to do so, see LeVeque's finite difference book (book web site), Strikwerda's finite difference book, etc. A previous SciComp.StackExchange question covers the basics:

• Dirichlet boundary conditions amount to setting the right-hand side of equations denoting boundary points to zero, and setting initial conditions on these equations to match a discrete version of the boundary conditions in the continuum problem.
• Neumann boundary conditions amount to replacing the ODEs governing degrees of freedom on those boundaries with (possibly) modified finite difference stencils to approximate the derivative. These modified stencils could be one-sided, or they could be the same stencil as on the interior, but using ghost points. You should consult a textbook (such as LeVeque or Strikwerda) for details.
• Thanks for the recommendation, the analytical solution was found for u=0 at x +/-\Infty .... So I guess I can use Dirichlet boundary conditions right? I'm also following LeVeque's book, and I have my own solver that I want to compare to scipy.odeint results. – ilciavo Jul 24 '15 at 14:31
• Yes, you should be using Dirichlet boundary conditions. However, you will only be able to solve over a finite domain, so you will either have to use the analytical solution to calculate compatible initial conditions on a finite domain or employ a coordinate transformation to map bijectively the infinite domain to a finite domain, then solve on the finite domain. – Geoff Oxberry Jul 24 '15 at 15:32

As for 1), if you want to verify convergence order etc., I'd start with prescribing the value of the exact solution at the boundary as a time-dependent Dirichlet BC. That is, at boundary $x_{\rm left}$, use $u(x_{\rm left},t)$ as Dirichlet value. Be careful to select $t$ correctly to match the time step you are computing. This way, you can verify the convergence order of your interior discretisation without having to worry about the BC.

Unfortunately, as for 2), I do not know how to do this in scipy.odeint so I can't help you there.