{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Nonlinear BVP using finite difference method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Conside the BVP\n",
"$$\n",
"y'' - 2 \\frac{(y')^2}{y} + y = 0, \\qquad x \\in (-1,+1)\n",
"$$\n",
"with boundary conditions\n",
"$$\n",
"y(-1) = y(1) = \\frac{1}{e + e^{-1}}\n",
"$$\n",
"The exact solution is given by\n",
"$$\n",
"y(x) = \\frac{1}{\\exp(x) + \\exp(-x)}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finite difference method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Make a uniform grid of $n+1$ points such that\n",
"$$\n",
"x_j = -1 + j h, \\qquad 0 \\le j \\le n, \\qquad h = 2/n\n",
"$$\n",
"Then we approximate the differential equation by finite differences\n",
"$$\n",
"\\phi_j = \\frac{y_{j-1}-2y_j + y_{j+1}}{h^2} - \\frac{2}{y_j} \\left(\\frac{y_{j+1} - y_{j-1}}{2h}\\right)^2 + y_j = 0, \\quad 1 \\le j \\le n-1\n",
"$$\n",
"with the boundary condition\n",
"$$\n",
"y_0 = y_n = \\frac{1}{e + e^{-1}}\n",
"$$\n",
"This is a non-linear system of equations. We can eliminate $y_0, y_n$ and form a problem to find the remaining unknowns. Define\n",
"$$\n",
"Y = [y_1, \\ldots, y_{n-1}], \\qquad \\phi = [\\phi_1, \\ldots, \\phi_{n-1}]\n",
"$$\n",
"Then $\\phi = \\phi(Y)$ is a map from $R^{n-1}$ to $R^{n-1}$. The Newton method is given by\n",
"$$\n",
"\\phi'(Y^k) \\delta Y^k = - \\phi(Y^k), \\qquad Y^{k+1} = Y^k + \\delta Y^k\n",
"$$\n",
"The Jacobian is given by\n",
"$$\n",
"[\\phi']_{jl} = \\frac{\\partial \\phi_j}{\\partial y_l}, \\quad 1 \\le j,l \\le n-1\n",
"$$\n",
"and it is a tri-diagonal matrix."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def yexact(x):\n",
" return 1.0/(np.exp(x)+np.exp(-x))\n",
"\n",
"# Computes the vector phi(y)\n",
"def phi(y):\n",
" n = len(y) - 1\n",
" h = 2.0/n\n",
" res = np.zeros(n-1)\n",
" for i in range(1,n):\n",
" dy = (y[i+1] - y[i-1])/(2.0*h)\n",
" res[i-1] = (y[i-1] - 2.0*y[i] + y[i+1])/h**2 - 2.0*dy**2/y[i] + y[i]\n",
" return res\n",
"\n",
"# Computes Jacobian d(phi)/dy\n",
"def dphi(y):\n",
" n = len(y) - 1\n",
" h = 2.0/n\n",
" res = np.zeros((n-1,n-1))\n",
" for i in range(1,n):\n",
" dy = (y[i+1] - y[i-1])/(2.0*h)\n",
" if i > 1:\n",
" res[i-1,i-2] = 1.0/h**2 + 4.0*dy/y[i]\n",
" res[i-1,i-1] = -2.0/h**2 + 1.0 - 2.0*(dy/y[i])**2\n",
" if i < n-1:\n",
" res[i-1,i] = 1.0/h**2 - 4.0*dy/y[i]\n",
" return res"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now solve the problem with Newton method."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of iterations = 13\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"n = 50\n",
"\n",
"# Initial guess for y\n",
"y = np.zeros(n+1)\n",
"y[:] = 1.0/(np.exp(1) + np.exp(-1))\n",
"\n",
"maxiter = 100\n",
"TOL = 1.0e-6\n",
"it = 0\n",
"while it < maxiter:\n",
" b = -phi(y)\n",
" if np.linalg.norm(b) < TOL:\n",
" break\n",
" A = dphi(y)\n",
" v = np.linalg.solve(A,b)\n",
" y[1:n] = y[1:n] + v\n",
" it += 1\n",
"print(\"Number of iterations = %d\" % it)\n",
"x = np.linspace(-1.0,1.0,n+1)\n",
"plt.plot(x,y,'o',x,yexact(x))\n",
"plt.legend((\"Numerical\",\"Exact\"));"
]
}
],
"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.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 1
}