Contents

Tutorial 11 - Boundary value problems

Tutorial 11 - Boundary value problems

Boundary value problems of ordinary differential equations, finite difference method, shooting method, finite element method.

import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg, integrate, optimize

Shooting method

The shooting method has its origin in artillery. When firing a cannon towards a target, the first shot is fired in the general direction of the target. If the cannon ball hits too far to the right, the cannon is pointed a little to the left for the second shot, and vice versa. This way, the cannon balls will hit ever closer to the target. Consider the following boundary value problem,

\[ \tag{1} \label{eq:shooting_bvp} y^{\prime \prime}(x) = f(x, y(x), y^{\prime}(x)), \quad x \in [a, b], \quad y(a) = \alpha, \quad y(b) = \beta. \]

Let \( y(x; \gamma) \) denote the solution of the initial value problem

\[ \tag{2} \label{eq:shooting_ivp} y^{\prime \prime}(x) = f(x, y(x), y^{\prime}(x)), \quad x \in [a, b], \quad y(a) = \alpha, \quad y^{\prime}(a) = \gamma. \]

Define the function \( F(\gamma) \) as the difference between \( y(b; \, \gamma) \) and the specified boundary value \( \beta \),

\[ F(\gamma) = y(b; \, \gamma) - \beta. \]

If \( \gamma^{*} \) is a root of \( F(\gamma) \) then the solution \( y(x; \gamma^{*}) \) of initial value problem \( \eqref{eq:shooting_ivp} \) is also a solution of boundary value problem \( \eqref{eq:shooting_bvp} \). The usual methods for finding roots may be employed here, such as the bisection method or Newton-Raphson method.

Note: The \( 2^{nd} \) order explicit ordinary differential equation

\[ y^{\prime \prime} = f(x, y(x), y^{\prime}(x)), \quad x \in [a, b], \quad y(a) = \alpha_1, \quad y^{\prime}(a) = \alpha_2 \]

can be expressed as a system of two \( 1^{st} \) order differential equations

\[\begin{split} \begin{align} & y^{\prime}(x) = z(x), \quad y(a) = \alpha_1 \\ & z^{\prime}(x) = f(x, y(x), z(x)), \quad z(a) = \alpha_2 \end{align} \end{split}\]

by defining a new variable \( z(x) = y^{\prime}(x) \).

Exercise 11.1: Implement the shooting method for the \( 2^{nd} \) order explicit ordinary differential equation.

def shooting_method(f, g, a, b, n, alpha, beta):
    """
    Solves 2nd order explicit ODE in the form y(x)'' = f(x, y(x), y(x)') in [a, b] with y(a) = alpha, y(b) = beta
    using shooting method.
    Args:
        f (function): function 
        g (function): function
        a (float): The left-most point of the interval
        b (float): The right-most point of the interval
        n (int): number of points dividing interval [a, b]
        alpha (float): initial value at x = a
        beta (float): initial value at x = b
    Returns:
        numpy.ndarray: Vector of grid-points
        numpy.ndarray: Vector of solution
    """

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_bvp function.

try:
    np.testing.assert_array_almost_equal(shooting_method(lambda x, y, z: x, lambda x, y, z: z, \
        0.0, 1.0, 1000, 0.0, 0.0)[1], integrate.solve_bvp(lambda x, y: np.vstack((y[1], x)), \
        lambda a, b: np.array([a[0], b[0]]), np.linspace(0.0, 1.0, 1000), np.zeros((2, 1000))).y[0], decimal=5)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 11.2: Solve the following ordinary differential equation,

\[ y^{\prime \prime} \left( x \right) = \frac{1}{2} \sqrt{1 + \left[ y^{\prime} \left( x \right) \right]^2}, \quad x \in [-1, \, 1], \quad y \left( -1 \right) = y \left( 1 \right) = 0 \]

using the shooting method with a partition of size \( N = 20 \). Plot the solution obtained by the shooting method as well as the analytical solution

\[ y \left( x \right) \approx 2 \cosh \left( \frac{x}{2} \right) - 2.255 \]

in the interval \( x \in [-1, \, 1] \).

# add your code here

Note: The solution of the equation in Exercise 11.2 is a catenary curve.

Finite difference method

The finite difference method is a numerical technique for solving differential equations by approximating derivatives with finite differences. Consider the following boundary value problem,

\[ \tag{3} \label{eq:fdm_bvp} p \left(x \right) y^{\prime \prime} \left( x \right) + q \left(x \right) y^{\prime} \left( x \right) + r \left(x \right) y \left( x \right) + s \left(x \right) = 0, \quad x \in [a, b], \quad y(a) = \alpha, \quad y(b) = \beta. \]

Let

\[ y_i^{\prime} \approx \frac{y_{i + 1} - y_{i}}{h} \quad \mathrm{and} \quad y_i^{\prime \prime} \approx \frac{y_{i + 1} - 2 y_{i} + y_{i - 1}}{h^2}, \quad i = 0, \dots, n - 1, \]

be the finite difference approximations of \( y^{\prime} \left( x_i \right) \) and \( y^{\prime \prime} \left( x_i \right) \), respectively, where \( h > 0 \) is a step size. Boundary value problem \( \eqref{eq:fdm_bvp} \) can be then rewritten as a system of linear algebraic equations

\[ \tag{4} \label{eq:fdm_system} \tilde{p_i} y_{i - 1} + \tilde{q_i} y_i + \tilde{r_i} y_{i + 1} = \tilde{s_i}, \quad y_0 = \alpha, \quad y_{n - 1} = \beta. \]

where \( \tilde{p_i} = p_i \), \( \tilde{q_i} = \left( -2 p_i - q_i h + r_i h^2 \right) \), \( \tilde{r_i} = \left(p_i + q_i h \right) \), and \( \tilde{s_i} = -s_i h^2 \). In the matrix form, system of equations \( \eqref{eq:fdm_system} \) reads as

\[\begin{split} \begin{pmatrix} 1 & & & & \\ \tilde{p}_1 & \tilde{q}_1 & \tilde{r}_1 & & \\ & \ddots & \ddots & \ddots & \\ & & \tilde{p}_{n - 2} & \tilde{q}_{n - 2} & \tilde{r}_{n - 2} \\ & & & & 1 \\ \end{pmatrix} \begin{pmatrix} y_0 \\ y_1 \\ \vdots\\ y_{n-2} \\ y_{n-1} \end{pmatrix} = \begin{pmatrix} \alpha \\ \tilde{s}_1 \\ \vdots \\ \tilde{s}_{n-2} \\ \beta \end{pmatrix}. \end{split}\]

Exercise 11.3: Implement the finite diference method for the \( 2^{nd} \) order linear ordinary differential equation.

def finite_difference_method(a, b, n, p, q, r, s, alpha, beta):
    """
    Solves 2nd order linear ODE in the form p(x)y(x)'' + q(x)y(x)' + r(x)y(x) + s(x) = 0 in [a, b] with y(a) = alpha,
    y(b) = beta using finite difference method.
    Args:
        a (float): The left-most point of the interval
        b (float): The right-most point of the interval
        n (int): Number of points dividing interval [a, b]
        p (function): Arbitrary differentiable function
        q (function): Arbitrary differentiable function
        r (function): Arbitrary differentiable function
        s (function): Arbitrary differentiable function
        alpha (float): Initial value at x = a
        beta (float): Initial value at x = b
    Returns:
        numpy.ndarray: Vector of grid-points
        numpy.ndarray: Vector of solution
    """

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_bvp function.

try:
    np.testing.assert_array_almost_equal(finite_difference_method(0.0, 1.0, 1000, lambda x: 1.0, lambda x: 0.0, \
        lambda x: 0.0, lambda x: -x, 0.0, 0.0)[1], integrate.solve_bvp(lambda x, y: np.vstack((y[1], x)), \
        lambda a, b: np.array([a[0], b[0]]), np.linspace(0.0, 1.0, 1000), np.zeros((2, 1000))).y[0], decimal=5)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 11.4: Solve the following ordinary differential equation,

\[ y^{\prime \prime} \left( t \right) + 2 y^{\prime} \left( t \right) + 100 y \left( t \right) = 0, \quad t \in [0, \, 2], \quad y \left( 0 \right) = 1, \quad y \left( 2 \right) = 0. \]

using the finite difference method with a partition of size \( N = 20 \). Plot the solution obtained by the finite difference method as well as the analytical solution

\[ y \left( x \right) = -\frac{\exp{\left( -t \right)} \sin \left[ 3 \sqrt{11} \left( t - 2 \right) \right]}{\sin \left( 6 \sqrt{11} \right) } \]

in the interval \( t \in [0, \, 2] \).

# add your code here

Note: The ordinary differential equation in Exercise 11.4 is the equation of damped harmonic oscillator.

Finite element method

The finite element method is a numerical technique for solving differential equations, commonly in weak formulation, by applying linear constraints determined by finite sets of basis functions. Consider the following boundary value problem,

\[ \tag{5} \label{eq:fem_bvp} y^{\prime \prime} \left( x \right) = f \left( x \right), \quad x \in \left[a, \, b \right], \quad y \left( a \right) = \alpha, \quad y \left( b \right) = \beta. \]

We are looking for the solution of \( \eqref{eq:fem_bvp} \) in the folowing form,

\[ \tilde{y} \left(x \right) = \sum_{j = 0}^{n - 1} q_j \phi_j \left( x \right), \]

where \( q_j \) are unknown constant coefficients and \( \phi_j \left( x \right) \) are basis functions. Let us define the weighted residual \( \Pi_i \left( \tilde{y} \right) \) for node \( i \) as follows,

\[ \Pi_i \left( \tilde{y} \right) = \int_a^b w_i \left( x \right) \left[ \tilde{y}^{\prime \prime} \left( x \right) - f \left( x \right) \right] dx, \quad i = 0, \dots, n - 1, \]

where \( w_i \left( x \right) \) is a test function. Using the Galerkin method, i.e., \( w_i \left( x \right) \) = \( \phi_i \left( x \right) \), and expanding the above integral one gets

\[\begin{split} \Pi_i \left( \tilde{y} \right) = \int_a^b \phi_i \left( x \right) \left[ \tilde{y}^{\prime \prime} \left( x \right) - f \left( x \right) \right] dx = \int_a^b \phi_i \left( x \right) \tilde{y}^{\prime \prime} \left( x \right) dx - \int_a^b \phi_i \left( x \right) f \left( x \right) dx = \\ = \left[ \phi_i \left( x \right) \tilde{y}^{\prime} \left( x \right) \right]^b_a - \int_a^b \phi_i^{\prime} \left( x \right) \tilde{y}^{\prime} \left( x \right) dx - \int_a^b \phi_i \left( x \right) f \left( x \right) dx. \end{split}\]

By defining a particular basis function \( \phi_i \) as

\[\begin{split} \phi_i(x) = \begin{cases} \frac{x - x_{i - 1}}{h} & \mathrm{if} \ x \in \left[x_{i - 1}, x_i \right) \\ \frac{x_{i + 1} - x}{h} & \mathrm{if} \ x \in \left[x_{i}, x_{i + 1} \right] \\ 0 & \mathrm{otherwise} \end{cases} \end{split}\]

where \( h > 0 \) is a step size, and setting each \( \Pi_i \left( \tilde{y} \right) \) to zero, one obtains a system of linear algebraic equations for unknowns \( q_i \),

\[ \tag{6} \label{eq:fem_system} q_{i - 1} - 2 q_{i} + q_{i + 1} = h^2 f \left( x_i \right), \quad q_0 = \alpha, \quad q_{n - 1} = \beta. \]

In the matrix form, system of equations \( \eqref{eq:fem_system} \) reads as

\[\begin{split} \begin{pmatrix} 1 & & & & \\ 1 & -2 & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & -2 & 1 \\ & & & & 1 \\ \end{pmatrix} \begin{pmatrix} q_0 \\ q_1 \\ \vdots\\ q_{n-2} \\ q_{n-1} \end{pmatrix} = \begin{pmatrix} \alpha \\ h^2 f \left( x_1 \right) \\ \vdots \\ h^2 f \left( x_{n - 2} \right) \\ \beta \end{pmatrix}. \end{split}\]

Exercise 11.5: Implement the finite element method for the \( 2^{nd} \) order ordinary differential equation in the form \( y^{\prime \prime} \left( x \right) = f \left( x \right) \).

def finite_element_method(f, a, b, n, alpha, beta):
    """
    Solves 2nd order ODE in the form y(x)'' = f(x) in [a, b] with y(a) = alpha,
    y(b) = beta using finite element method.
    Args:
        f (function): Funtion
        a (float): The left-most point of the interval [a, b]
        b (float): The right-most point of the interval [a, b]
        n (int): Number of points dividing interval [a, b]
        alpha (float): Initial value at x = a
        beta (float): Initial value at x = b
    Returns:
        numpy.ndarray: Vector of grid-points
        numpy.ndarray: Vector of solution
    """

    # add your code here

You may evaluate the correctness of your implementation using the scipy.integrate.solve_bvp function.

try:
    np.testing.assert_array_almost_equal(finite_element_method(lambda x: x, 0.0, 1.0, 1000, 0.0, 0.0)[1], integrate.solve_bvp(lambda x, y: np.vstack((y[1], x)), \
        lambda a, b: np.array([a[0], b[0]]), np.linspace(0.0, 1.0, 1000), np.zeros((2, 1000))).y[0], decimal=5)
except AssertionError as E:
    print(E)
else:
    print("The implementation is correct.")

Exercise 11.6: Solve the following ordinary differential equation,

\[ T^{\prime \prime} \left( x \right) = -50 e^x, \quad x \in \left[ -1, \, 1 \right], \quad T \left( -1 \right) = 0, \quad T \left( 1 \right) = 10, \]

using the finite element method with a partition of size \( N = 20 \). Plot the solution obtained by the finite element method as well as the analytical solution

\[ T \left( x \right) \approx -50 e^x + 63.76 x + 82.15 \]

in the interval \( x \in \left[ -1, \, 1 \right] \).

# add your code here

Note: The ordinary differential equation in Exercise 11.6 is the Poisson’s equation.

Note: The finite diference and finite element methods lead to system of linear algebraic equations with banded matrix. These systems can be solved, e.g., using the Thomas algorithm.

previous

Tutorial 10 - Initial value problems

next

Solutions to tutorial exercises