{
"cells": [
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2021-08-27 Finite Difference 2\n",
"\n",
"## Last time\n",
"* Evaluating derivatives\n",
"* Taylor series and truncation error\n",
"* Activity on stability of computing derivatives\n",
"\n",
"## Today\n",
"\n",
"* Stability\n",
"* Derivatives as matrices\n",
"* Boundary conditions\n",
"* Discrete eigenvalues/eigenvectors"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "009bd630",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"using Plots"
]
},
{
"cell_type": "markdown",
"id": "3d1ef4b2",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Discuss activity on stability and accuracy computing derivatives"
]
},
{
"cell_type": "markdown",
"id": "cf428730",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# A note on stable computation"
]
},
{
"cell_type": "code",
"execution_count": 93,
"id": "136c9b3e",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x = 1.0e-15\n",
"log(1 + x) = 1.110223024625156e-15\n"
]
},
{
"data": {
"text/plain": [
"1.110223024625156e-15"
]
},
"execution_count": 93,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = 1e-15\n",
"@show x\n",
"@show log(1 + x)"
]
},
{
"cell_type": "code",
"execution_count": 95,
"id": "d5a1ad88",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 95,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot([x -> log(1 + x), log1p], xlims=(-1e-15, 1e-15))\n",
"plot!(x -> x)"
]
},
{
"cell_type": "markdown",
"id": "7625613a",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# A few methods on grids"
]
},
{
"cell_type": "code",
"execution_count": 99,
"id": "716028d5",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [],
"source": [
"diff1l(x, u) = x[2:end], (u[2:end] - u[1:end-1]) ./ (x[2:end] - x[1:end-1])\n",
"diff1r(x, u) = x[1:end-1], (u[2:end] - u[1:end-1]) ./ (x[2:end] - x[1:end-1])\n",
"diff1c(x, u) = x[2:end-1], (u[3:end] - u[1:end-2]) ./ (x[3:end] - x[1:end-2])\n",
"difflist = [diff1l, diff1r, diff1c]\n",
"\n",
"n = 20\n",
"h = 2 / (n - 1)\n",
"x = LinRange(-3, 3, n)\n",
"u = sin.(x)\n",
"fig = plot(cos, xlims=(-3, 3))\n",
"for d in difflist\n",
" xx, yy = d(x, u)\n",
" plot!(fig, xx, yy, marker=:circle, label=d)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 100,
"id": "e6660042",
"metadata": {
"cell_style": "split",
"scrolled": false,
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 100,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig"
]
},
{
"cell_type": "markdown",
"id": "680ffb31",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Measuring error on grids"
]
},
{
"cell_type": "code",
"execution_count": 103,
"id": "8da2b252",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"refinement_error (generic function with 1 method)"
]
},
"execution_count": 103,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using LinearAlgebra\n",
"\n",
"grids = 2 .^ (2:10)\n",
"hs = 1 ./ grids\n",
"function refinement_error(f, fprime, d)\n",
" error = []\n",
" for n in grids\n",
" x = LinRange(-3, 3, n)\n",
" xx, yy = d(x, f.(x))\n",
" push!(error, norm(yy - fprime.(xx), Inf))\n",
" end\n",
" error\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 105,
"id": "2df531e7",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 105,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig = plot(xscale=:log10, yscale=:log10)\n",
"for d in difflist\n",
" error = refinement_error(sin, cos, d)\n",
" plot!(fig, hs, error, marker=:circle, label=d)\n",
"end\n",
"plot!(fig, hs, hs .^ 2)"
]
},
{
"cell_type": "markdown",
"id": "ff3425ff",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Stability\n",
"\n",
"Are there \"rough\" functions for which these formulas estimate $u'(x_i) = 0$?"
]
},
{
"cell_type": "code",
"execution_count": 107,
"id": "64bad794",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 107,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = LinRange(-1, 1, 9)\n",
"f_rough(x) = cos(.1 + 4Ï€*x)\n",
"fp_rough(x) = -4Ï€*sin(.1 + 4Ï€*x)\n",
"\n",
"plot(x, f_rough, marker=:circle)\n",
"plot!(f_rough)"
]
},
{
"cell_type": "code",
"execution_count": 108,
"id": "8509f935",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 108,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig = plot(fp_rough, xlims=(-1, 1))\n",
"for d in difflist\n",
" xx, yy = d(x, f_rough.(x))\n",
" plot!(fig, xx, yy, label=d, marker=:circle)\n",
"end\n",
"fig"
]
},
{
"cell_type": "markdown",
"id": "f298a8e0",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"If we have a solution $u(x)$, then $u(x) + f_{\\text{rough}}(x)$ is indistinguishable to our FD method."
]
},
{
"cell_type": "markdown",
"id": "93b9924e",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Second derivatives\n",
"\n",
"We can compute a second derivative by applying first derivatives twice."
]
},
{
"cell_type": "code",
"execution_count": 113,
"id": "bc169f00",
"metadata": {
"cell_style": "split"
},
"outputs": [],
"source": [
"function diff2a(x, u)\n",
" xx, yy = diff1c(x, u)\n",
" diff1c(xx, yy)\n",
"end\n",
"\n",
"function diff2b(x, u)\n",
" xx, yy = diff1l(x, u)\n",
" diff1r(xx, yy)\n",
"end\n",
"\n",
"diff2list = [diff2a, diff2b]\n",
"n = 10\n",
"x = LinRange(-3, 3, n)\n",
"u = - cos.(x);"
]
},
{
"cell_type": "code",
"execution_count": 114,
"id": "17f737fd",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 114,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig = plot(cos, xlims=(-3, 3))\n",
"for d2 in diff2list\n",
" xx, yy = d2(x, u)\n",
" plot!(fig, xx, yy, marker=:circle, label=d2)\n",
"end\n",
"fig"
]
},
{
"cell_type": "markdown",
"id": "c889f810",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# How fast do these approximations converge?"
]
},
{
"cell_type": "code",
"execution_count": 111,
"id": "11b40669",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"refinement_error2 (generic function with 1 method)"
]
},
"execution_count": 111,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"grids = 2 .^ (3:10)\n",
"hs = 1 ./ grids\n",
"function refinement_error2(f, f_xx, d2)\n",
" error = []\n",
" for n in grids\n",
" x = LinRange(-3, 3, n)\n",
" xx, yy = d2(x, f.(x))\n",
" push!(error, norm(yy - f_xx.(xx), Inf))\n",
" end\n",
" error\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 112,
"id": "bbf49a0c",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 112,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig = plot(xscale=:log10, yscale=:log10)\n",
"for d2 in diff2list\n",
" error = refinement_error2(x -> -cos(x), cos, d2)\n",
" plot!(fig, hs, error, marker=:circle, label=d2)\n",
"end\n",
"plot!(fig, hs, hs .^ 2) "
]
},
{
"cell_type": "markdown",
"id": "52e5295b",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"* Both methods are second order accurate.\n",
"* The `diff2b` method is more accurate than `diff2a` (by a factor of 4)\n",
"* The `diff2b` method can't compute derivatives at points adjacent the boundary.\n",
"* We don't know yet whether either is stable"
]
},
{
"cell_type": "markdown",
"id": "60810431",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Differentiation matrices\n",
"\n",
"All our `diff*` functions thus far have been linear in `u`, therefore they can be represented as matrices.\n",
"$$\\frac{u_{i+1} - u_i}{x_{i+1} - x_i} = \\begin{bmatrix} -1/h & 1/h \\end{bmatrix} \\begin{bmatrix} u_i \\\\ u_{i+1} \\end{bmatrix}$$"
]
},
{
"cell_type": "code",
"execution_count": 115,
"id": "f916afef",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"diff1_mat (generic function with 1 method)"
]
},
"execution_count": 115,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function diff1_mat(x)\n",
" n = length(x)\n",
" D = zeros(n, n)\n",
" h = x[2] - x[1]\n",
" D[1, 1:2] = [-1/h 1/h]\n",
" for i in 2:n-1\n",
" D[i, i-1:i+1] = [-1/2h 0 1/2h]\n",
" end\n",
" D[n, n-1:n] = [-1/h 1/h]\n",
" D\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 119,
"id": "0f7cca47",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 119,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = LinRange(-3, 3, 10)\n",
"plot(x, diff1_mat(x) * sin.(x), marker=:circle)\n",
"plot!(cos)"
]
},
{
"cell_type": "markdown",
"id": "1898e284",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# How accurate is this derivative matrix?"
]
},
{
"cell_type": "code",
"execution_count": 120,
"id": "af598cf6",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 120,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig = plot(xscale=:log10, yscale=:log10, legend=:topleft)\n",
"error = refinement_error(sin, cos, (x, u) -> (x, diff1_mat(x) * u))\n",
"plot!(fig, hs, error, marker=:circle)\n",
"plot!(fig, hs, hs, label=\"\\$h\\$\")\n",
"plot!(fig, hs, hs .^ 2, label=\"\\$h^2\\$\")"
]
},
{
"cell_type": "markdown",
"id": "96a6d2d9",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Can we study it as a matrix?"
]
},
{
"cell_type": "code",
"execution_count": 121,
"id": "b7b7a95f",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 121,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"D = diff1_mat(x)\n",
"spy(D, marker=(:square, 10), c=:bwr)"
]
},
{
"cell_type": "code",
"execution_count": 133,
"id": "79534945",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"10-element Vector{Float64}:\n",
" 2.268133218393964\n",
" 2.2674392839412794\n",
" 1.4265847744427302\n",
" 1.3683737968309653\n",
" 1.2135254915624207\n",
" 1.0228485194005281\n",
" 0.8816778784387094\n",
" 0.5437139466339259\n",
" 0.46352549156242107\n",
" 3.204643550177702e-17"
]
},
"execution_count": 133,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"svdvals(D)"
]
},
{
"cell_type": "markdown",
"id": "6bea73e3",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Second derivative with Dirichlet boundary conditions"
]
},
{
"cell_type": "markdown",
"id": "e737d769",
"metadata": {
"cell_style": "split"
},
"source": [
"The left endpoint in our example boundary value problem has a Dirichlet boundary condition,\n",
"$$u(-1) = a . $$\n",
"With finite difference methods, we have an explicit degree of freedom $u_0 = u(x_0 = -1)$ at that endpoint.\n",
"When building a matrix system for the BVP, we can implement this boundary condition by modifying the first row of the matrix,\n",
"$$ \\begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\\\ \\\\ & & A_{2:,:} & & \\\\ \\\\ \\end{bmatrix} \\begin{bmatrix} u_0 \\\\ \\\\ u_{2:} \\\\ \\\\ \\end{bmatrix} = \\begin{bmatrix} a \\\\ \\\\ f_{2:} \\\\ \\\\ \\end{bmatrix} . $$\n",
"\n",
"* This matrix is not symmetric even if $A$ is."
]
},
{
"cell_type": "code",
"execution_count": 73,
"id": "ac080794",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"laplacian_dirichlet (generic function with 1 method)"
]
},
"execution_count": 73,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function laplacian_dirichlet(x)\n",
" n = length(x)\n",
" D = zeros(n, n)\n",
" h = x[2] - x[1]\n",
" D[1, 1] = 1\n",
" for i in 2:n-1\n",
" D[i, i-1:i+1] = (1/h^2) * [-1, 2, -1]\n",
" end\n",
" D[n, n] = 1\n",
" D\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "b4d0be45",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Laplacian as a matrix"
]
},
{
"cell_type": "code",
"execution_count": 123,
"id": "db33a2cf",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 123,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L = laplacian_dirichlet(x)\n",
"spy(L, marker=(:square, 10), c=:bwr)"
]
},
{
"cell_type": "code",
"execution_count": 124,
"id": "e95b53e2",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"text/plain": [
"36.554572047954174"
]
},
"execution_count": 124,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cond(L)"
]
},
{
"cell_type": "markdown",
"id": "718abdde",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Discrete \"Green's functions\""
]
},
{
"cell_type": "code",
"execution_count": 81,
"id": "33decdc7",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(x, inv(L)[:, 2])"
]
},
{
"cell_type": "code",
"execution_count": 82,
"id": "aaa0112d",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Ln = copy(L)\n",
"Ln[n, n-1] = -1\n",
"plot(x, inv(Ln)[:, 4])"
]
},
{
"cell_type": "markdown",
"id": "b0b9e6b6",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Discrete eigenfunctions"
]
},
{
"cell_type": "code",
"execution_count": 125,
"id": "28404d95",
"metadata": {
"cell_style": "split",
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 125,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Lambda, V = eigen(L)\n",
"plot(Lambda, marker=:circle)"
]
},
{
"cell_type": "code",
"execution_count": 131,
"id": "450e6556",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 131,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(x, V[:, end])"
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Julia 1.6.2",
"language": "julia",
"name": "julia-1.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.6.2"
},
"rise": {
"enable_chalkboard": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}