Using the BVP solver in scipy

AMath 585, Winter Quarter 2020 at the University of Washington. Developed by R.J. LeVeque and distributed under the BSD license. You are free to modify and use as you please, with attribution.

These notebooks are all available on Github.

It is valuable to understand how to go about setting up a BVP solver yourself in order to understand the various issues that arise, and the ways that things might go wrong. But for many real problems the best way to actually solve it is to use high-quality software that has been designed by experts and that has features such as adaptive choice of grids built in.

This notebook briefly illustrates how to use the scipy.integrate.solve_bvp, which uses a collocation method and adaptively chooses the grid nodes based on a specified initial grid.

The collocation method used does not return approximate solution values at a set of points (as the finite difference method does). Instead it returns a representation of a cubic sline, a piecewise cubic function that is a cubic polynomial on each interval between nodes. A spline has the property that it is continuous and also has continuous first and second derivatives at each node, and hence is quite smooth. To plot the solution you have to evaluate this function at a set of points.

Specifying the BVP

The solve_bvp function is written in a way that it can be used to solve BVPs with arbitrary order derivatives, e.g. both $u''(x) = f(x)$ and $u''''(x) = f(x)$ (and much more complicated nonlinear problems). It does so by actually solving a system of ODEs of the form $$ y'(x) = f(x,y) $$ where $y(x)$ is a vector with $n$ components. As described at the start of Chapter 5 in the textbook, a higher-order scalar equation can generally be rewritten as a first-order system.

For example, for $u''(x) = f(x)$, let $y_0(x) = u(x)$ and $y_1(x) = u'(x)$. Then the system of equations is \begin{align*} y_0'(x) &= y_1(x),\\ y_1'(x) &= f(x). \end{align*} If the original problem had Dirichlet boundary conditions $u(a)=\alpha$, $u(b)=\beta$ then the system has boundary conditions \begin{align*} y_0(a) &= \alpha,\\ y_0(b) &= \beta. \end{align*} A Neumann boundary condition $u'(a) = \alpha_1$ could be expressed as $y_1(a) = \alpha_1$ in place of the first equation above.

The solve_bvp function has one argument that is the function f and another argument that is a function bc that takes two n-vectors ya, yb as input and returns $n$ values that should be zero, e.g. ya[0]-alpha and yb[0]-beta for the Dirichlet case in this second-order ($n=2$) equation. You also have to specify the initial set of nodes (the first and last of which are the end points of the interval), and an initial guess for the solution values at these nodes. The software will adjust the nodes to try to achieve some estimated error.

There are many other optional arguments that can be specified, see the documentation.

In [1]:
%matplotlib inline
In [2]:
from pylab import *
In [3]:
from scipy.integrate import solve_bvp
from scipy.interpolate import interp1d  # used for continuation

Beam equation

First let's solve $u''''(x) = f(x)$ as in Homework 2, with $$ u(0)=\alpha_0,\quad u'(0)=\alpha_1,\quad u(L)=\beta_0,\quad u'(L)=\beta_1. $$

In this case we introduce $$ y_0(x) = u(x), \quad y_1(x) = u'(x), \quad y_2(x) = u''(x), \quad y_3(x) = u'''(x), $$ and the system of ODEs is then \begin{align*} y_0'(x) &= y_1(x),\\ y_1'(x) &= y_2(x),\\ y_2'(x) &= y_3(x),\\ y_3'(x) &= f(x). \end{align*}

In [4]:
L = 5.
gamma = -0.01
ax = 0.
bx = L
alpha0 = 0.
alpha1 = 0.
beta0 = 0
beta1 = 0

def f(x,y):
    fval = zeros(y.shape)
    fval[0,:] = y[1,:]       # y0' = y1
    fval[1,:] = y[2,:]       # y1' = y2
    fval[2,:] = y[3,:]       # y2' = y3
    fval[3,:] = gamma        # y3' = gamma
    return fval

def bc(ya, yb):
    bcval = array([ya[0]-alpha0,
                   ya[1]-alpha1,
                   yb[0]-beta0,
                   yb[1]-beta1])
    return bcval

Choose an initial set of only 10 nodes and initial guess that is identically zero:

In [5]:
x_initial = linspace(ax, bx, 10)
y_initial = zeros((4,len(x_initial)))

bvp_output = solve_bvp(f, bc, x_initial, y_initial,
                       max_nodes=1000, tol=1e-3)

print('niter = ', bvp_output.niter)
print('status = ', bvp_output.status)
print('message = ', bvp_output.message)
if not bvp_output.success:
    print('*** Warning, did not succeed!')

niter =  1
status =  0
message =  The algorithm converged to the desired accuracy.

Note that it says convergence was obtained in 1 iteration. For a general nonlinear problem a "damped Newton" method is used and for this linear problem it converges in 1 iteration from any initial guess.

Evaluating the solution

Note that bvp_output.sol is the solution, as a function that could be evaluated at any point in the interval:

In [6]:
bvp_output.sol
Out[6]:
<scipy.interpolate.interpolate.PPoly at 0x119eb0888>
In [7]:
bvp_output.sol(2.5)
Out[7]:
array([-1.62735609e-02,  3.93023288e-19,  1.04166667e-02, -8.19174975e-19])

The first component of this is the approximation to $y_0(2.5) = u(2.5)$.
Note that at this point, the midpoint of the interval, the first and third derivatives are close to zero.

Plotting the solution

If there are lots of nodes you could simply evaluate the function at these nodes and plot the resulting points. But recall the plot function connects points with linear segments, so in this case we get something that looks piecewise linear:

In [8]:
x = bvp_output.x
y = bvp_output.sol(x)  # evaluate the sol function returned at the x points
u = y[0,:]  # select the y_0 = u component
plot(x, u, 'rx-')
title('Solution returned with %i nodes' % len(x));

This might not look very accurate, but here we have only plotted the solution at the 10 nodes in the solution returned. If we evaluate the function returned on a finer grid we see that it agrees well with the expected exact solution:

In [9]:
utrue = lambda x: gamma/24. * (x**4 - 2*L * x**3 + L**2 * x**2)
xfine = linspace(ax, bx, 101)

figure(figsize=(12,6))
plot(x, u, 'r+', markersize=10, label='nodes') # plot nodes

# evaluate the numerical solution on the xfine grid:
y = bvp_output.sol(xfine)
ufine = y[0,:]  # select the y_0 = u component
plot(xfine, ufine, 'r-', linewidth=2, label='numerical')

# plot the true solution on the xfine grid:
plot(xfine, utrue(xfine), 'b--', label='exact')
legend()

# compute the max error at these xfine points:
err = abs(ufine - utrue(xfine)).max()
print('Maximum error on fine grid is %g' % err)

Maximum error on fine grid is 2.48073e-06

Cantilever beam

We can easily change the boundary conditions to model a cantilever beam, that is rigidly supported at the left boundary but unsupported at the right (e.g. a beam supporting a balcony sticking out from a building). In this case the proper BCs at the free boundary are $u''(L) = u'''(L) = 0$, which become $y_2(L)=y_3(L)=0$ in the system.

In [10]:
def bc(ya, yb):
    bcval = array([ya[0]-alpha0,
                   ya[1]-alpha1,
                   yb[2],
                   yb[3]])
    return bcval

x_initial = linspace(ax, bx, 10)
y_initial = zeros((4,len(x_initial)))

bvp_output = solve_bvp(f, bc, x_initial, y_initial,
                       max_nodes=1000, tol=1e-3)

print('niter = ', bvp_output.niter)
print('status = ', bvp_output.status)
print('message = ', bvp_output.message)
if not bvp_output.success:
    print('*** Warning, did not succeed!')

niter =  1
status =  0
message =  The algorithm converged to the desired accuracy.
In [11]:
x = bvp_output.x
y = bvp_output.sol(x)  # evaluate the sol function returned at the x points
u = y[0,:]  # select the y_0 = u component
plot(x, u, 'rx')
# evaluate the numerical solution on the xfine grid:
y = bvp_output.sol(xfine)
ufine = y[0,:]  # select the y_0 = u component
plot(xfine, ufine, 'r-', linewidth=2, label='numerical')
title('Solution returned with %i nodes' % len(x));

Singular perturbation problem

Solve the nonlinear BVP $$ \epsilon u''(x) + u(x)(u'(x) - 1) = f(x) $$ with Dirichlet boundary conditions, in the case where $\epsilon > 0$ is very small, the singular perturbation problem discussed in Section 2.17 of the textbook.

First try a large epsilon for which the solution is smooth:

In [12]:
ax = 0.; alpha = -1.;
bx = 1.; beta = 1.5; 
epsilon = 0.1
In [13]:
def f(x,y):
    fval = zeros(y.shape)
    fval[0,:] = y[1,:]                          # y0' = y1
    fval[1,:] = -y[0,:]*(y[1,:] - 1) / epsilon  # y1' = -y0(y1-1)/epsilon
    return fval

def bc(ya, yb):
    return array([ya[0]-alpha, yb[0]-beta])
In [14]:
x_initial = linspace(ax, bx, 10)
y0_initial = alpha + x_initial * (beta-alpha) / (bx-ax)
y1_initial = (beta-alpha) / (bx-ax)  * ones(x_initial.shape)
y_initial = vstack((y0_initial, y1_initial))

bvp_output = solve_bvp(f, bc, x_initial, y_initial,
                       max_nodes=10000, tol=1e-3)

print('niter = ', bvp_output.niter)
print('status = ', bvp_output.status)
print('message = ', bvp_output.message)
if not bvp_output.success:
    print('*** Warning, did not succeed!')
    
x = bvp_output.x
y = bvp_output.sol(x)
u = y[0,:]
plot(x, u, 'rx-')
title('Solution returned with %i nodes' % len(x));

niter =  4
status =  0
message =  The algorithm converged to the desired accuracy.
In [15]:
epsilon = 0.005
def f(x,y):
    fval = zeros(y.shape)
    fval[0,:] = y[1,:]                          # y0' = y1
    fval[1,:] = -y[0,:]*(y[1,:] - 1) / epsilon  # y1' = -y0(y1-1)/epsilon
    return fval

    
x_initial = linspace(ax, bx, 500)
y0_initial = alpha + x_initial * (beta-alpha) / (bx-ax)
y1_initial = (beta-alpha) / (bx-ax)  * ones(x_initial.shape)
y_initial = vstack((y0_initial, y1_initial))

bvp_output = solve_bvp(f, bc, x_initial, y_initial,
                       max_nodes=10000, tol=1e-3)

print('niter = ', bvp_output.niter)
print('status = ', bvp_output.status)
print('message = ', bvp_output.message)
if not bvp_output.success:
    print('*** Warning, did not succeed!')

niter =  4
status =  0
message =  The algorithm converged to the desired accuracy.
In [16]:
x = bvp_output.x
y = bvp_output.sol(x)
u = y[0,:]
plot(x, u, 'rx-')
title('Solution returned with %i nodes' % len(x));

Note that in this case it refined the grid substantially. Mostly it added small cells near the interior layer, as we can see by plotting the cell widths:

In [17]:
dx0 = x_initial[1] - x_initial[0]
print('Initial interval length: %g' % dx0)
dx = diff(x)
xmid = 0.5*(x[1:] + x[:-1])
plot(xmid, dx, 'b')
title('Final interval lengths');

Initial interval length: 0.00200401

Convergence failure

Newton's method can fail to converge if the initial guess isn't good enough, or doesn't have enough nodes:

In [18]:
epsilon = 0.005
def f(x,y):
    fval = zeros(y.shape)
    fval[0,:] = y[1,:]                          # y0' = y1
    fval[1,:] = -y[0,:]*(y[1,:] - 1) / epsilon  # y1' = -y0(y1-1)/epsilon
    return fval

    
x_initial = linspace(ax, bx, 100)
y0_initial = alpha + x_initial * (beta-alpha) / (bx-ax)
y1_initial = (beta-alpha) / (bx-ax)  * ones(x_initial.shape)
y_initial = vstack((y0_initial, y1_initial))

bvp_output = solve_bvp(f, bc, x_initial, y_initial,
                       max_nodes=10000, tol=1e-3)

print('niter = ', bvp_output.niter)
print('status = ', bvp_output.status)
print('message = ', bvp_output.message)
if not bvp_output.success:
    print('*** Warning, did not succeed!')

niter =  5
status =  1
message =  The maximum number of mesh nodes is exceeded.
*** Warning, did not succeed!
In [19]:
x = bvp_output.x
y = bvp_output.sol(x)
u = y[0,:]
plot(x, u, 'rx-')
title('Solution returned with %i nodes' % len(x));

A better initial guess

Rather than choosing an initial guess that's just linear, for this problem we know roughly what the solution looks like. Start with a discontinuity at the correct location:

In [20]:
x_initial = linspace(ax, bx, 200)
y0_initial = where(x_initial<0.25, x_initial-1., x_initial+0.5)
y1_initial = ones(x_initial.shape)
y_initial = vstack((y0_initial, y1_initial))

plot(x_initial, y0_initial, 'bx-')
title('Initial guess')
Out[20]:
Text(0.5, 1.0, 'Initial guess')
In [21]:
bvp_output = solve_bvp(f, bc, x_initial, y_initial,
                       max_nodes=10000, tol=1e-3)

print('niter = ', bvp_output.niter)
print('status = ', bvp_output.status)
print('message = ', bvp_output.message)
if not bvp_output.success:
    print('*** Warning, did not succeed!')

niter =  4
status =  0
message =  The algorithm converged to the desired accuracy.
In [22]:
x = bvp_output.x
y = bvp_output.sol(x)
u = y[0,:]
plot(x, u, 'rx-')
title('Solution returned with %i nodes' % len(x));
In [23]:
dx0 = x_initial[1] - x_initial[0]
print('Initial interval length: %g' % dx0)
dx = diff(x)
xmid = 0.5*(x[1:] + x[:-1])
plot(xmid, dx, 'b')
title('Final interval lengths');

Initial interval length: 0.00502513

Note that there are many small intervals near $x=0.25$, clearer if we plot dx vs. index i rather than vs. x[i]:

In [24]:
plot(dx, 'b')
title('Final interval lengths vs. index');