Numerical Analysis

I was planning on covering all the topics in numerical methods in the Lecture notes though the pace of this class is about $10$ times faster than I expected, and so covering everything would be redundant and tedious. Thus, this post exists. Ideally I’ll be documenting about $3$ weeks of content into a single post, every $3$ weeks. Many of the topics in numerical methods tend to be very similar, which is probably partly responsible for the rapid pace. Recognising these similarities will make it unreasonably easier to digest the content. Rejoice classmates.

Taylor Polynomials

There’s lots of different notational ways to display taylor series. A good starting point would be

\begin{equation} f(x) \approx f(a)+f'(a)(x-a)+\frac{f''(a)}{2!}(x-a)^2+\dots+\frac{f^{(n)}(a)}{n!}(x-a)^n. \end{equation}

Notice that the first term $f(a)$ and $f'(a)$ follow the same pattern as the rest of the series when starting at $n=0,$ so it simplifies nicely while the third term $(n=2)$ keeps the redundant $2!$ to display the reoccuring pattern. Now writing this out can get tedious so we’ll use summation notation and recognise that eqref(1) is equal to

\begin{equation} \sum_{n=0}^{M} \frac{f^{(n)}(a)}{n!}(x-a)^n, \end{equation}

where $f^{(n)}(a)$ is the $n_{th}$ derivative of the function and $M \in \mathbb{Z}^+$ . With this, we can find polynomials of degree $n$ that approximate $f(x).$ To do this, we need to choose a function we want to approximate, some element $a\in \mathbb{R}$ $($ or $\mathbb{C},$ though I won’t be covering the complex case $)$ and the degree of the polynomial. We’ll denote the degree of the polynomial to be $P_{n}(x),$ hence, for example,

\begin{equation} P_{2}(x) = f(a)+f'(a)(x-a)+\frac{f''(a)}{2!}(x-a)^2. \end{equation}

Given the data, the standard is to first check the degree of the polynomial you want to use to approximate the function. Say we chose a degree of order 4, i.e. $P_4(x).$ Then we would look at the given function and get it’s derivatives up to four. That is, $f'(a),f''(a),f'''(a),f^{(4)}(a).$ Next we would plug in our values into $P_4(x),$ specifically, the value of $a$ and the derivatives. After evaluating $P_4(x)$ we note that

\begin{equation} P_4(x) \approx f(x) \end{equation}

Meaning we can approximate any value of $f(x)$ (that the function allows), such as $f(3)$ by computing $P_4(3)$ instead. And so we are done. If you happened to choose $a=0$ you may have noticed that it’s the Maclaurin Series. Looking back at eqref(1), there is a small subletly in the $\approx$ symbol. If the function we are “approximating” is a polynomial, then we can always get the exact function after finitely many steps; this isn’t the case for non-polynomial functions such as $e^x$ for example.

Taylor’s error

Recall that

\begin{equation} f(x)= \sum_{n=0}^\infty \frac{f^{(n)}(a)}{n!}(x-a)^n. \end{equation}

This gives us the exact function, though we can’t compute an infinite amount of terms on a computer. So how do we get know when our approximation has enough accuracy? though without an error term, we don’t really know how well our polynomial approximates the point of a function. Thus we have $R_n(x),$ the remainder/error term

\begin{equation} R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}(x-a)^{n+1} \end{equation}

where $f^{(n+1)}(\xi)$ is the $(n+1)th$ derivative of $f$ at $\xi.$ Also $x\leq \xi \leq a .$ Equivalently $\xi$ is a point $\xi \in [x,a]$ . $\newline$ At first glance this formula might seem like it’s come out of nowhere, but it’s a term added onto $P_n(x)$ following the same pattern which can be seen below.

\begin{equation} \begin{aligned} f(x)={}&f(a)+f'(a)(x-a)+\frac{f''(a)}{2!}(x-a)^2+\dots+\frac{f^{(n)}(a)}{n!}(x-a)^n \\ &+\frac{f^{(n+1)}(\xi)}{(n+1)!}(x-a)^{n+1} \\ ={}&\sum_{k=0}^n \frac{f^{(k)}(a)}{k!}(x-a)^k + R_n(x). \end{aligned} \end{equation}

If you’re familiar with factorials, you may have realised that $R_n(x)$ can be written in the terms of the gamma function:

\begin{equation} R_n(x)=\frac{f^{(n+1)}(\xi)}{\Gamma(n+2)}(x-a)^{n+1}, \end{equation}

because $\Gamma(n+2)=(n+1)!$ for $n \in \mathbb{Z}^+$ . Why would you want to use this? I don’t know, but you can if you want. Actually obtaining the exact value of $R_n(x)$ pretty much never happens in practice, rather, we use it to find the maximum error $|R_n(x)|,$ by choosing the smallest value of $\xi$ possible, but since we don’t know $\xi$ we can’t do that, which is why the interval $\xi \in [x,a]$ is how we end up choosing a suitable $\xi.$

Lagrange polynomials

In Lagrange form, we can find a polynomial $l_k(x)$ of degree $n$ by

\begin{equation} l_k(x) = \frac{(x-x_0) \dots (x-x_{k-1})(x-x_{k+1}) \dots (x-{x_k})}{{(x_k-x_0) \dots (x_k-x_{k-1})(x_k-x_{k+1}) \dots (x_k-x_n)}} \end{equation}

so that

\begin{equation} l_k(x_i) = \begin{cases} 0 &\text{if } i \neq k \\ 1 &\text{if } i=k \end{cases}. \end{equation}

The interpolating polynomial

\begin{equation} \begin{aligned} P_n(x)={}&l_0(x)y_0+l_1(x)y_1+ \dots + l_k(x)y_k + \dots + l_n(x)y_n \\ &= \sum_{k=0}^n l_k(x)y_k. \end{aligned} \end{equation}

Note that $P_n(x_j)=y_j.$ because $P_n(x_j)=y_j = 0+ \dots + 0 + l_j(x_j)y_j+0+ \dots + 0 = y_j.$ Now to find the interpolating polynomial we need $x$ and $y$ values. For some questions you won’t get both $x$ and $y$ values, but rather a function and just $x-$ values, for those cases, you evaluate the functions by plugging in each $x_i$ value to get all $y_i$ values before continuing. Now the general method for finding $P_n(x)$ using lagrange polynomials can get tedious by hand, luckily there is a pattern which can speed up the process. Setting $\pi_n(x) = (x-x_0)(x-x_1) \dots (x-x_n)$ and $\pi_{n}'(x_k) = \prod_{i=0,i\neq k} (x_k-x_i),$ for $k=0,\dots,n,$ let’s us rewrite $l_k(x)$ as $\frac{\pi_n(x)}{(x-x_k)\pi_n'(x_k)}.$ This is nice because $(x-x_k)$ in the denominator of $l_k(x)$ always cancels out. So $l_0(x),l_1(x),\dots, l_n(x)$ will always be of the form

\begin{equation} \begin{gathered} l_0(x)=\frac{(x-x_1)(x-x_2)(x-x_3)(x-x_4)\dots(x-x_n)}{\pi_n'(x_0)} \\[6pt] l_1(x)=\frac{(x-x_0)(x-x_2)(x-x_3)(x-x_4)\dots(x-x_n)}{\pi_n'(x_1)} \\[6pt] l_2(x)=\frac{(x-x_0)(x-x_1)(x-x_3)(x-x_4)\dots(x-x_n)}{\pi_n'(x_2)} \\[6pt] l_3(x)=\frac{(x-x_0)(x-x_1)(x-x_2)(x-x_4)\dots(x-x_n)}{\pi_n'(x_3)} \\[6pt] \vdots \\[6pt] l_n(x)=\frac{(x-x_0)(x-x_1)(x-x_2)(x-x_3)(x-x_4)\dots(x-x_{n-1})}{\pi_n'(x_n)}. \end{gathered} \end{equation}

This let’s us obtain $P_n(x)$ without getting lost in the details.

Newtons divided differences

This part is a bit boring and tedious to write all out, so I’ll explain a geometric way of thinking about it, if you care for the details Wikipedia gives the gist. We want to create a divided difference table to evaluate $P_n(x)$ by using a set of points, while making sure our $P_n(x)$ fits through the points.

\begin{equation} \begin{array}{|c|c|c|c|c|c|c|} \hline i\vphantom{\vdots} & x_i\vphantom{\vdots} & y_i\vphantom{\vdots} & D_1\vphantom{\vdots} & D_2\vphantom{\vdots} & \cdots\vphantom{\vdots} & D_n\vphantom{\vdots} \\ \hline 0\vphantom{\vdots} & \bullet & \bullet & \times & \times & \cdots & \times \\ \hline 1\vphantom{\vdots} & \bullet & \bullet & \bullet & \times & \cdots & \times \\ \hline 2\vphantom{\vdots} & \bullet & \bullet & \bullet & \bullet & \cdots & \times \\ \hline 3 & \bullet & \bullet & \bullet & \bullet & \dots & \vdots \\ \hline \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ \hline n\vphantom{\vdots} & \bullet & \bullet & \bullet & \bullet & \cdots & \bullet \\ \hline \end{array} \end{equation}

Suppose we were given the points $x=(0,1,2,3,4,5)$ and $y=(0,1,4,15,40,85).$ So $x_0=0,x_1=1,x_2=2,$ etc. and $y_0=0,y_1=1,y_2=4,$ etc. Then using the newtons’s divided difference table we could find a polynomial which passes through the points $(x,y)=((0,0),(1,1),(2,4),(3,15),(4,40),(5,85)).$ There are many different ways to reach this polynomial. One way is to just compute

\begin{equation} \begin{aligned} P_n(x) ={}& f[x_0] \\ &+ (x-x_0)f[x_0,x_1] \\ &+ (x-x_0)(x-x_1)f[x_0,x_1,x_2] \\ &+ \cdots \\ &+ (x-x_0)(x-x_1)\cdots(x-x_{n-1})f[x_0,x_1,\dots,x_n]. \end{aligned} \end{equation}

Though as you can see, this can get quite messy, and while we know the values of $f[x_0],$ we still don’t know $f[x_0,x_1], f[x_0,x_1,x_2]$ etc., so before we can compute $P_n(x)$ we need to find those values. Though the notation is getting a bit annoying. We can simplify this by assigning each $f[x_0,\dots,x_n]$ with a variable $a$ , so that $f[x_0] = a_0,$ $f[x_0,x_1]=a_1,$ $f[x_0,x_1,x_2]=a_2$ and so forth:

\begin{equation} \begin{aligned} P_n(x) &= a_0 +(x-x_0)a_1 +(x-x_0)(x-x_1)a_2 + \dots +(x-x_0)(x-x_1)\dots(x-x_{n-1})a_n \\ &= a_0 + \sum_{i=1}^{n} a_i \prod_{j=0}^{i-1}(x-x_j) \end{aligned} \end{equation}

Below is an example of how I would geometrically solve the divided difference table, to avoid the influx of notation.The hollow-yellow boxes are the values $a_0,a_1,\dots,a_6$ respectively. The full-blue box is the value currently being evaluated. The full-yellow boxes are the values used to compute the full-blue value. Explaining the procedure is pretty awkward so just look at the animation and see if you can see what I’m doing. Otherwise solve the newton’s divided difference table however you like, there’s plenty of different setups and tables you can find online that do this.

Now that the table has been complete, we know that $a = (0,1,1,1,0,0),$ so we can plug these values along with all the $x_i$ values into $P_n(x)$ to get the polynomial. In general you won’t get nice values like the one in the animation, those values were chosen specifically for simplicity sake.You may be worried about what to do if the procedure leads to a division by zero. For such cases you take the derivative of the function instead of the standard algorithmic process. This means that just the $(x,y)$ points aren’t enough for such cases, since you won’t be able to get a derivative of a function with just points. In fact, hermite interpolation is a special case of this.

Hermite Interpolation

Given the function $f(x)=x^4+x^3-x$ and the points $(x,y)=((-1,1),(0,0),(1,1)),$ find $P_5(x)$ . Below can be used for reference for how I would complete the divided difference table.

Notice that for hermite interpolation, it’s pretty much the exact same as the above example, just that we need a function to compute the derivatives because we have to repeat each $x_i$ and $y_i$ twice in their respective columns. This is because the hermite polynomial approximation is

\begin{equation} P_{2n+1}=a_0+\sum_{i=1}^{2n+1}a_i \prod_{j=0}^{i-1}(x-x_j). \end{equation}

Similarly to lagrange’s error term, the hermite polynomials error term is

\begin{equation} R_{2n+1}(x)=f(x)-P_{2n+1}(x)=(\pi(x))^2\frac{f^{(2n+1)}(\xi)}{(2n+1)!}. \end{equation}

Again, you can also write $(2n+1)!$ in terms of the gamma function if you want to for whatever reason. $\pi(x)$ here again refers to $(x-x_0)\dots(x-x_n).$ To find a bound using the error term, you first find the hermite polynomial, and plot it specify the max for the function, then use that to get the an inequality, that is the bound of the error, similar to how we obtain Taylor’s error bound. $\newline$ And that’s it for now. Bye.