{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"# Tutorial 11 - Boundary value problems"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Boundary value problems of ordinary differential equations, finite difference method, shooting method, finite element method."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy import linalg, integrate, optimize"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Shooting method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The [shooting method](https://en.wikipedia.org/wiki/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,\n",
"\n",
"$$\n",
"\\tag{1} \\label{eq:shooting_bvp}\n",
"y^{\\prime \\prime}(x) = f(x, y(x), y^{\\prime}(x)), \\quad x \\in [a, b], \\quad y(a) = \\alpha, \\quad y(b) = \\beta.\n",
"$$\n",
"\n",
"Let $ y(x; \\gamma) $ denote the solution of the initial value problem\n",
"\n",
"$$\n",
"\\tag{2} \\label{eq:shooting_ivp}\n",
"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.\n",
"$$\n",
"\n",
"Define the function $ F(\\gamma) $ as the difference between $ y(b; \\, \\gamma) $ and the specified boundary value $ \\beta $,\n",
"\n",
"$$\n",
"F(\\gamma) = y(b; \\, \\gamma) - \\beta.\n",
"$$\n",
"\n",
"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](https://en.wikipedia.org/wiki/Bisection_method) or [Newton-Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Note:** The $ 2^{nd} $ order explicit ordinary differential equation \n",
"\n",
"$$ \n",
"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\n",
"$$\n",
"\n",
"can be expressed as a system of two $ 1^{st} $ order differential equations\n",
"\n",
"$$\n",
"\\begin{align}\n",
"& y^{\\prime}(x) = z(x), \\quad y(a) = \\alpha_1 \\\\\n",
"& z^{\\prime}(x) = f(x, y(x), z(x)), \\quad z(a) = \\alpha_2\n",
"\\end{align}\n",
"$$\n",
"\n",
"by defining a new variable $ z(x) = y^{\\prime}(x) $. \n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 11.1:** Implement the shooting method for the $ 2^{nd} $ order explicit ordinary differential equation.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def shooting_method(f, g, a, b, n, alpha, beta):\n",
" \"\"\"\n",
" 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\n",
" using shooting method.\n",
" Args:\n",
" f (function): function \n",
" g (function): function\n",
" a (float): The left-most point of the interval\n",
" b (float): The right-most point of the interval\n",
" n (int): number of points dividing interval [a, b]\n",
" alpha (float): initial value at x = a\n",
" beta (float): initial value at x = b\n",
" Returns:\n",
" numpy.ndarray: Vector of grid-points\n",
" numpy.ndarray: Vector of solution\n",
" \"\"\"\n",
"\n",
" # add your code here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_bvp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"try:\n",
" np.testing.assert_array_almost_equal(shooting_method(lambda x, y, z: x, lambda x, y, z: z, \\\n",
" 0.0, 1.0, 1000, 0.0, 0.0)[1], integrate.solve_bvp(lambda x, y: np.vstack((y[1], x)), \\\n",
" lambda a, b: np.array([a[0], b[0]]), np.linspace(0.0, 1.0, 1000), np.zeros((2, 1000))).y[0], decimal=5)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 11.2:** Solve the following ordinary differential equation,\n",
"\n",
"$$\n",
"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\n",
"$$\n",
"\n",
"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\n",
"\n",
"$$\n",
"y \\left( x \\right) \\approx 2 \\cosh \\left( \\frac{x}{2} \\right) - 2.255\n",
"$$\n",
" \n",
"in the interval $ x \\in [-1, \\, 1] $.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Note:** The solution of the equation in Exercise 11.2 is a [catenary](https://en.wikipedia.org/wiki/Catenary) curve.\n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finite difference method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The [finite difference method](https://en.wikipedia.org/wiki/Finite_difference_method) is a numerical technique for solving differential equations by approximating derivatives with finite differences. Consider the following boundary value problem,\n",
"\n",
"$$\n",
"\\tag{3} \\label{eq:fdm_bvp}\n",
"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.\n",
"$$\n",
"\n",
"Let \n",
"\n",
"$$\n",
"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,\n",
"$$\n",
"\n",
"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\n",
"\n",
"$$\n",
"\\tag{4} \\label{eq:fdm_system}\n",
"\\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.\n",
"$$\n",
"\n",
"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\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
"1 & & & & \\\\\n",
"\\tilde{p}_1 & \\tilde{q}_1 & \\tilde{r}_1 & & \\\\\n",
" & \\ddots & \\ddots & \\ddots & \\\\\n",
" & & \\tilde{p}_{n - 2} & \\tilde{q}_{n - 2} & \\tilde{r}_{n - 2} \\\\\n",
" & & & & 1 \\\\\n",
"\\end{pmatrix} \n",
"\\begin{pmatrix}\n",
"y_0 \\\\\n",
"y_1 \\\\\n",
"\\vdots\\\\\n",
"y_{n-2} \\\\\n",
"y_{n-1}\n",
"\\end{pmatrix} =\n",
"\\begin{pmatrix}\n",
"\\alpha \\\\\n",
"\\tilde{s}_1 \\\\\n",
"\\vdots \\\\\n",
"\\tilde{s}_{n-2} \\\\\n",
"\\beta\n",
"\\end{pmatrix}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 11.3:** Implement the finite diference method for the $ 2^{nd} $ order linear ordinary differential equation.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def finite_difference_method(a, b, n, p, q, r, s, alpha, beta):\n",
" \"\"\"\n",
" 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,\n",
" y(b) = beta using finite difference method.\n",
" Args:\n",
" a (float): The left-most point of the interval\n",
" b (float): The right-most point of the interval\n",
" n (int): Number of points dividing interval [a, b]\n",
" p (function): Arbitrary differentiable function\n",
" q (function): Arbitrary differentiable function\n",
" r (function): Arbitrary differentiable function\n",
" s (function): Arbitrary differentiable function\n",
" alpha (float): Initial value at x = a\n",
" beta (float): Initial value at x = b\n",
" Returns:\n",
" numpy.ndarray: Vector of grid-points\n",
" numpy.ndarray: Vector of solution\n",
" \"\"\"\n",
"\n",
" # add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_bvp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"try:\n",
" np.testing.assert_array_almost_equal(finite_difference_method(0.0, 1.0, 1000, lambda x: 1.0, lambda x: 0.0, \\\n",
" lambda x: 0.0, lambda x: -x, 0.0, 0.0)[1], integrate.solve_bvp(lambda x, y: np.vstack((y[1], x)), \\\n",
" lambda a, b: np.array([a[0], b[0]]), np.linspace(0.0, 1.0, 1000), np.zeros((2, 1000))).y[0], decimal=5)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 11.4:** Solve the following ordinary differential equation,\n",
" \n",
"$$\n",
"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.\n",
"$$\n",
" \n",
"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\n",
" \n",
"$$\n",
"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) }\n",
"$$\n",
"\n",
"in the interval $ t \\in [0, \\, 2] $. \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"# add your code here\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Note:** The ordinary differential equation in Exercise 11.4 is the equation of [damped harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator#Damped_harmonic_oscillator). \n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finite element method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The [finite element method](https://en.wikipedia.org/wiki/Finite_element_method) is a numerical technique for solving differential equations, commonly in [weak formulation](https://en.wikipedia.org/wiki/Weak_formulation), by applying linear constraints determined by finite sets of basis functions. Consider the following boundary value problem,\n",
"\n",
"$$\n",
"\\tag{5} \\label{eq:fem_bvp}\n",
"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.\n",
"$$\n",
"\n",
"We are looking for the solution of $ \\eqref{eq:fem_bvp} $ in the folowing form,\n",
"\n",
"$$\n",
"\\tilde{y} \\left(x \\right) = \\sum_{j = 0}^{n - 1} q_j \\phi_j \\left( x \\right),\n",
"$$\n",
"\n",
"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,\n",
"\n",
"$$\n",
"\\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,\n",
"$$\n",
"\n",
"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\n",
"\n",
"$$\n",
"\\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.\n",
"$$\n",
"\n",
"By defining a particular basis function $ \\phi_i $ as\n",
"\n",
"$$ \n",
"\\phi_i(x) =\n",
"\\begin{cases} \n",
"\\frac{x - x_{i - 1}}{h} & \\mathrm{if} \\ x \\in \\left[x_{i - 1}, x_i \\right) \\\\\n",
"\\frac{x_{i + 1} - x}{h} & \\mathrm{if} \\ x \\in \\left[x_{i}, x_{i + 1} \\right] \\\\\n",
"0 & \\mathrm{otherwise}\n",
"\\end{cases}\n",
"$$\n",
"\n",
"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 $,\n",
"\n",
"$$\n",
"\\tag{6} \\label{eq:fem_system}\n",
"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.\n",
"$$\n",
"\n",
"In the matrix form, system of equations $ \\eqref{eq:fem_system} $ reads as\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
"1 & & & & \\\\\n",
"1 & -2 & 1 & & \\\\\n",
" & \\ddots & \\ddots & \\ddots & \\\\\n",
" & & 1 & -2 & 1 \\\\\n",
" & & & & 1 \\\\\n",
"\\end{pmatrix} \n",
"\\begin{pmatrix}\n",
"q_0 \\\\\n",
"q_1 \\\\\n",
"\\vdots\\\\\n",
"q_{n-2} \\\\\n",
"q_{n-1}\n",
"\\end{pmatrix} =\n",
"\\begin{pmatrix}\n",
"\\alpha \\\\\n",
"h^2 f \\left( x_1 \\right) \\\\\n",
"\\vdots \\\\\n",
"h^2 f \\left( x_{n - 2} \\right) \\\\\n",
"\\beta\n",
"\\end{pmatrix}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**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) $.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def finite_element_method(f, a, b, n, alpha, beta):\n",
" \"\"\"\n",
" Solves 2nd order ODE in the form y(x)'' = f(x) in [a, b] with y(a) = alpha,\n",
" y(b) = beta using finite element method.\n",
" Args:\n",
" f (function): Funtion\n",
" a (float): The left-most point of the interval [a, b]\n",
" b (float): The right-most point of the interval [a, b]\n",
" n (int): Number of points dividing interval [a, b]\n",
" alpha (float): Initial value at x = a\n",
" beta (float): Initial value at x = b\n",
" Returns:\n",
" numpy.ndarray: Vector of grid-points\n",
" numpy.ndarray: Vector of solution\n",
" \"\"\"\n",
"\n",
" # add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.integrate.solve_bvp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"try:\n",
" 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)), \\\n",
" lambda a, b: np.array([a[0], b[0]]), np.linspace(0.0, 1.0, 1000), np.zeros((2, 1000))).y[0], decimal=5)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 11.6:** Solve the following ordinary differential equation,\n",
" \n",
"$$\n",
"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,\n",
"$$\n",
" \n",
"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\n",
" \n",
"$$\n",
"T \\left( x \\right) \\approx -50 e^x + 63.76 x + 82.15\n",
"$$\n",
"\n",
"in the interval $ x \\in \\left[ -1, \\, 1 \\right] $. \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Note:** The ordinary differential equation in Exercise 11.6 is the [Poisson's equation](https://en.wikipedia.org/wiki/Poisson%27s_equation).\n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Note:** The finite diference and finite element methods lead to system of linear algebraic equations with [banded matrix](https://en.wikipedia.org/wiki/Band_matrix). These systems can be solved, e.g., using the [Thomas algorithm](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm).\n",
" \n",
"
"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}