. Precision, Interpolation, and Numerical Differentiation

.	Finite and Infinite Precision

. . 	Exact Inverse of a Hilbert Matrix

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

\[LeftBracketingBar](Subscript[ℋ, 10])\[RightBracketingBar]

(Subscript[ℐ, 10]) = 
      Inverse[(Subscript[ℋ, 10])]

(Subscript[ℐ, 10]) . (Subscript[ℋ, 10])

% - IdentityMatrix[Length[%]]

. . 	Matrix Norm

Table[Random[], {5}, {5}]

\(TraditionalForm\`%\^2\)

Apply[Plus, %]

Apply[Plus, %]

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

\[LeftDoubleBracketingBar]ℳ_\[RightDoubleBracketingBar]:=Apply[Plus,Apply[Plus,(ℳ^2)]]

\[LeftDoubleBracketingBar]Table[
            0, {5}, {5}]\[RightDoubleBracketingBar]

\[LeftDoubleBracketingBar]IdentityMatrix[
            5]\[RightDoubleBracketingBar]

. . 	Numerical Inverse of a Hilbert Matrix

(Subscript[𝒩, 10]) = 
      N[(Subscript[ℋ, 10])]

Map[Precision, %, {2}]

(Subscript[ℐ, 10]) = 
      Inverse[(Subscript[𝒩, 10])]

(Subscript[ℐ, 10]) . (Subscript[𝒩, 10])

\[LeftDoubleBracketingBar]% - 
        IdentityMatrix[Length[%]]\[RightDoubleBracketingBar]

(Subscript[𝒩, 10]) = 
        N[(Subscript[ℋ, 10]), 30];

Map[Precision, %, {2}]

(Subscript[ℐ, 10]) = 
        Inverse[(Subscript[𝒩, 10])];

(Subscript[ℐ, 10]) . (Subscript[𝒩, 10])

\[LeftDoubleBracketingBar]% - 
        IdentityMatrix[Length[%]]\[RightDoubleBracketingBar]

.	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. Entering !!rt.data should allow you to view this file.

Data Structure

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

rθ = 
      ReadList["rt.data", {Number, Number}]

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

θs = N[θs °]

θr = {θs, rs}^T

..	Interpolation

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

𝓇 = Interpolation[θr]

𝓇[0.2]

ParametricPlot[{𝓇(θ) cos(θ), 𝓇(θ) sin(θ)}, {θ, 0, 
          2 π}, AspectRatio -> Automatic, 
        GridLines -> Automatic];

..	Area of a parametrically defined region

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

(∫_0^(f[θ]))r \[DifferentialD]r

From this result, the definition of the (numerical) area of a shape parametrized by f is

𝒜(𝒻_) := 
      (1/2) NIntegrate[((𝒻(θ))^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[θ]\)\ \(r\^2\) \(cos(θ)\)\ dr\ dθ\)\)\/𝒜\)\)

(∫_0^(f[θ]))(r^2) \[DifferentialD]r

(Subscript[x, CM])(𝒻_) :=  
      NIntegrate[((𝒻(θ))^3) cos(θ), {θ, 0, 2 π}]/(3  𝒜[𝒻])

(Subscript[y, CM])(𝒻_) := 
      NIntegrate[((𝒻(θ))^3) sin(θ), {θ, 0, 2 π}]/(3  𝒜[𝒻])

(Subscript[x, CM])(𝓇)

(Subscript[y, CM])(𝓇)

Options[NIntegrate]

NIntegrate[((𝓇(θ))^3) cos(θ), {θ, 0, 2 π}, AccuracyGoal -> 3]/(3  𝒜[𝓇])

.	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?

((f(h + x) - f(x))/h) + ((O(h))^3)

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.

((f(h + x) - f(x - h))/(2 h)) + ((O(h))^3)

Suppose that you can easily compute (or know) the values f[x], f[x + h] and f[x - h].  Form an arbitrary linear combination of these terms, say

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

ℓ + ((O(h))^2)

Solve[% == (f^′)(
          x), {α, β, γ}]

ℓ /. First[%] // Simplify

You should have found the best "two-point" approximation to (f^′)[x] which is (( f[h + x] -  f[x - h])/(2  h)).

((f(h + x) - f(x - h))/(2 h)) + ((O(h))^3)

exact = (sin^′)(x) /. 
        x -> 0.2

unsymmetric = ((sin(h + x) - sin(x))/h) /. {x -> 
            0.2, h -> 0.1}

twopoint = ((sin(h + x) - sin(x - h))/(2 h)) /. {x -> 0.2, h -> 0.1}

exact - unsymmetric

exact - twopoint

(ℓ = ((α f(x) + β f(x - h) + γ f(x + h) + δ f(x - 2  h) + ϵ f(x + 2  h))/h)) + ((O(h))^4)

Solve[% == (f^′)(
          x), {α, β, γ, δ, ϵ}]

ℓ /. First[%] // Simplify

% + ((O(h))^5)

(ℓ = ((α f(x) + β f(x - h) + γ f(x + h))/(h^2))) + O(h)

Solve[% == (f^″)(
          x), {α, β, γ}]

ℓ /. First[%] // Simplify

% + ((O(h))^3)

vals = Table[(x^3), {x, 0, 3, 0.1}]

((RotateLeft[vals] + RotateRight[vals] - 2  vals)/(0.1^2))

Drop[Rest[%], -1]

Table[Evaluate[(((∂^2) (x^3))/(∂(x^2)))],{x, 0, 3, 0.1}]

. . 	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, \(TraditionalForm
      \`y\_n = y\_\(n - 1\) + h\ \(f(x\_\(n - 1\), 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.

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

(Subscript[y, n])_ := (Subscript[y, n]) = 
        h f((Subscript[x, n]), (Subscript[y, n - 1])) + (Subscript[y, n - 1])

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

h = 0.05;

(Subscript[y, 0]) = 1;

f(x_, y_) = x (y^2);

Table[{(Subscript[x, n]), (Subscript[y, n])}, {n, 30}]

p(1) = 
        ListPlot[%, PlotJoined -> True];

Remove[f, y, h, Subscript]

DSolve[{y'(x)==x ((y(x))^2),y(0) == 1},y(x),x]

p(2) = 
        Plot[-(2/((x^2) - 2)), {x, 0, 1.4}, 
          PlotStyle -> Hue[1]];

Show[p(1), p(2)];

. . 	Taylor series expansion

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

𝒯 = 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)), etc, and so, for specific f we can express 𝒯 in terms of y(x) and x.

𝒯 = 
      Normal[𝒯]

f(x_, y_) := x (y^2)

𝒯=𝒯//.{y'(x)->f(x, y(x)),y''(x)->((∂f(x, y(x)))/(∂x))}

Note that, for higher orders, you would need to include analogous expressions for higher derivatives of y(x). To solve this equation numerically we replace y(x) by (Subscript[y, n]) and x by (Subscript[x, n]).

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

(Subscript[y, n])_ := (Subscript[y, n]) = 𝒯 /. {y(x) -> (Subscript[y, n - 1]), x -> 
              (Subscript[x, n - 1])}

h = 0.05;

(Subscript[y, 0]) = 1;

Table[{(Subscript[x, n]), (Subscript[y, n])}, {n, 30}]

p(3) = 
      ListPlot[%, PlotJoined -> True, 
        PlotStyle -> Hue[(1/3)]];

Show[p(1), p(2), p(3)];

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.

Remove[f, y, h, Subscript]

When you’ve finished...