Boundary value problems#
KEYWORDS: scipy.integrate.solve_bvp, numpy.polyfit
Solving nonlinear BVPs by finite differences#
Adapted from Example 8.7 in Numerical Methods in Engineering with Python by Jaan Kiusalaas.
We want to solve \(y''(x) = -3 y(x) y'(x)\) with \(y(0) = 0\) and \(y(2) = 1\).
This is a boundary value problem not an initial value problem. First we consider using a finite difference method. We discretize the region and approximate the derivatives as:
\(y''(x) \approx \frac{y_{i-1} - 2 y_i + y_{i+1}}{h^2} \)
\(y'(x) \approx \frac{y_{i+1} - y_{i-1}}{2 h} \)
We define a function \(y''(x) = F(x, y, y')\). At each node in our discretized region, we will have an equation that looks like \(y''(x) - F(x, y, y') = 0\), which will be nonlinear in the unknown solution \(y\). The set of equations to solve is:
\begin{eqnarray} y_0 - \alpha &=& 0 \ \frac{y_{i-1} - 2 y_i + y_{i+1}}{h^2} + (3 y_i) (\frac{y_{i+1} - y_{i-1}}{2 h}) &=& 0 \ y_L - \beta &=&0 \end{eqnarray}
Since we use a nonlinear solver, we will have to provide an initial guess to the solution. We will in this case assume a line. In other cases, a bad initial guess may lead to no solution.
import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
x1 = 0.0
x2 = 2.0
alpha = 0.0
beta = 1.0
We need to specify a grid of points to discretize the solution on. We will start with a small grid because it is easy to visualize, but note that the grid spacing determines how good the approximation to the derivative is, so we will have to return here to see what the impact of our spacing is.
N = 100
X, h = np.linspace(x1, x2, N, retstep=True)
X, h
(array([0. , 0.02020202, 0.04040404, 0.06060606, 0.08080808,
0.1010101 , 0.12121212, 0.14141414, 0.16161616, 0.18181818,
0.2020202 , 0.22222222, 0.24242424, 0.26262626, 0.28282828,
0.3030303 , 0.32323232, 0.34343434, 0.36363636, 0.38383838,
0.4040404 , 0.42424242, 0.44444444, 0.46464646, 0.48484848,
0.50505051, 0.52525253, 0.54545455, 0.56565657, 0.58585859,
0.60606061, 0.62626263, 0.64646465, 0.66666667, 0.68686869,
0.70707071, 0.72727273, 0.74747475, 0.76767677, 0.78787879,
0.80808081, 0.82828283, 0.84848485, 0.86868687, 0.88888889,
0.90909091, 0.92929293, 0.94949495, 0.96969697, 0.98989899,
1.01010101, 1.03030303, 1.05050505, 1.07070707, 1.09090909,
1.11111111, 1.13131313, 1.15151515, 1.17171717, 1.19191919,
1.21212121, 1.23232323, 1.25252525, 1.27272727, 1.29292929,
1.31313131, 1.33333333, 1.35353535, 1.37373737, 1.39393939,
1.41414141, 1.43434343, 1.45454545, 1.47474747, 1.49494949,
1.51515152, 1.53535354, 1.55555556, 1.57575758, 1.5959596 ,
1.61616162, 1.63636364, 1.65656566, 1.67676768, 1.6969697 ,
1.71717172, 1.73737374, 1.75757576, 1.77777778, 1.7979798 ,
1.81818182, 1.83838384, 1.85858586, 1.87878788, 1.8989899 ,
1.91919192, 1.93939394, 1.95959596, 1.97979798, 2. ]),
np.float64(0.020202020202020204))
Now, we can define functions for the differential equation, and for the nonlinear equations.
def residuals(y):
"When we have the right values of y, this function will be zero."
res = np.zeros(y.shape) # we need a zero for each node
res[0] = y[0] - alpha # this is the boundary value y(alpha) = 0
for i in range(1, N - 1):
x = X[i] # This is not actually used
# Approximation of y'' from the current point
YPP = (y[i - 1] - 2 * y[i] + y[i + 1]) / h**2
# Approximation of y'
YP = (y[i + 1] - y[i - 1]) / (2 * h)
# y'' + 3 * y * y' = 0
res[i] = YPP + 3 * y[i] * YP
res[-1] = y[-1] - beta # y(beta) = 0
return res
We need a guess, and here we guess a line. It is always a good idea to plot your guess just to make sure it does what you want. Here, we want it to have the right boundary values.
X
array([0. , 0.02020202, 0.04040404, 0.06060606, 0.08080808,
0.1010101 , 0.12121212, 0.14141414, 0.16161616, 0.18181818,
0.2020202 , 0.22222222, 0.24242424, 0.26262626, 0.28282828,
0.3030303 , 0.32323232, 0.34343434, 0.36363636, 0.38383838,
0.4040404 , 0.42424242, 0.44444444, 0.46464646, 0.48484848,
0.50505051, 0.52525253, 0.54545455, 0.56565657, 0.58585859,
0.60606061, 0.62626263, 0.64646465, 0.66666667, 0.68686869,
0.70707071, 0.72727273, 0.74747475, 0.76767677, 0.78787879,
0.80808081, 0.82828283, 0.84848485, 0.86868687, 0.88888889,
0.90909091, 0.92929293, 0.94949495, 0.96969697, 0.98989899,
1.01010101, 1.03030303, 1.05050505, 1.07070707, 1.09090909,
1.11111111, 1.13131313, 1.15151515, 1.17171717, 1.19191919,
1.21212121, 1.23232323, 1.25252525, 1.27272727, 1.29292929,
1.31313131, 1.33333333, 1.35353535, 1.37373737, 1.39393939,
1.41414141, 1.43434343, 1.45454545, 1.47474747, 1.49494949,
1.51515152, 1.53535354, 1.55555556, 1.57575758, 1.5959596 ,
1.61616162, 1.63636364, 1.65656566, 1.67676768, 1.6969697 ,
1.71717172, 1.73737374, 1.75757576, 1.77777778, 1.7979798 ,
1.81818182, 1.83838384, 1.85858586, 1.87878788, 1.8989899 ,
1.91919192, 1.93939394, 1.95959596, 1.97979798, 2. ])
# we need an initial guess
init = alpha + (beta - alpha) / (x2 - x1) * X
plt.plot(X, init);
We should check our residuals function. We mostly want to see that it runs, and produces the right shaped output.
len(residuals(init))
100
Now, we solve the BVP.
Y, info, status, msg = fsolve(residuals, init, full_output=1)
print(msg)
plt.plot(X, Y)
plt.xlabel('x')
plt.ylabel('y');
The solution converged.
np.allclose(residuals(Y), 0.0, atol=1e-7), np.mean(np.abs(residuals(Y)))
(True, np.float64(1.5879661186615786e-08))
np.min(residuals(Y)), np.max(residuals(Y))
(np.float64(-5.32498396488279e-08), np.float64(4.9007423896796354e-08))
The solution is has some apparent discontinuities because we only used about 10 points. How can you tell if the solution is correct? We can estimate the derivatives, and see how well they fit the equation. We look for:
\(y'' + 3 y y' = 0\) for all \(x\).
yp = np.gradient(Y, X, edge_order=2)
ypp = np.gradient(yp, X, edge_order=2)
plt.plot(X, ypp + 3 * Y * yp)
plt.xlabel('x')
plt.ylabel('residuals');
This result doesn’t look great at the origin, but remember:
we used a coarse grid, so the derivative approximations are probably not that accurate
Numerical derivatives at the end-points are less accurate than in the middle.
exercise Go back and repeat this for a finer grid, e.g. with 50, 100 points.
The approach described here is pretty general. Here, we were able to solve a second-order BVP by discretizing it, approximating the derivatives at the points, and solving the corresponding nonlinear algebra equations. This approach can be extended in a variety of ways, including to systems of equations, and to 2D or 3D systems (where this approach is called finite-element). You will see these kinds of problems extensively in the spring semster in the Transport class.
As we have seen before, however, there are functions in scipy that can help solve these problems.
Introduction to solve_bvp#
from scipy.integrate import solve_bvp
solve_bvp?
A worked bvp problem#
In the pressure driven flow of a fluid with viscosity \(\mu\) between two stationary plates separated by distance \(d\) and driven by a pressure drop \(\Delta P/\Delta x\), the governing equations on the velocity \(u\) of the fluid are (assuming flow in the x-direction with the velocity varying only in the y-direction):
with boundary conditions \(u(y=0) = 0\) and \(u(y=d) = 0\), i.e. the no-slip condition at the edges of the plate.
we convert this second order BVP to a system of ODEs by letting \(u_1 = u\), \(u_2 = u_1'\) and then \(u_2' = u_1''\). This leads to:
\(\frac{d u_1}{dy} = u_2\)
\(\frac{d u_2}{dy} = \frac{1}{\mu}\frac{\Delta P}{\Delta x}\)
with boundary conditions \(u_1(y=0) = 0\) and \(u_1(y=d) = 0\).
for this problem we let the plate separation be d=0.1, the viscosity \(\mu = 1\), and \(\frac{\Delta P}{\Delta x} = -100\).
import numpy as np
d = 0.1
mu = 1
deltaPdeltax = -100
The function defining the BVP has to return an array that has a row for each equation, and a column for each value in the grid.
def bvp(y, U):
u1, u2 = U
du1dy = u2 # this will be an array of values of u2(y)
du2dy = np.ones(y.shape) / mu * deltaPdeltax
# Both rows must be arrays.
return [du1dy, du2dy]
The boundary condition function will get the whole numeric solution at each boundary. We want \(u1(a) = 0\) and \(u1(b)=0\).
def bc(Ua, Ub):
u1a, u2a = Ua
u1b, u2b = Ub
return [u1a, u1b] # These will be zero at the solution
Next, we need an initial guess for u1 and u2 on a grid of points. You have to make some decisions here. You need a guess that is reasonably close, but not hard to construct. Here, we anticipate a solution that looks parabolic, and that goes through the points: (0, 0), (d, 0), and some point at (d / 2, ?), where ? represents the point of maximum velocity in middle. We can easily get this polynomial with np.polyfit.
We don’t know what the maximum velocity is, so we make a guess, say 0.5. Then, we get the parameters, and apply them to an array of y values.
pars = np.polyfit([0, d / 2, d], # x-points
[0, 0.5, 0], # y-points
2) # polynomial order
pars # parabola coefficients for U1
array([-2.00000000e+02, 2.00000000e+01, -3.84592537e-16])
Now, we can define a Y grid and define the guess for the first U1.
Y = np.linspace(0, d)
U1 = np.polyval(pars, Y)
plt.plot(Y, U1)
plt.xlabel('Y')
plt.ylabel('U1')
Y.shape, U1.shape
((50,), (50,))
We also need a guess for U2, and in this case we know that \(u2 = u1'\), so we just use that.
du1 = np.polyder(pars)
U2 = np.polyval(du1, Y)
plt.plot(Y, U2)
#U2 = np.gradient(U1, Y, edge_order=2)
# Combine the two rows into one initial guess for U
U = np.array([U1, U2])
print(U.shape)
(2, 50)
You should always visualize the guess to make sure it does what you want. It is hard to make these!
import matplotlib.pyplot as plt
plt.plot(Y, U[0], label='u1')
plt.gca().tick_params('y', colors='b')
plt.ylabel('u1')
plt.twinx() # makes graphs share the x-axis, and adds a second y-axis
plt.plot(Y, U[1], 'r', label='u2')
plt.gca().tick_params('y', colors='r')
plt.ylabel('u2');
Now, we are ready to solve the BVP.
from scipy.integrate import solve_bvp
sol = solve_bvp(bvp, bc, Y, U)
print(sol.message)
plt.plot(sol.x, sol.y[0])
plt.xlabel('y')
plt.ylabel('U');
The algorithm converged to the desired accuracy.
exercise Try using different guesses, e.g. lines, or triangle shapes, etc. What else looks like this shape? Half a cycle of a sin wave? A semi-circle?
exercise How can you show this is a solution to the BVP?
Concentration profile in a particle#
Another typical boundary value problem in chemical engineering is the concentration profile inside a catalyst particle. Here is the dimensionless equation for a second order reaction in a slab. Note here we have a boundary condition on the derivative at the origin. This kind of condition means either there is no flux at this position, or that the slab is symmetric about this position.
\(\frac{d^2c}{dx^2} = \Phi^2 c^2\)
with \(c'(0)\) = 0 and \(c(1) = 1\)
We again convert this to a system of first order differential equations like this:
Let c1 = c, c1’ = c’, and c2 = c1’, so c2’ = c1’’ = c’’
Then we have:
\(c1' = c2\)
\(c2' = \Phi^2 c1^2\)
with boundary conditions \(c1'(0) = 0\) and \(c1(1) = 1\).
We begin with the required functions:
Phi = 100 # Constant given for the problem
def bvp(x, C):
c1, c2 = C # these are two rows for all the values of x
dc1dx = c2
dc2dx = Phi**2 * c1**2
return [dc1dx, dc2dx]
def bc(Ca, Cb):
c1a, c2a = Ca
c1b, c2b = Cb
# Now, evaluate the derivatives at the first boundary condition
c1prime, c2prime = bvp(0, [c1a, c2a])
return [c1prime, # will all equal zero. Alternatively use c2a
c1b - 1] # c1(b) = 1
We need an initial guess. We make a naive one, that \(c(x) = 1\) in the slab, i.e. there is no reaction. As usual, we visualize the guess to be sure it does what we intended.
X = np.linspace(0, 1)
C1 = X**1200#np.ones(X.shape)
C2 = np.gradient(C1, X)
plt.plot(X, C1);
Now we solve the system.
C = [C1, C2]
sol = solve_bvp(bvp, bc, X, C)
print(sol.message)
len(sol.x)
The algorithm converged to the desired accuracy.
227
plt.plot(sol.x, sol.y[0])
plt.xlabel('x')
plt.ylabel('C')
plt.xlim([0, 1])
plt.ylim([0, 1]);
plt.plot(X, X**8) # Alternative initial guess
len(X)
50
c = sol.y[0]
len(c)
227
You can see the solution looks nothing like our initial guess. In this case, a high thiele modulus means most of the reaction happens near the catalyst surface, and the interior of the slab has hardly any reactant in it. This solution is consistent with that.
The effectiveness factor for this system is defined by:
\(E = \int_0^1 c^2 dx\)
We can estimate this with the trapezoid or Simpson’s method (remember that the solution is a vector of numbers).
c = sol.y[0]
print(np.trapz(c**2, sol.x))
from scipy.integrate import simps
print(simps(c**2, sol.x))
0.008249311811255714
/tmp/ipykernel_2389/1627362779.py:2: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
print(np.trapz(c**2, sol.x))
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In[26], line 4
1 c = sol.y[0]
2 print(np.trapz(c**2, sol.x))
----> 4 from scipy.integrate import simps
5 print(simps(c**2, sol.x))
ImportError: cannot import name 'simps' from 'scipy.integrate' (/opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/scipy/integrate/__init__.py)
Or, we can use the dense_output of the solution with quad.
from scipy.integrate import quad
def integrand(x):
c1, c2 = sol.sol(x)
return c1**2
quad(integrand, 0, 1)
(0.008164991938713655, 1.1791754240830934e-08)
excercise Repeat this example for different values of Φ.
exercise Try different kinds of guesses. Think of a guess that has the properties of the boundary conditions, e.g. c’(0) = 0, and c(1) = 1.
exercise Evaluate the quality of the solution based on the equations.
Summary#
Today, we leveraged the ability to solve systems of nonlinear algebraic equations to solve boundary value problems by discretizing them on a grid, approximating them at the grid points, and then solving the resulting nonlinear equations.
We also learned about the solve_bvp function, which is in scipy.integrate to solve systems of first-order boundary value problems.
Next time, we will return to nonlinear algebra to see how the algorithms can be used to find minima and maxima.