. Precision, Interpolation, and Numerical Differentiation

Paul Abbott

Pre-requisites

To improve the appearance of graphics, TraditionalForm will be used as the default format type for text that appears in graphics,

$FormatType = TraditionalForm;

and the default style for text in graphics is set as follows.

$TextStyle = {FontSize -> 10};

The following code is used to display two approaches to numerical differentiation.

ln[δ_][x_, y_, m_] = 
        Line[{{x - (δ/(Sqrt[(m^2) + 1])), 
              y - ((δ m)/(Sqrt[(m^2) + 1]))}, {x + (δ/(Sqrt[(m^2) + 1])), y + ((δ m)/(Sqrt[(m^2) + 1]))}}];

NDPlot[f_, {x_, a_, b_}, n_:  5] := 
      Module[{δ = ((b - a)/n), dat, i1, xs, ys, pts}, 
        dat = Table[{x, f}, {x, a, b, δ}]; 
        i1 = Interpolation[dat, InterpolationOrder -> 1]; {xs, ys} = 
          Transpose[dat]; 
        pts = Drop[
            Rest[{xs, 
                  ys, ((RotateLeft[ys] - RotateRight[ys])/(2 δ))} // 
                Transpose], -1]; 
        Plot[{f, i1(x)}, {x, a, b}, PlotStyle -> {Hue[(2/3)], Hue[(1/3)]}, 
          Epilog -> {AbsolutePointSize[3], Hue[1], Point /@ dat, 
              GrayLevel[0], Apply[ln[δ], pts, 1]}]]

In this section, some examples of finite or machine arithmetic (see §§3.1.4, 3.1.6 of the Mathematica book) and infinite or arbitrary-precision numbers (see §§1.1.4, 3.1.5) are presented.

.. Exact Inverse of a Hilbert Matrix

A Hilbert matrix (see http://mathworld.wolfram.com/HilbertMatrix.html) is a popular example for testing numerical matrix methods because they are difficult to invert numerically.

(Subscript[ℋ, n])_ := 
      Table[(1/(i + j + 1)),  ... ]

.. Matrix Norm

(i)
(ii)

Congratulations, you have just computed the matrix norm of your matrix! Can you see how this worked?

\[LeftDoubleBracketingBar]ℳ_\[RightDoubleBracketingBar] := Apply[Plus, …]

.. Numerical Inverse of a Hilbert Matrix

Now we repeat the analysis of §.. numerically.

. Interpolation

.. Reading in data from an external file

Getting the Data

Get a copy of the file rt.data from the URL ftp://ftp.physics.uwa.edu.au/pub/Computational/rt.data. Save the resulting file as rt.data. You can then use Get File Path... under the Input menu to access this file. You can use !! (show file command, §A.2.11) to view this file.

Data Structure

You should be able to see the structure of the data file. (Hint: the data was taken by mounting an object on a turn-table and measuring the radius (r) every 5 degrees (θ)).

{rs, θs} = r(θ^T)

.. Interpolation

Interpolation can be used to give the radius, r, for any value of θ:

.. Area of a parametrically defined region

The area of an object whose radius r is expressed as a function of θ, f(θ), is given by \(TraditionalForm\`𝒜 = ∫\_0\%\(2  π\)\(∫\_0\%\(f(θ)\)r\ dr\ dθ\)\).

𝒜[f_] := 
      (1/2) NIntegrate[((f(θ))^2), {θ, 0, 
            2 π}]

.. Centre of mass

The x-coordinate of the centre of mass of an object whose radius r is expressed as a function of θ, is given by the following integral average:

\(TraditionalForm\`x\_CM = \(\(∫\_0\%\(2  π\)\(∫\_0\%\(f(θ)\)x\ r\ dr\ dθ\)\)\/\(∫\_0\%\(2  π\)\(∫\_0\%\(f(θ)\)\ r\ dr\ dθ\)\) = \(\(\(∫\_0\%\(2  π\)\(∫\_0\%\(f(θ)\)x\ r\ dr\ dθ\)\)\/𝒜\)\(.\)\)\)\)

. Numerical Differentiation

.. Numerical derivatives

To solve differential equations numerically, a method for computing derivatives numerically is required. How can the derivative of a function be computed numerically?

This result, with modification is the basis of techniques for integrating differential equations. For computing the derivative of a function f numerically, relying on this definition of differentiation is not optimal. As h -> 0, the error term in the Taylor series is first order in h.

Here is a plot of x (1 - x) (1 + 2  x) (—) sampled at a set of points (•) showing the numerical derivatives computed using ((f(x + h)  -  f(x))/h) (—) and ((f(x + h) -  f(x - h))/(2 h)) (—).

NDPlot[
        x (1 - x) (1 + 2  x), {x, -0.6, 1.1}];

Suppose that you can easily compute (or know) the values f(x), f(x + h) and f(x - h), i.e., the values of the function at three points, x and x ± h.

lc = ((α f(x) + β f(h + x)  + γ f(x - h) )/h)

You should have found the best “two-point” approximation to (f^′)(x) which is ((f(x + h) -  f(x - h))/(2 h)). Because α -> 0 this expression does not depend upon f(x).

.. Euler's method

Consider the completely general first-order differential equation, (y^′) = f(x, y). If you know y(x) at some initial point, say y((Subscript[x, 0])) = (Subscript[y, 0]), then you can compute f((Subscript[x, 0]), (Subscript[y, 0])) and the value of y((Subscript[x, 0]) + h) = (Subscript[y, 1]) can be obtained approximately from Taylor's theorem, i.e.,  (Subscript[y, 1]) = (Subscript[y, 1]) + h f((Subscript[x, 0]), (Subscript[y, 0])).  In general, (Subscript[y, n]) = 
        (Subscript[y, n - 1]) + h f((Subscript[x, n - 1]), (Subscript[y, n - 1])) where (Subscript[x, n]) = n h. This method is known as Euler's method. The error in each step is of ((O(h))^2).

x /: (Subscript[x, n])_ = n h;

y /: 
      (Subscript[y, n_ /; n > 0]) := y /: (Subscript[y, n]) = 
        (Subscript[y, n - 1]) + h f((Subscript[x, n]), (Subscript[y, n - 1]))

Consider the differential equation (y^′) = f(x, y) = x (y^2).

Conclude this section by clearing the values assigned to x, y, and h.

Clear[h, x, y]

.. Taylor series expansion

Euler's method only uses the first order term of the Taylor series.  The expansion to next order is

𝒯 = 
      Normal[y(h + x) + ((O(h))^3)]

From the differential equation, y'(x)=f(x, y(x)), we can also obtain exact expressions for higher derivatives, e.g., y''(x)=((∂f(x, y(x)))/(∂x)), Derivative[3][y](x)=(((∂^2) f(x, y(x)))/(∂(x^2))), etc., and so, for specific f, we can express 𝒯 in terms of y(x) and x.

The method of this section is an example of a numeric-symbolic hybrid. We used exact symbolic differentiation applied to the right hand side of our differential equation so that we could substitute this into the Taylor series of the solution to obtain a numerical method whose error depended on the order of the Taylor series expansion. Note that we could easily extend the method to higher orders.