{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# p01: Convergence of fourth order finite differences"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us compute the derivative of\n",
"\n",
"$$\n",
"u(x) = \\exp(\\sin(x)), \\qquad x \\in [-\\pi,\\pi]\n",
"$$\n",
"\n",
"using fourth order finite difference scheme\n",
"\n",
"$$\n",
"u'(x_j) \\approx w_j = \\frac{1}{h} \\left( \\frac{1}{12} u_{j-2} - \\frac{2}{3} u_{j-1} + \\frac{2}{3} u_{j+1} - \\frac{1}{12} u_{j+2} \\right)\n",
"$$\n",
"\n",
"using periodic boundary conditions."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%config InlineBackend.figure_format='svg'\n",
"from scipy.sparse import coo_matrix\n",
"from numpy import arange,pi,exp,sin,cos,ones,inf\n",
"from numpy.linalg import norm\n",
"from matplotlib.pyplot import figure,loglog,semilogy,text,grid,xlabel,ylabel,title"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The differentiation matrix is anti-symmetric; we first construct the upper triangular part $U$ and get\n",
"\n",
"$$\n",
"D = \\frac{1}{h} ( U - U^\\top)\n",
"$$\n",
"\n",
"The matrix is stored in [coo_format](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html) which is a sparse matrix format.\n",
"\n",
"The error in the derivative approximation is measured as\n",
"\n",
"$$\n",
"E_N = \\max_{1 \\le j \\le N} |w_j - u_j'|\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Nvec = 2**arange(3,13)\n",
"for N in Nvec:\n",
" h = 2*pi/N\n",
" x = -pi + arange(1,N+1)*h\n",
" u = exp(sin(x)) # function values\n",
" uprime = cos(x)*u # Exact derivative\n",
" e = ones(N)\n",
" e1 = arange(0,N)\n",
" e2 = arange(1,N+1); e2[N-1]=0\n",
" e3 = arange(2,N+2); e3[N-2]=0; e3[N-1]=1;\n",
" D = coo_matrix((2*e/3,(e1,e2)),shape=(N,N)) \\\n",
" - coo_matrix((e/12,(e1,e3)),shape=(N,N))\n",
" D = (D - D.T)/h\n",
" error = norm(D@u-uprime,inf)\n",
" loglog(N,error,'or')\n",
" \n",
"loglog(Nvec,100*Nvec**(-4.0),'--')\n",
"text(105,5e-8,'$N^{-4}$')\n",
"grid(True); xlabel('N'); ylabel('error')\n",
"title('Convergence of fourth-order finite difference');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With about 4000 points, the error is of order $10^{-12}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise\n",
"\n",
"Find sixth and eigth order finite differences, implement them and make error plot as above.\n",
"\n",
"## Exercise\n",
"\n",
"Construct the differentiation matrix using [scipy.toeplitz](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.toeplitz.html) function. You need to pass the first column and first row."
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.12.8"
}
},
"nbformat": 4,
"nbformat_minor": 4
}