{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Numerical Differentiation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To approximate the derivative of a function at a point, we apply finite difference formulas. We examine the application of interpolation and extrapolation for this problem."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Finite Difference Formulas"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The straighforward application of difference formulas leads to an accumulation of roundoff errors and very inaccurate results, when the step size becomes too small."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In our experiment we apply the forward difference formula on the exponential function."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"fwdexp (generic function with 1 method)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fwdexp(x, h) = (exp(x+h) - exp(x))/h"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Suppose we want to approximate the derivative of $\\exp$ at one. As the derivative of $\\exp$ is $\\exp$ again, we approximate $\\exp(1)$, or the constant $e$. In the application of the forward difference formula, we start with $h = 0.1$ and we divide $h$ each time by 10, sixteen times, till we end up with $h = 10^{-16}$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Printf\n",
"xval = 1.0 # the value for x"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Approximating 2.7182818284590451e+00 :\n",
" 1.000e-01 2.8588419548738830e+00 1.406e-01\n",
" 1.000e-02 2.7319186557871245e+00 1.364e-02\n",
" 1.000e-03 2.7196414225332255e+00 1.360e-03\n",
" 1.000e-04 2.7184177470829241e+00 1.359e-04\n",
" 1.000e-05 2.7182954199567173e+00 1.359e-05\n",
" 1.000e-06 2.7182831874306141e+00 1.359e-06\n",
" 1.000e-07 2.7182819684057331e+00 1.399e-07\n",
" 1.000e-08 2.7182818218562939e+00 6.603e-09\n",
" 1.000e-09 2.7182820439008983e+00 2.154e-07\n",
" 1.000e-10 2.7182833761685279e+00 1.548e-06\n",
" 1.000e-11 2.7183144624132174e+00 3.263e-05\n",
" 1.000e-12 2.7187141427020829e+00 4.323e-04\n",
" 1.000e-13 2.7178259642823828e+00 4.559e-04\n",
" 1.000e-14 2.7089441800853815e+00 9.338e-03\n",
" 1.000e-15 3.1086244689504379e+00 3.903e-01\n",
" 1.000e-16 0.0000000000000000e+00 2.718e+00\n"
]
}
],
"source": [
"exact = exp(xval)\n",
"strexact = @sprintf(\"%.16e\", exact)\n",
"println(\"Approximating $strexact :\")\n",
"hval = 0.1 # the start value for h\n",
"for i=1:16\n",
" approx = fwdexp(xval, hval)\n",
" err = abs(exact - approx)\n",
" strhval = @sprintf(\"%.3e\", hval)\n",
" strapprox = @sprintf(\"%.16e\", approx)\n",
" strexact = @sprintf(\"%.16e\", exact)\n",
" strerr = @sprintf(\"%.3e\", err)\n",
" println(\" $strhval $strapprox $strerr\")\n",
" hval = hval/10.0\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The three columns above are respectively the value for $h$, the approximation, and then in the last column the error of the approximation. That the last approximation is zero is because $h$ is less than the machine precision and therefore $\\exp(1.0 + h) = \\exp(1.0)$. The second important observation is that the smallest error occurs at $h = 10^{-8}$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider a log plot of the errors:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"(process:86700): GLib-GIO-WARNING **: 09:16:38.261: Unexpectedly, UWP app `Clipchamp.Clipchamp_2.5.1.0_neutral__yxz26nhyzhsrt' (AUMId `Clipchamp.Clipchamp_yxz26nhyzhsrt!App') supports 41 extensions but has no verbs\n",
"\n",
"(process:86700): GLib-GIO-WARNING **: 09:16:38.277: Unexpectedly, UWP app `Evernote.Evernote_10.46.9.0_x64__q4d96b2w5wcc2' (AUMId `Evernote.Evernote_q4d96b2w5wcc2!Evernote') supports 1 extensions but has no verbs\n"
]
}
],
"source": [
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"aloghvals = [abs(log(10,10.0^(-p))) for p=1:16]\n",
"alogerror = [log(10,abs(fwdexp(1.0, 10.0^(-p)) - exp(1.0))) for p=1:16]\n",
"scatter(aloghvals, alogerror, marker=(:red), label=\"log(10, errors)\")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"savefig(\"forwarderrors.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The plot confirms the computed errors in the loop above."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Differentiation by Interpolation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Taking derivatives of polynomials is straightforward. What if we would interpolate the function and then use the derivative of the interpolating polynomial as the approximation for the derivative?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In our first experiment, we take a cubic polynomial and approximate with a quadratic, using Chebyshev interpolation at three points."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"df (generic function with 1 method)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(x) = x^3 + 2*x^2 - 3*x + 4\n",
"df(x) = 3*x^2 + 4*x - 3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We interpolate ``f`` at the three Chebyshev points over ``[0,2]``."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1.0 1.8660254037844388 3.482050807568878; 1.0 1.0 1.0; 1.0 0.1339745962155613 0.01794919243112269]\n",
"[11.863620667976086, 4.0, 3.6363793320239166]\n",
"[4.25, -5.249999999999999, 4.999999999999999]\n",
"[11.863620667976086, 4.0, 3.636379332023917]\n",
"x1 : 1.0\n",
"the exact derivative : 4.0\n",
"the approximation : 4.749999999999999\n"
]
}
],
"source": [
"x0 = 1 + cos((2*1-1)*pi/6)\n",
"x1 = 1 + cos((2*2-1)*pi/6)\n",
"x2 = 1 + cos((2*3-1)*pi/6)\n",
"X = [1 x0 x0*x0; 1 x1 x1*x1; 1 x2 x2*x2]\n",
"println(X)\n",
"F = [f(x0); f(x1); f(x2)]\n",
"println(F)\n",
"p = X\\F\n",
"println(p)\n",
"pol(x) = p[1] + p[2]*x + p[3]*x^2\n",
"dpol(x) = p[2] + 2*p[3]*x\n",
"P = [pol(x0); pol(x1); pol(x2)]\n",
"println(P)\n",
"df1 = df(x1)\n",
"dpol1 = dpol(x1)\n",
"println(\"x1 : \", x1)\n",
"println(\"the exact derivative : \", df1)\n",
"println(\"the approximation : \", dpol1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This does not look so accurate, a visualization plots the tangent lines."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tf(x) = df1*(x - x1) + f(x1)\n",
"tp(x) = dpol1*(x - x1) + f(x1)\n",
"xr = -0.1:0.1:2.1\n",
"xv = [xr[i] for i=1:length(xr)]\n",
"fv = [f(xr[i]) for i=1:length(xr)]\n",
"pv = [pol(xr[i]) for i=1:length(xr)]\n",
"tfv = [tf(xr[i]) for i=1:length(xr)]\n",
"tpv = [tp(xr[i]) for i=1:length(xr)]\n",
"plot(xv, fv, label=\"f(x) = x^3 + 2*x^2 - 3*x + 4\")\n",
"plot!(xv, pv, label=\"p(x) interpolates f(x)\")\n",
"plot!(xv, tfv, label=\"tangent line to y = f(x)\")\n",
"plot!(xv, tpv, label=\"tangent line to y = p(x)\")\n",
"scatter!([x0, x1, x2], [f(x0), f(x1), f(x2)], marker=(:red),\n",
" label=\"Chebyshev points\")"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"savefig(\"interpoldiff1.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Three points is surely too little, let us see if we take more points, interpolating a quintic with a quartic."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"df (generic function with 1 method)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(x) = x^5 + 2*x^2 - 3*x + 4\n",
"df(x) = 5*x^4 + 4*x - 3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We interpolate ``f`` at the five Chebyshev points over ``[0,2]``."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1.0 1.9510565162951536 3.806621529777781 7.426933740742366 14.490367470967735; 1.0 1.5877852522924731 2.5210620073974725 4.002905075460565 6.355753645142974; 1.0 1.0 1.0 1.0 1.0; 1.0 0.412214747707527 0.1699209982275801 0.07004394141459308 0.028873145638657286; 1.0 0.04894348370484647 0.0023954645971665718 0.00011724238247695871 5.738250636278406e-6]\n",
"[34.03159938841303, 14.370340162879668, 4.0, 3.115099689777541, 3.8579607589297704]\n",
"[4.062500000000011, -4.562500000000013, 8.250000000000009, -8.749999999999996, 4.999999999999997]\n",
"[34.03159938841303, 14.370340162879675, 4.000000000000007, 3.115099689777548, 3.85796075892978]\n",
"x2 : 1.0\n",
"The exact derivative : 6.0\n",
"The approximation : 5.6875000000000036\n"
]
}
],
"source": [
"x0 = 1 + cos((2*1-1)*pi/10)\n",
"x1 = 1 + cos((2*2-1)*pi/10)\n",
"x2 = 1 + cos((2*3-1)*pi/10)\n",
"x3 = 1 + cos((2*4-1)*pi/10)\n",
"x4 = 1 + cos((2*5-1)*pi/10)\n",
"X = [1 x0 x0^2 x0^3 x0^4;\n",
" 1 x1 x1^2 x1^3 x1^4;\n",
" 1 x2 x2^2 x2^3 x2^4;\n",
" 1 x3 x3^2 x3^3 x3^4;\n",
" 1 x4 x4^2 x4^3 x4^4]\n",
"println(X)\n",
"F = [f(x0); f(x1); f(x2); f(x3); f(x4)]\n",
"println(F)\n",
"p = X\\F\n",
"println(p)\n",
"pol(x) = p[1] + p[2]*x + p[3]*x^2 + p[4]*x^3 + p[5]*x^4\n",
"dpol(x) = p[2] + 2*p[3]*x + 3*p[4]*x^2 + 4*p[5]*x^3\n",
"P = [pol(x0); pol(x1); pol(x2); pol(x3); pol(x4)]\n",
"println(P)\n",
"df2 = df(x2)\n",
"dpol2 = dpol(x2)\n",
"println(\"x2 : \", x2)\n",
"println(\"The exact derivative : \", df2)\n",
"println(\"The approximation : \", dpol2)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tf(x) = df2*(x - x2) + f(x2)\n",
"tp(x) = dpol2*(x - x2) + f(x2)\n",
"xr = -0.1:0.1:2.1\n",
"xv = [xr[i] for i=1:length(xr)]\n",
"fv = [f(xr[i]) for i=1:length(xr)]\n",
"pv = [pol(xr[i]) for i=1:length(xr)]\n",
"tfv = [tf(xr[i]) for i=1:length(xr)]\n",
"tpv = [tp(xr[i]) for i=1:length(xr)]\n",
"plot(xv, fv, label=\"f(x) = x^5 + 2*x^2 - 3*x + 4\")\n",
"plot!(xv, pv, label=\"p(x) interpolates f(x)\")\n",
"plot!(xv, tfv, label=\"tangent line to y = f(x)\")\n",
"plot!(xv, tpv, label=\"tangent line to y = p(x)\")\n",
"scatter!([x0, x1, x2, x3, x4], [f(x0), f(x1), f(x2), f(x3), f(x4)],\n",
" marker=(:red), label=\"Chebyshev points\")"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"savefig(\"interpoldiff2.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Taking more Chebyshev interpolation points helps, but our problem is a local one, we would do better to work with more values closer to the point of interest. Instead of interpolation, we must use extrapolation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Richardson Extrapolation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The main idea is to consider the Taylor expansions of forward difference formulas at two different values of $h$, typically at $h$ and $h/2$. By a linear combination of the approximated values we can eliminate one term in the Taylor expansion and obtain a higher order approximation. This idea is coded in the function below."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"richardson"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" richardson(f::Function,z::Float64,h::Float64,n::Int)\n",
"\n",
"returns the triangular table of numerical approximations of\n",
"the derivative of f at z, for the initial value h and n steps.\n",
"\"\"\"\n",
"function richardson(f::Function,z::Float64,h::Float64,n::Int)\n",
" R = zeros(n, n)\n",
" w = h\n",
" for i=1:n\n",
" R[i, 1] = (f(z+w) - f(z))/w\n",
" w = w/2\n",
" end\n",
" for j=2:n\n",
" for i=j:n\n",
" w = 2^(j-1) - 1\n",
" R[i, j] = (2^(j-1)*R[i, j-1] - R[i-1, j-1])/w\n",
" end\n",
" end\n",
" return R\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the experiment, we take again $\\exp$ at 1.0 and use $h = 0.1$ as the initial value."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"using Printf\n",
"Base.show(io::IO, f::Float64) = @printf(io, \"%.8e\", f)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Richardson extrapolation on exp(x) at x = 1.0 :\n",
"4×4 Matrix{Float64}:\n",
" 2.85884195e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
" 2.78738579e+00 2.71592963e+00 0.00000000e+00 0.00000000e+00\n",
" 2.75254528e+00 2.71770478e+00 2.71829649e+00 0.00000000e+00\n",
" 2.73534210e+00 2.71813892e+00 2.71828363e+00 2.71828179e+00"
]
}
],
"source": [
"n = 4\n",
"R = richardson(exp,1.0,0.1,n)\n",
"println(\"Richardson extrapolation on exp(x) at x = 1.0 :\")\n",
"show(stdout, \"text/plain\", R)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us compare with the exact value."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The table with errors :\n",
"4×4 Matrix{Float64}:\n",
" 1.40560126e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
" 6.91039636e-02 2.35219917e-03 0.00000000e+00 0.00000000e+00\n",
" 3.42634558e-02 5.77051997e-04 1.46637268e-05 0.00000000e+00\n",
" 1.70602718e-02 1.42912242e-04 1.80100989e-06 3.65210884e-08\n"
]
}
],
"source": [
"E = zeros(n, n)\n",
"for i=1:n\n",
" for j=1:i\n",
" E[i, j] = abs(R[i, j] - exp(1.0))\n",
" end\n",
"end\n",
"println(\"The table with errors :\")\n",
"show(stdout, \"text/plain\", E); println(\"\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Observe the exponents of the errors on the last row. At the left, we have an error $\\approx 10^{-2}$ and at the right the error $\\approx 10^{-8}$, four times less in magnitude."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 4. The Complex Step Derivative"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With complex arithmetic we can approximate derivatives accurately, avoiding cancellation errors caused by subtraction."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the function\n",
"\n",
"$$\n",
" g(z) = \\frac{\\exp(z)}{(\\cos(z))^3 + (\\sin(z))^3}.\n",
"$$\n",
" \n",
"With ``SymPy`` we can compute the derivatives of $g(z)$\n",
"(although the expressions are swelling)."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"using SymPy"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"g(z) = exp(z)/(sin(z)^3 + cos(z)^3)\n",
"dg(z) = (-3*sin(z)^2*cos(z) + 3*sin(z)*cos(z)^2)*exp(z)/(sin(z)^3 + cos(z)^3)^2 + exp(z)/(sin(z)^3 + cos(z)^3)\n",
"dg(1.0) = 1.64087713599607\n"
]
}
],
"source": [
"z = Sym(\"z\")\n",
"g(z) = exp(z)/((cos(z))^3 + (sin(z))^3)\n",
"dg = diff(g(z),z)\n",
"println(\"g(z) = \", g(z))\n",
"println(\"dg(z) = \", dg)\n",
"dg1 = dg(1.0)\n",
"println(\"dg(1.0) = \", dg1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In ``dg1`` we have the value of the derivative correct within the 64-bit floating-point precision. This value will be used to compare the accuracy of the numerical approximations."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"the complex step derivative\n",
"h = 1.00e-01 : 1.5914476349582889e+00 : 4.94e-02\n",
"h = 1.00e-02 : 1.6403947367308929e+00 : 4.82e-04\n",
"h = 1.00e-03 : 1.6408723131660323e+00 : 4.82e-06\n",
"h = 1.00e-04 : 1.6408770877678895e+00 : 4.82e-08\n",
"h = 1.00e-05 : 1.6408771355137928e+00 : 4.82e-10\n",
"h = 1.00e-06 : 1.6408771359912524e+00 : 4.82e-12\n",
"h = 1.00e-07 : 1.6408771359960252e+00 : 4.93e-14\n",
"h = 1.00e-08 : 1.6408771359960737e+00 : 8.88e-16\n",
"h = 1.00e-09 : 1.6408771359960739e+00 : 6.66e-16\n",
"h = 1.00e-10 : 1.6408771359960739e+00 : 6.66e-16\n"
]
}
],
"source": [
"i = complex(0.0,1.0)\n",
"csdg(a,h) = imag(g(a+i*h))/h\n",
"println(\"the complex step derivative\")\n",
"for k=1:10\n",
" hk = 10.0^(-k)\n",
" strhk = @sprintf(\"%.2e\", hk)\n",
" approx = csdg(1.0,hk)\n",
" strapp = @sprintf(\"%.16e\", approx)\n",
" error = abs(dg1 - approx)\n",
" strerr = @sprintf(\"%.2e\", error)\n",
" println(\"h = \", strhk, \" : \", strapp, \" : \", strerr)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We observe the $O(h^2)$ in the error as the exponents of the error are twice the exponents of $h$,\n",
"up to $k=8$, when the error is very close to the machine precision. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Moreover, in contrast to the finite differences, the error stays small, as $h$ goes smaller than $10^{-8}$."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.8.0",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}