Day 16, Mat360: Numerical Analysis

Today:

• Announcements:
  • I still haven't finished grading the "take-home" portion of your exam (although I have begun).

    I hope to get them done by Friday, so I don't have to do them over Spring break. By the way, have a nice break!

  • You have a reading assignment over break: please read chapter 5. It will form the basis for your first project, which will be to create your own very special computer-drawn signature. We won't cover it in detail, but I'll summarize what you need for this project.

• Let's recall what we did last time:
  • We went over the exam, and
  • we got started on polynomial interpolation.
  • Observations:
    • Interpolation is about fitting a few points exactly.
    • Muller's method's quadratic is easily fit using these special Lagrange polynomials, which we do using three points.
    • An alternative root-finding method, based on fitting more information at a point (three pieces of information, in particular \(f(a)\), \(f'(a)\), and \(f''(a)\), is a generalization of Newton's method, and leads to the scheme \[ x_{n+1} = x_n - \frac{ f'(x_n)}{f''(x_n)} + \text{sign}{(f'(x_n))} \frac{ \sqrt{f'(x_n)^2-2 f(x_n) f''(x_n)}}{f''(x_n)} \] Instead of using a tangent line, this method makes use of a tangent parabola to make its estimates for the next root (where the sign of \(f'(x_n)\) multiplies the square root term, to choose the root closest to \(x_n\).

• Today:
  • Back to the Lagrange interpolator.
    • Here is how our author lays out the problem:

      ---

      The main use for interpolating polynomials in numerical analysis is to approximate non-polynomial functions in the following way.

      Suppose we know the value of $f$ at a selection of points. That is, we know $f(x_{0})=y_{0},f(x_{1})=x_{1},\ldots,f(x_{n})=y_{n}$ and perhaps not much more. The interpolating polynomial of least degree \(P_n\) passing through the $n+1$ points \[ (x_{0},y_{0}),(x_{1},y_{1}),\ldots,(x_{n},y_{n}) \] will, by construction, agree with $f$ at $x_{0},x_{1},\ldots,x_{n}$ and we can say with some precision how closely this interpolating polynomial \(P_n\) agrees with $f$ at other points as well. The values of the interpolating polynomial at these "other points" are what we refer to as approximations of the non-polynomial function.

      Then \[ L_{n}(x)=\sum_{i=0}^{n}\frac{p_{i}(x)}{p_{i}(x_{i})}y_{i}, \] where \[ p_{i}(x)=\prod_{\substack{j=0\\j\neq i}}^{n}(x-x_{j}) = (x-x_{0})\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_{n}) \]

      "We will adhere to the practice of calling it the interpolating polynomial of least degree, or use the notation $P_{n}$, when the form is unimportant and will add the phrase Lagrange form, or use the notation $L_{n}$, when it is."

      ---

    • Start simple -- linear

      We will write the linear Lagrange interpolator as a sum of two linear functions, each of which "takes care of one point, and gets out of the way at the other". What does this mean? These "basis" functions will have that property that

      \(l(x,x_1,x_2)\) will be 1 at \(x_1\), and 0 at \(x_2\).

    • We can then recursively define the higher order Lagrange interpolators, based on these linear ones -- simple! How do we get quadratics that fit three points, while reusing our linear functions?

    • Lagrange's form is computationally silly, however: wasteful.

      Neville's method may be a better way to compute the interpolating polynomial values.

  • One of the most interesting facets of Lagrange interpolation is the form of the error term engendered by Lagrange interpolation.

    ---

    Then if $f$ has $n+1$ derivatives on $(a,b)$ and $f,f',f'',\ldots,f^{(n)}$ are all continuous on $[a,b]$, then there is a value $\xi_{x}\in(a,b)$ such that \[ f(x)-P_{n}(x)=\frac{f^{(n+1)}(\xi_{x})}{(n+1)!}(x-x_{0})(x-x_{1})\cdots(x-x_{n}). \]

    ---

  • One of the most important things to understand is the uniqueness of the interpolating polynomial of least degree:

    Theorem: The polynomial, $P_{n}$ , of least degree interpolating the data $(x_{0},y_{0}),(x_{1},y_{1}),\ldots,(x_{n},y_{n})$ exists and is unique. Moreover, any interpolating polynomial of degree at most $n$ is equal to $P_{n}$.

  • Why does Neville's method work? Let's work out the geometric details for the first few cases, \(n=1,2,3\), and hopefully we'll get the idea.

  • Sidi's method, which is another variation on Newton's method, is one that I've never encountered -- but looks like a really interesting and sensible idea:

    ---

    Sidi's method is summarized by the formula \[ x_{n+1}=x_{n}-\frac{g(x_{n})}{p'_{n,k}(x_{n})} \] where $p_{n,k}$ is the interpolating polynomial passing through the points \[ (x_{n},g(x_{n})),(x_{n-1},g(x_{n-1})),\ldots,(x_{n-k},g(x_{n-k})). \] When $k=1$, this is exactly the secant method. When $k=2$, this method uses the same parabola as does Müller's method, but in a different way. In Müller's method, the next iteration is found by locating a root of the interpolating polynomial. In this method, the next iteration is found by locating a root of a tangent line to the interpolating polynomial.

    ---

  • Here's my Mathematica code to solve several of our problems from the section.

• Links:
  • Our Textbook: Tea Time Numerical Analysis (Leon Q. Brin)
  • Homework solutions:
    • 1.1 (Mathematica)
    • 1.2 (Mathematica)
    • 1.3 (Mathematica)
    • 2.1 (Mathematica)
    • 2.2 (Mathematica)
    • 2.3 (Mathematica)
    • 2.4 (Mathematica)

  • My Mathematica code for Section 2.6. Muller code, too.
  • My Octave Project for 2.6 homework.
  • Octave project for 2.5
  • Divided Differences. I use these in my Muller code, and they will prove useful when it comes to interpolating polynomials.
  • Author's octave code for these complex fractal convergence diagrams
  • Newton's method in Mathematica
  • fixed point with Aitken's acceleration in Mathematica
  • Some nice notes (and code) snippets on fixed point iteration, including these acceleration techniques.
  • fixed point in Mathematica
  • Check out how my bisection works, and how it is applied to arcsine....
  • Various sins....
    • Sins in Mathematica
    • Sins in Octave
    • ArcSins in Mathematica