Numerical Analysis: A Nonlinear Story

Shortcut to this page: ntrllog.netlify.app/nonlinear

Notes provided by Professor Li-Tien Cheng

This is the nonlinear version of the linear algebra story. The linear story was looking at how to solve problems of the form `Ax = b` i.e. `Ax - b = 0`, which is linear. This story looks at how to solve problems of the form `f(x) = g(x)` where `f(x), g(x)` are nonlinear. For example, they could be `x^2, x^3, cos(x), e^(x), ...`

A Numbers Store Stores Numbers

How does a computer store numbers? Vaguely, we could describe it as putting each digit in a space somewhere in the computer's memory. So `100` would need `3` spaces`color(red)(text(*))`. This seems like an easy and obvious idea until we get to numbers like `pi`, which has an infinite number of digits. Does this mean we need an infinite amount of memory to store `pi`? But computers don't have infinite memory, so what happens then?

`color(red)(text(*)`Computers work in binary, so things are a little bit different than this.

Computers store a number by keeping as much of it as they can and then throwing away the rest. So if a computer can only store at most `8` digits, then it would store `pi` as `3.1415926`. But complications arise when we're dealing with really big or really small numbers. For example, `0.00000000000000045` would be stored as `0.0000000`, which is just `0`. This is a problem because `0.00000000000000045` isn't equal to `0`, but it would be treated as such. A similar problem occurs for `1000000000000000`, which would be stored as `10000000`, which is much much smaller.

This suggests needing a smarter way to represent the numbers. Enter scientific notation. With scientific notation, we can represent `0.00000000000000045` as `0.45 xx 10^(-16)` and `1000000000000000` as `1 xx 10^15`. Then we would only need to store the `45` and the `-16` for `0.45 xx 10^(-16)`. And the `1` and the `15` for `1 xx 10^15`.

Scientific notation helps make it easier to store numbers, but some complications still pop up. Consider a computer that can only store at most `3` digits in the mantissa. Let's say we wanted to store the number `0.002146`. In scientific notation, it would be `0.2146 xx 10^(-2)`, which has `4` digits in the mantissa. Our computer can only store `3` digits in the mantissa, so how should it be stored? There are 2 options: chopping and rounding.

Under chopping, `0.2146 xx 10^(-2)` would become `0.214 xx 10^(-2)`. As the name implies, the rest of the digits (after the `3^(rd)`) are "chopped off". Under rounding, `0.2146 xx 10^(-2)` would become `0.215 xx 10^(-2)`.

From this, we can see that what the computer stores isn't always going to be exactly equal to the original number. The version of the number that the computer stores can be denoted by `fl(x)`, where `x` is the original number.

So the difference between the original number and the stored number is going to be the error. There are `2` types of errors:

• absolute error ` = abs(x - fl(x))`

• relative error ` = abs(x - fl(x))/abs(x)`

Errors are bad, but unavoidable. So it would be helpful to know how bad the relative error could be for any number we want to store.

Let's say we want to store the number `x` (which is not `0`) on an `n`-digit chopping computer (a computer that stores at most `n` digits and uses the chopping method). Written in scientific notation:

Since this computer uses chopping, it would only store the first `n` digits of the mantissa.

We want to compute the relative error, so we need `abs(x - fl(x))` and `abs(x)`.

To find out the worst our error could be, we want to make it as large as possible. We do that by making the numerator as big as possible (by letting `d_(n+1) = d_(n+2) = ... = 9`) and making the denominator as small as possible (by letting `d_1 = 1` `color(red)(text(*))` and `d_2 = d_3 = ... = 0`).

`color(red)(text(*))``d_1` cannot be `0` because of the way numbers are written in scientific notation. We wouldn't write something like `0.083 xx 10^3`; it would be `0.83 xx 10^2`.

`color(red)(text(*))``0.999... = 1`

So the relative error — which we can denote by `epsilon` — can be at most `10^(1-n)`. This is the largest error possible for an `n`-digit chopping computer.

Errors in Arithmetic Operations

Let's suppose we wanted to add `0.42312` and `0.63546` on a `3`-digit rounding computer. First, the computer would store them as `0.423` and `0.635`. Then, the computer would add those two to get `1.058``color(red)(text(*))`. Since this is a `3`-digit rounding computer, it would say that the answer is `1.06`. In math notation:

Notice how this differs from the actual answer, which is `1.05858`.

`color(red)(text(*))` This is a `4`-digit number, but the computer can only store `3`-digit numbers. It turns out that computers temporarily store results of arithmetic operations using more digits if needed.

We've established that there are errors in storing numbers. And performing operations with "wrong" numbers will lead to numbers that are "more wrong". But how much error could there be when doing arithmetic operations?

Let's suppose we wanted to multiply two numbers, `x` and `y`. In the computer, they're stored as `fl(x)` and `fl(y)`. Well, we can write:

In words: `fl(text(some number))` is the actual value with some error attached to it.

So multiplying those two values would result in:

So how big is the relative error?

In the previous section, we defined `epsilon` to be the largest possible relative error. So all of these relative errors are smaller than `epsilon`, i.e. `abs(delta_1), abs(delta_2), abs(delta_3) le epsilon`

So multiplying two numbers will lead to error, but it will be at most `3` times as bad as the worst error in storing a number.

Now, let's suppose we wanted to subtract two numbers `x` and `y`.

But `abs(x/(x-y))` and `abs(y/(x-y))` can be very large when `x` is close to `y`. When this happens, `abs(delta)` will be very large. This is called catastrophic cancellation.

1) Errors become slightly larger every arithmetic operation.

2) Catastrophic cancellation can occur when subtracting numbers.

Solving Nonlinear Equations: Roots

How do we solve something like `cos(x) = x`? Let's say we didn't have a calculator, so we can't find the exact answer. How can we find a (good) approximation? In general, how do we find an approximation for `f(x) = g(x)`?

Getting to the Root of the Problem

One way is to change the equation to the form `f(x) = 0`. Then, finding the roots of `f(x)` will be the solutions. Of course, it isn't going to be obvious what the roots are just from looking at something like `cos(x) - x = 0`. But we can try to trap a root in a small interval (i.e. find a small interval between which a root lies).

This can be done with an iterative algorithm. We start with a root trapped in a large interval and reduce the size of the interval at every step. But how do we find an interval in which there is a root? Without knowing this, we can't start the algorithm!

We do this by using the Intermediate Value Theorem, which says if `f` is a continuous function in `[a,b]`, then `f` assumes every value between `f(a)` and `f(b).`

In this example, `f` achieves every value between the dotted lines between `a` and `b`.

So how does this help us find roots? We can apply it to the situation where `f(a)` and `f(b)` have different signs (e.g. one is positive and the other is negative). If `f(a)` and `f(b)` have different signs, then `f` has a root in `[a, b]`.

Again, `f` achieves every value between the dotted lines between `a` and `b`. Here, the `x`-axis lies between the dotted lines, so it is guaranteed that `f` will cross the `x`-axis, i.e. `f(c) = 0` at some point `c` between `a` and `b`.

This suggests that we start the algorithm by picking `a_0, b_0` such that `f(a_0), f(b_0)` have different signs (e.g. `f(a_0)` is positive and `f(b_0)` is negative).

Now we have to reduce the size of the interval. We can split it in half; then, the root has to either be on the left half or the right half. So let's split it at `c_0` (the midpoint) and suppose that `f(c_0)` is positive.

Now we have two choices for what our new interval `[a_1, b_1]` could be: `[a_0, c_0]` or `[c_0, b_0]`. The IVT doesn't tell us if there's a root between `a_0` and `c_0`, but it does guarantee that there is a root between `c_0` and `b_0`. So the new interval `[a_1, b_1]` should be `[c_0, b_0]`.

And this process of dividing the interval in half repeats until the interval is small enough. This process is called the bisection method.

Error of the Bisection Method

The thing is, we're looking for a number, but we end up with an interval (two numbers). So how do we get our answer? We could just pick some number between that interval. At this point, the interval is really small, so our approximation should be close to the actual solution. But how much error is there?

Let's pick some number `d_n` and use that as our approximation. In this picture, there would be the most error if the actual solution was at `b_n`. In general,

We usually pick as our approximation the midpoint of `a_n` and `b_n`. That actually gives us a better approximation because the most error we could have then is half of the previous error. (This would be when the exact solution is either at `a_n` or at `b_n`.)

`abs(text(exact)-text(approximate)) = abs(r - c_n) le abs(b_n - a_n)/2`

So to summarize, the algorithm for finding a root of a polynomial is:

1. Start with an interval `[a_0,b_0]` where `f(a_0)` and `f(b_0)` have different signs

2. Cut the interval in half at `c_0`, which is the midpoint of `[a_0,b_0]`

3. Pick a new (smaller) interval such that either `f(a_0)` and `f(c_0)` have different signs or `f(c_0)` and `f(b_0)` have different signs

4. Repeat the steps for that new (smaller) interval

Convergence of the Bisection Method

The bisection method generates a sequence of approximations:

If the actual solution is `r`, then what we want to say is that `color(blue)(lim_(n->oo) c_n = r)`. In words, the more times we do the algorithm, the more approximations we will get and the closer those approximations will be to the actual solution `r`.

Another way to say that is that the difference (i.e. error) between the approximation we get and the actual solution gets closer and closer to `0`. Generally, it might not be guaranteed to happen, but for the bisection method, it is true that the more times we do the algorithm, the closer it will get to `0`.

Previously, we found that:

Well that means:

But what is `abs(b_1 - a_1)`? Thinking back to the steps of the algorithm, `[a_1,b_1]` represented the new interval after splitting `[a_0,b_0]` in half. So `[a_1,b_1]` is really half of `[a_0,b_0]`. Which means

Putting all that together:

Generalizing that idea, we can see that:

Now, as `n` gets larger and larger (i.e. as we do more steps of the algorithm), the right-hand side gets smaller and smaller. When `n` goes to `oo`, the right-hand side goes to `0`.

This proves that what we want to say is true: `color(blue)(lim_(n->oo) abs(r - c_n) = 0 => lim_(n->oo)c_n = r)`

When this happens, we say that the bisection method converges.

So in theory, the error gets really small when `n` is large. But how large should `n` be? How many times should we keep dividing intervals in half until we get a satisfactory answer?

If we know that the worst error is `1/2^(n+1) (b_0 - a_0)`, then we can decide how small the error could be.

function [c] = bisection(a,b,n)
    c = (a+b)/2;
for step = 1:n
    % f is defined as another function in another .m file
    if f(a)*f(c) <= 0
        b = c;
    else
        a = c;
    end
    c = (a+b)/2;
 end
end

Speed of the Bisection Method

Since we're dealing with nonlinear equations, we can say that computing `f(x)` is relatively hard (i.e. time-consuming). Of course, a computer can do it almost instantly, but it takes more time to compute powers and roots than additions and divisions. So we can measure the speed of the bisection method by counting how many times we compute `f(x)`.

In the code above, we're calcuating `f(a)` and `f(c)` every iteration. That's `2` costly operations. It turns out we can reduce it to only `1` costly operation per iteration if we store their values.

function [c] = bisectionOptimized(a,b,n)
    c = (a+b)/2;
fa = f(a);
for step = 1:n
    fc = f(c);
    if fa*fc <= 0
        b = c;
    else
        a = c;
        fa = fc;
    end
    c = (a+b)/2;
end
end

We only need to calculate `f(a)` once at the beginning to use IVT. After that, `a` becomes `c`. At this point, we've also already calculated `f(c)`, so when we're trying to find new intervals, we don't need to calculate `f(c)` again.

Straight to the Point

Another method we could use is Newton's method. In this method, we use tangent lines to approximate the root.

The idea is to start at an initial guess `x_0` and repeatedly draw tangent lines until we (hopefully) get close to the root.

Using the point-slope formula for the equation of a line, we can get the equation of the tangent line at `x_0`:

Our goal is to get `y = 0`, so we have:

Solving for `x`, we get:

This is the intuition behind the iterative formula for Newton's method:

There are some cases where Newton's method won't work.

Here, each approximation results in a tangent line that just keeps going to the right.

Here, the approximations will just keep bouncing back and forth between `x_0` and `x_1`.

Here, the tangent line doesn't cross the `x`-axis, so we can't get an approximation.

So unlike the bisection method, Newton's method won't always work. However, it's still worth looking at Newton's method, because it turns out that it is generally faster than the bisection method.

Error of Newton's Method

To look at the error, we need an equation with `abs(x_(n+1) - r)`. Since `x_(n+1) = x_n - f(x_n)/(f'(x_n))`, we need `x_n`, `f(x_n)`, and `f'(x_n)` to get `x_(n+1)`. We can get those `3` parts by looking at the Taylor series expansion of `f`.

Notice that the first `2` terms of the Taylor series are the equation of a line for `f(x)` (and thus, are a linear approximation of `f(x)`). Since we're looking at approximations using lines, we'll consider the first `2` terms and collect the rest of the terms as some error.

Since we're interested in finding the root of `f(x)`, we can use `x = r` where `r` is a root to get:

Well, `x_n - f(x_n)/(f'(x_n)) = x_(n+1)` from the iterative formula for Newton's method, so we have:

Then, we can let `e_n = x_n - r` be the error between the approximation and the root (so `e_(n+1) = x_(n+1)-r`).

This means that the error between our `(n+1)^(st)` approximation and the root is roughly equal to the square of the error between the `n^(th)` approximation and the root. So if the error was small for the `n^(th)` approximation, the error will be even smaller for the `(n+1)^(st)` approximation.

Speed of Newton's Method

The iterative formula for Newton's method is `x_(n+1) = x_n - f(x_n)/(f'(x_n))`. From this, we can see that it requires `2` function evaluations per iteration. The bisection method required `2` function evaluations at the beginning, but then it only required `1` function evaluation per iteration. This means that Newton's method does roughly twice as many function evaluations per iteration as the bisection method. Doing more function evaluations means that Newton's method is slower than the bisection method.

Convergence of Newton's Method

But maybe Newton's method will converge faster than the bisection method will? Let's look back at the error equation `e_(n+1) = (f''(xi_n))/(2f'(x_n))(e_n)^2`. If we start with an initial guess close enough to `r`, then `(f''(xi_n))/(2f'(x_n)) ~~ (f''(r))/(2f'(r))`, which is a constant, so we can call it `c`. Then the error for Newton's method can be represented as `e_(n+1) ~~ c(e_n)^2`.

From this, we can see that Newton's method (generally) converges faster than the bisection method. The convergence of Newton's method is described as "quadratic convergence" because of the power `2`.

But this doesn't mean that Newton's method always converges faster than the bisection method. It's possible that Newton's method won't converge at all (when `f' = 0` or `f''` doesn't exist). However, there are some theorems that guarantee the convergence of Newton's method.

Let's suppose `f'(x)`,`f''(x)` are continuous around `x = r` and `f'(r) != 0`. Then it turns out that if we start with an initial guess close enough to `r`, we can get `x_n -> r` (i.e. the `n^(th)` approximation will approach `r`).

Another convergence result is that if `f''(x) > 0` (`f` is concave up) in the interval `[r, oo)` and `f'(x) > 0` (`f` is increasing) in the interval `[r,oo)`, then `x_0 in [r,oo)` implies Newton's method will converge to `r`.

This is an example where `f'' gt 0` and `f' gt 0`.

We can see that if we continue, we will eventually hit `r`. This is what the second convergence result is saying. It's important to note that we have to start at an `x_0` to the right of `r`, i.e. `x_0 gt r`.

We proved that Newton's method converges on the interval `[r,oo)` when `f'' gt 0` and `f' gt 0`. We can take a similar approach to prove that Newton's method converges:

• on the interval `[r,oo)` when `f'' lt 0` and `f' lt 0`

• on the interval `(-oo,r]` when `f'' gt 0` and `f' lt 0`

• on the interval `(-oo,r]` when `f'' lt 0` and `f' gt 0`

Oh See Can't You See

So Newton's method converges faster than the bisection method does, but it requires more function evaluations. Is there a method that converges faster than the bisection method and requires less function evaluations than Newton's method does?

Newton's method looked at tangent lines. Let's look instead at secant lines.

The secant method is a variation of Newton's method. It still uses the iterative formula `x_(n+1) = x_n - f(x_n)/(f'(x_n))`, but the difference is that `f'(x_n)` will be the slope of the secant line, not the tangent line.

The derivative can be calculated by `f'(a) = lim_(x->a) (f(x)-f(a))/(x-a)`. Without the limit, it can be approximated by `f'(a) ~~ (f(x)-f(a))/(x-a)`.

One way to get `f'(x_n)` is to consider the slope of the secant line through `x_(n-1)` and `x_n`. Then we have `f'(x_n) ~~ (f(x_n)-f(x_(n-1)))/(x_n-x_(n-1))`. This leads to the iterative formula for the secant method:

The idea is that we're getting our next approximation by drawing the secant line through our previous two approximations.

Notice how using the previous two approximations means we only need `1` function evaluation per iteration.

Error of Secant Method

To look at the error, we can start by rewriting the iterative formula:

Now we can look at `e_(n+1) = x_(n+1) - r`:

Using the fact that `e_(n-1) = x_(n-1)-r` and `e_n = x_n - r`:

Multiplying by `1 = (x_n-x_(n-1))/(x_n-x_(n-1))`:

Combining the results above, we get that `x_(n+1)-r = e_(n+1) ~~ Ce_(n-1)e_n` where `C = (f''(r))/(2f'(r))`. If we want to compare this error to the error in Newton's method (which was `e_(n+1) = c(e_n)^2`), then we should rewrite `e_(n+1) = Ce_(n-1)e_n` in the form of `A(e_n)^(alpha)` and find out what `alpha` is. If `alpha lt 2`, then the secant method converges slower than Newton's method. If `alpha > 2`, then the secant method converges faster than Newton's method.

Suppose `e_(n+1) = Ce_(n-1)e_n = A(e_n)^(alpha)`.

From `e_(n+1) = A(e_n)^(alpha)`, we get:

Rearranging terms, we get:

Plugging `e_(n-1)` and `e_n` into `Ce_(n-1)e_n = A(e_n)^(alpha)`:

Using `e_n = A(e_(n-1))^(alpha)`, we can find `(e_(n-1))^(alpha)`:

Now we substitute `e_n/A` for `(e_(n-1))^(alpha)` in `C(e_n)^(1/alpha)/A^(1/alpha)A(e_(n-1))^(alpha)`:

So `Ce_(n-1)e_n = C/A^(1/alpha)(e_n)^(1+1/alpha) = A(e_n)^(alpha)`, which means `alpha = 1 + 1/alpha`

Using the quadratic formula, we get:

So `x_(n+1)-r = e_(n+1) = A(e_n)^(1.62)`, which means the secant method converges slower than Newton's method.

Speed of Secant Method

Recall that the secant method only requires `1` function evaluation per iteration, and that Newton's method requires `2` function evaluations per iteration. This means that we can perform `2` steps of the secant method in the time it takes to perform `1` step of Newton's method. So to compare them on a level playing field, we can compare the error when performing `2` steps of the secant method vs. the error in performing `1` step of Newton's method.

For the secant method:

Comparing `Cabs(e_n)^2.62` to Newton's `cabs(e_n)^2`, we can see that the secant method can be considered faster.

Solving Nonlinear Equations: Fixed Point Method

Switching gears, let's look at ways to solve nonlinear equations without finding roots.

A Fixed Point of View

With Newton's method, we had: `x_(n+1) = x_n - f(x_n)/(f'(x_n))`. Well, we could consider the right hand side as just one big function of `x_n`, call it `F(x_n)`. Then we have `x_(n+1) = F(x_n)`. So we plug in `x_n` and get out `x_(n+1)`. Let's suppose that we do this enough times so that `x_n -> s`. Then:

This is saying that if we plug in a value, we get back that same value. This leads to another way of solving nonlinear equations: solving `x = F(x)` (as opposed to solving `f(x) = 0`).

The points `x` that satisfy `F(x) = x` are called fixed points. The idea is that if we apply `F` to `x` and we get back the same `x`, then that `x` is not changing, i.e. fixed.

For this method to work, we need `F` to be contractive. `F` is contractive if there exists a positive `lambda lt 1` such that `abs(F(x)-F(y)) le lambdaabs(x-y)` for all `x,y`. It's "contractive" because the distance between `x` and `y` (represented by `abs(x-y)`) is getting smaller after applying `F` (because `lambda lt 1`).

Note: If `F` is contractive, then `F` is continuous.

So how do we know what to use for `lambda`? It turns out that we can answer this by looking at the Mean Value Theorem, which says that if `f` is continuous on `[a,b]` and differentiable on `(a,b)`, then there exists a point `c` between `a` and `b` such that `f'(c) = (f(b)-f(a))/(b-a)`.

For our situation, we can adjust that equation a little bit to get:

Notice how similar this is to the contractive property: `abs(F(x)-F(y)) le lambdaabs(x-y)`. This means we can use `F'(xi)` for `lambda`.

This is the graphical interpretation of looking at fixed points. Notice that the point is on the line `y = x` and `y = F(x)`. So `x = F(x)`.

Convergence of the Fixed Point Method

Let's suppose the fixed point iterations converge to `s`. To look at the convergence, let's look at how the error compares to the error of Newton's method.

So let's start with `e_(n+1) = x_(n+1) - s`. From the formula `x_(n+1) = F(x_n)` and the definition of fixed point, we get:

To introduce `x_n - s = e_n`, we can look at the Taylor series expansion of `F(x_n)`:

So `e_(n+1) = F'(s)(e_n) + (F''(s))/2(e_n)^2 + ...` Notice that the error depends on the differentiability of `F`.

• `color(blue)(F'(s) != 0 => e_(n+1) ~~ F'(s)e_n)`

• `color(green)(F'(s) = 0\ text(and)\ F''(s) != 0 => e_(n+1) ~~ (F''(s))/2(e_n)^2)`

• `color(purple)(F'(s) = F''(s) = 0\ text(and)\ F'''(s) != 0 => e_(n+1) ~~ (F'''(s))/(3!)(e_n)^3)`

Recall that the error for Newton's method is `e_(n+1) = c(e_n)^2` where `c` is a constant.

• Fixed point method vs. Newton's method

• `color(blue)(e_(n+1) ~~ F'(s)e_n\ text(vs.)\ e_(n+1) = c(e_n)^2)`

• `color(blue)(text(Fixed point converges slower than Newton's method does))`

• `color(green)(e_(n+1) ~~ (F''(s))/2(e_n)^2\ text(vs.)\ e_(n+1) = c(e_n)^2)`

• `color(green)(text(Fixed point converges at the same rate as Newton's method does))`

• `color(purple)(e_(n+1) ~~ (F'''(s))/(3!)(e_n)^3\ text(vs.)\ e_(n+1) = c(e_n)^2)`

• `color(purple)(text(Fixed point converges faster than Newton's method does))`

We say a sequence has order of convergence `q` if `abs(e_(n+1))/abs(e_n)^q` converges to some nonzero number. (Sometimes the nonzero number is referred to as the "asymptotic error constant".) Higher orders of convergence means the sequence converges faster.

Solving Nonlinear Equations: Polynomial Interpolation

Going back to the root-finding problem, we've been able to find just one root. However, it would be nice to be able to use the algorithms to find all the roots. It turns out that we actually can do this with polynomials because of their special properties.

Polynomial Theorems

`p(x) = a_0 + a_1x + a_2x^2 + ... + a_nx^n` (with `a_n != 0`) is a polynomial. The degree of a polynomial is the highest power, which in this case is `n`, denoted as `deg\ p(x) = n`.

The Fundamental Theorem of Algebra says that every nonconstant polynomial has at least one root in the complex field. In fact, an extension of this theorem states something more powerful: if the degree of a polynomial is `n`, then there are exactly `n` roots in the complex field. What we would like is to use our algorithms to find all these (real) roots.

We can divide polynomials just like we can divide numbers. If we divide a polynomial `p(x)` (which has degree of at least `1`) by a linear polynomial `(x-s)`, then the result is a quotient `q(x)` (which is a polynomial) and a remainder `r`.

The Remainder Theorem says that if we divide a polynomial `p(x)` by a linear polynomial `(x-s)`, then `p(s)` is the remainder `r`.

In fact, if `s` is a root, then `p(x) = (x-s)q(x)` (notice there's no remainder.)

A special case of the Remainder Theorem — called the Factor Theorem — combines the two points above: `p(s) = 0 iff p(x) = (x-s)q(x)`. In words, if `s` is a root of `p(x)`, then `p(x)` can be divided evenly by `(x-s)`, which means `p(x) = (x-s)q(x)` (and if `p(x)` can be divided evenly by `(x-s)`, then `p(s) = 0`).

Divide, Rinse, & Repeat (To a Certain Degree)

So we can use the algorithms to find one root of `p(x)`, call it `r_1`, and then divide `p(x)` by `(x-r_1)` to get `p(x) = (x-r_1)q(x)`. Then we can use the algorithms to find the root of `q(x)`, call it `r_2`, and divide `q(x)` by `(x-r_2)` to get `p(x) = (x-r_1)(x-r_2)u(x)`. Wait, it seems like we can just keep on dividing forever. There's actually a limit to the number of times we can divide polynomials, and it involves looking at the degrees of the polynomials.

Recall that when we multiply like terms with exponents, we add the exponents. Like `x^2 * x^3 = x^5`.

If `p(x)` has degree `n`, then `(x-s)q(x)` also needs to have degree `n` (since they're equal to each other). Well, `(x-s)` has degree `1`, which means `q(x)` needs to have degree `n-1` (notice that gives us `1 + (n-1) = n`). Also notice that the degree of `q(x)` is `1` less than the degree of `p(x)`. In a very rough sense, the degree is decreasing by `1` every time we divide.

We can only divide polynomials `p(x)` by linear polynomials `(x-s)` as long as `p(x)` has degree greater than `0`. Since the degree is "decreasing by `1` every time we divide", if `p(x)` has degree `n`, then we can only divide it `n` times. The result is that:

Giving us the `n` roots stated in the extension of the Fundamental Theorem of Algebra. So to find all the roots of `p(x)` we:

1. Find one root of `p(x)`

2. Remove the root by dividing it out (setting `p(x)` equal to `(p(x))/(x-text(root))`)

3. Repeat

Getting the Points

So we've seen that it's possible to find all the roots of polynomials. But what about finding roots for complicated functions that aren't polynomials? It turns out that complicated functions can be replaced or approximated by polynomials.

Let's suppose we had a table of `n+1` data points `(x_i, y_i)` for `i = 0, ..., n` obtained from some complicated function `f(x)`. Then there exists a unique polynomial with degree less than or equal to `n` that passes through (interpolates) the points. The polynomial interpolating the data points is called the interpolation polynomial. Each `x_i` is called a node.

The polynomial obtained using this formula is called the the Lagrange form of the interpolation polynomial.

Let's say we have a polynomial that interpolates the data. However, we later find out that there are more data points we now need to consider. (Maybe we missed some the first time? Nah, let's say someone handed us more data.) Our current interpolation polynomial won't work with these new points, so we have to change our polynomial. How hard is it to change the interpolation polynomial to fit these new data points along with the old data points?

Let's suppose `q(x)` was our original polynomial interpolating the old data. What we want is to get a new polynomial `p(x) = q(x) + r(x)` where `r(x)` is a polynomial interpolating just the new data. So let's see what `r(x)` is.

Since `q(x)` interpolates the data, `q(x_0) = y_0`

Similarly, `r(x_1) = r(x_2) = 0`. So `r(x) = C(x-x_0)(x-x_1)(x-x_2)`. Let's find out what `C` is.

So `p(x) = q(x) + C(x-x_0)(x-x_1)(x-x_2)` interpolates all of the data. Polynomials obtained using this method are called the Newton form of the interpolation polynomial.

Error of Polynomial Interpolation

So given a set of points `(x_i, f(x_i))`, it's possible to find a polynomial (of lowest possible degree) that interpolates those points. However, the polynomial is just an approximation, so there's going to be some error. How well does an interpolation polynomial approximate `f(x)`?

"Interpolation der Runge-Funktion (n=5)" by Cuvwb is licensed under CC BY-SA 4.0

The blue line is the actual function `f(x) = 1/(1+25x^2)` (known as the Runge function), and the red line is the degree `5` interpolation polynomial. Notice that while they have the same `y`-values at `x = -5, -3, -1, 1, 3, 5`, the `y`-values differ in between those `x`-values.

Suppose `f(x) in C^(n+1)[a,b]color(red)(text(*))` for `x in [a,b]`. Suppose `p(x)` is the interpolation polynomial of degree `le n`. Then:

Note that if `x` is far from the data points `x_i` then `abs(x-x_i)` is large, and `prod_(i=0)^n(x-x_i)` is large, making the error large.

If we have a polynomial interpolating the points on the right, then it's not very likely that it will fit well with the point all the way on the left.

Rolle's Theorem: If `f in C^1[a,b]` and `f(a) = f(b)`, then there exists a `c in (a,b)` such that `f'(c) = 0`

Notice that if `a` and `b` are roots of `f`, then there is a root of `f'` between the roots of `f`. This generalizes for higher-order derivatives. Between any two roots of `f'`, there is a root of `f''`. And so on. This idea is used in the proof of the error.

There are several ways to get better error. One way is to connect all the data points with straight lines. Then use only the data points we need instead of all the points. This concept is known as piecewise linear interpolation.

Since we're looking at only `2` data points (`n=1`), the error of piecewise linear interpolation is:

Ex: Consider `f(x) = e^(2x)` in the interval `[-1,1]` with the `n+1` data points evenly spaced and `x_0 = -1`, `x_n = 1`. Let `h` be the stepsize `= x_(i+1) - x_i = (1-(-1))/n = 2/n`. Fix `x in [-1,1]`, so `x in [x_i,x_(i+1)]` for some `i`. Let `p` be the piecewise linear interpolant. Then the error at `x` is:

`f(x)-p(x) = (f''(xi_x))/2(x-x_i)(x-x_(i+1))` for some `xi_x in [x_i,x_(i+1)]`

Putting an upper bound on the error, we get:

`abs(f(x)-p(x)) le (max_(z in [x_i,x_(i+1)])abs(4e^(2z)))/2abs((x-x_i)(x-x_(i+1)))`

So for `x in [-1,1]`:

`abs(f(x)-p(x)) le (max_(z in [-1,1])abs(4e^(2z)))/2abs((x+1)(x-1))`

`= (4e^2)/2abs((x+1)(x-1))`

We can look at the error of `abs((x+1)(x-1))` graphically. `abs((x-x_i)(x-x_(i+1)))` is a degree `2` polynomial, so it is a parabola with roots at `x_i` and `x_(i+1)`.

So `abs((x-x_i)(x-x_(i+1))) le (h/2)(h/2) = h^2/4`, which means:

`(4e^2)/2abs((x+1)(x-1)) le (4e^2)/2*h^2/4`

`= (e^2h^2)/2`

Notice that the error goes to `0` as `n->oo` or `h->0`. In words, we can get better error if we use more data points or decrease the spacing between each data point.

Chebyshev Polynomials

Let's say instead of working with a table of data points, we're given the function `f(x)`. It turns out that there's a way to optimally choose our nodes so the error is small. (We're essentially creating our own table of data points.) Determining which nodes to choose involves something called Chebyshev polynomials.

There are `2` definitions of Chebyshev polynomials. The first one is a recursive definition:

The second definition of a Chebyshev polynomial is:

The second definition tells us a lot about the properties of Chebyshev polynomials.

• `cos(x)` ranges from `-1` to `1`, so `abs(T_n(x)) le 1`

• `T_n(cos((2j-1)/(2n)pi)) = 0` for `j = 1, ..., n`

• — this tells us about the roots of `T_n(x)`

• `T_n(cos((jpi)/n)) = (-1)^j` for `j = 0, ..., n`

• — this tells us about the min and max of `T_n(x)`

Recall that the error of polynomial interpolation is `f^((n+1))(xi_x)/((n+1)!) prod_(i=0)^n(x-x_i)` for some `xi_x in (a,b)`. We can make the error small by making `prod_(i=0)^n(x-x_i) = (x-x_0)...(x-x_n)` small. Well, `(x-x_0)...(x-x_n)` is a monic polynomial, and there's a property of monic polynomials that allows us to put a bound on the error of them.

Notice that `abs((T_n(x))/2^(n-1)) le 1/2^(n-1) => max_(x in [-1,1]) abs((T_n(x))/2^(n-1)) = 1/2^(n-1)`. Combining this with the result above means that

Also notice that `abs(p(x))` (a monic polynomial of degree `n`) is smallest when `p(x) = (T_n(x))/2^(n-1)`. So `(x-x_0)...(x-x_n)` (a monic polynomial of degree `n+1`) is smallest when `(x-x_0)...(x-x_n) = (T_(n+1)(x))/2^n`.

The nodes `x_0`, ..., `x_n` make the left-hand side equal to `0`. This means the nodes `x_0`, ..., `x_n` are the roots of the Chebyshev polynomial. This suggests that the minimum error is achieved when we pick nodes to be the roots of the Chebyshev polynomial, i.e. `x_i = cos((2i+1)/(2n+2)pi)` for `i = 0, ..., n`.

Newton's Form, Revisited

The process for generating the interpolation polynomial of Newton's form involved generating the interpolation polynomial for `1` data point, and then adding more data points.

So we can say that `c_i` is the coefficient of `x_i` for `q_i(x)` of

The problem right now is that `c_i` is pretty hard to calculate (refer back to first mention of Newton's form). It turns out that there's an easier way to find out what the coefficients are without having to generate the interpolation polynomial piece by piece.

Since `c_i` depends on the values of `f(x)` at `x_0, ..., x_i`, there is a notation for `c_i`:

`f[x_0, ..., x_i]` is called a divided difference of order `i`.

`f[x_0, ..., x_i]` is the coefficient of the `i^(th)` degree term. Also, `f[x_i] = f(x_i)`.

From this, we can see that the coefficients of `x_i` are going to lie along the top-most path.

We can even construct the interpolation polynomial at just some of the points. (From here on, the leftmost column contains the `x`-values and the column immediately to the right of it contains the `y`-values.)

The interpolation polynomial for the points `0,1,2` is `p(x) = 4-4(x-0) + 3/2(x-1)(x-0)`

Properties of Divided Differences

Switching the order of the values in the table doesn't change what the divided differences are. Let's say we had a table of values and the interpolation polynomial `p(x) = a_0 + a_1x + ... + a_3x^3`.

Now let's form another table with the order of the values changed and have an interpolation polynomial `q` for this table.

Since `q` interpolates the same points as `p`, `q(x) = a_0 + a_1x + ... + a_3x^3 = p(x)`. This means the divided differences of `q` are the same as the divided differences of `p`. So `f[text(any permutation of)\ x_0, ..., x_n] = f[x_0, ..., x_n]`.

Earlier, we said that the coefficients lie along the top-most path. However, because of this result about switching order, we can go along any path and still get the interpolation polynomial.

Here, the interpolation polynomial is `p(x) = 0-1(x-1)+3/2(x-1)(x-2)+1/6(x-1)(x-2)(x-0)-1/3(x-1)(x-2)(x-0)(x-3)`

This is more easily seen after rearranging the table:

We can also look at the error of polynomial interpolation in terms of divided differences.

Suppose we have an interpolation polynomial `p` interpolating the data points `(x_i, f(x_i))` for `i = 0, ..., n`. Now, let's add a new data point `(t,f(t))` and look at the error at that point. To interpolate that new point, we would do

At `x = t`,

Recall that the error of polynomial interpolation is (at `x = t`)

Comparing the two equations, we see that

Hermite Interpolation

It's possible that we're given more information about our data points than just their values. We could know something about their first derivatives, second derivatives, and so on. Finding the interpolation polynomial for these situations is called Hermite interpolation.

Hermite interpolation will only work if enough information is given. Formally, for `x_0, ..., x_n`, we need to be given (along with their values `f(x_i)`) `f^((j))(x_i)` for `j = 0, ..., k_i-1`. For example, we could do Hermite interpolation for a table that looks like this:

And we can't do Hermite interpolation for a table that looks like this:

For Hermite interpolation, there exists a unique polynomial of degree `le m` that interpolates the data where `m = (text(total amount of information given)) - 1 = sum_(i=0)^n k_i -1`

Error of Hermite Interpolation

Suppose `f in C^(2n+2)[a,b]` and `p` is the interpolation polynomial of degree `le 2n+1` such that `p(x_i) = f(x_i)` and `p'(x_i) = f'(x_i)` for `i = 0, ..., n`. Then

Just like for regular polynomial interpolation we could use piecewise interpolation and Chebyshev polynomials to reduce the error. For piecewise interpolation, instead of piecewise linear interpolation, it would be piecewise cubic interpolation because of the extra information given through the derivatives.

To actually construct the interpolation polynomial, we follow a process similar to constructing the polynomial using the Newton form (this will actually be considered the Newton form for constructing a Hermite polynomial.) We suppose we have an interpolation polynomial `q` interpolating `f` and `f'` at `x_0, ..., x_(n-1)` and add new data `f` and `f'` at `x_n`.

So we have `p(x) = q(x) + r(x)` where `q(x)` interpolates the old data at `x_0, ..., x_(n-1)` and `r(x)` interpolates just the new data at `x_n`. This is achieved by doing:

But this polynomial only gives us the value of `f(x_n)`. We also need the polynomial to give us the value of `f'(x_n)`. So we end up with the final interpolation polynomial:

In fact,

To find the values of `a_0, ..., a_(2n+2)`, we use divided differences, except we repeat each `x_i` twice (because we have two pieces of information on each `x_i`).

However, notice that `f[x_0,x_0] = (f[x_0]-f[x_0])/(x_0-x_0) = 0/0`. To avoid this, we use a property of divided differences we showed earlier:

So `f[x_0,x_0] = (f'(xi))/(1!)` (for some `xi` between `x_0` and `x_0`) `= f'(x_0)`.

Notice that we repeated `0` three times because we had three pieces of information on `0`.

Also notice that `f[0,0,0] = (1-1)/(0-0) = 0/0`. Like we did earlier, we use `f[0,0,0] = (f''(xi))/(2!)` (for some `xi` between `[0,0]`) `= (f''(0))/2 = 4/2 = 2`

Spline Interpolation

Interpolation polynomials are smooth (have continuous derivatives), but they sometimes have bad error. To get better error, we said that using piecewise interpolation would help, but they are not smooth. Splines try to get the best of both: smoothness and better error.

A spline function is a piecewise function of polynomials on subintervals from `i = 0, ..., n`. The subintervals lie between `t_0`, ..., `t_n`, which are called knots. The knots `t_1, ..., t_(n-1)` are called junctions. A spline function of degree `k` means that

• all polynomial pieces have degree `le k`

• all polynomial pieces have a continuous `(k-1)^(st)` derivative

Usually we choose the knots `t_0, ..., t_n` to be the nodes `x_0, ..., x_n`.

A spline of degree `0` consists of `n` constant polynomials. It is also known as a step function.

A spline of degree `1` consists of `n` linear polynomials. It is the classic piecewise linear interpolation polynomial.

Splines of degree `0` usually aren't continuous and splines of degree `1` aren't smooth. So the splines of interest are splines of degree `ge 2`.

Splines of a certain degree need to have a specific form and need to satisfy certain conditions (e.g. continuity, continuity of first derivatives, ...). Because of this, there are going to be a number of equations and unknowns. What we would like to know is if it is possible to set up a spline of degree `k` that interpolates the data and has continuous `k-1` derivatives.

Splines of degree `2` have this form:

From this, we have `3` unknowns per equation. There are a total of `n` equations from `t_0` to `t_n`. So the total number of unknowns is `3n`.

The number of equations comes from the conditions on interpolation and continuity. For each subinterval `[t_i,t_(i+1)]`, we need `S(t_i) = y_i` and `S(t_(i+1)) = y_(i+1)`. That's `2` equations per subinterval, and there are `n` subintervals. So the total number of equations is `2n`.

We also need continuous `1^(st)` derivatives, i.e. `S'_i(t_i) = S'_(i+1)(t_i)` for `i = 1, ..., n-1`. So that's another `n-1` equations.

The total number of unknowns is `3n` and the total number of equations is `2n + n-1 = 3n-1`. There are more unknowns than equations, so it is possible to set up a quadratic spline that satisfies interpolation and continuity. More formally, quadratic splines interpolate and are `C^1`.

Splines of degree `3` have this form:

From this, we have `4` unknowns per equation, and `n` equations total. So the total number of unknowns is `4n`.

For the number of equations, the same reasoning for the degree `2` case applies to the degree `3` case. The only difference is we now need to consider the continuity of `2^(nd)` derivatives, which is another `n-1` equations (`S''_i(t_i) = S''_(i+1)(t_i)` for `i = 1, ..., n-1`).

So there are more unknowns (`4n`) than equations (`4n-2`), which means it is possible to set up a cubic spline. More formally, cubic splines interpolate and are `C^2`.

Solving Nonlinear Equations: Best Approximation

Instead of finding the polynomial that interpolates the data, sometimes we want to find the polynomial that best approximates the data. These polynomials and functions are going to live inside vector spaces.

Suppose `E` is a normed vector space. Let `G sube E` be a finite-dimensional subspace. For `f in E`, we want to find a `g in G` that best approximates `f` i.e., we want to minimize

If we let `E` be the set of continuous functions and `G` be the set of polynomials, then this models the situation we want: finding a polynomial that best approximates a function.

Finding the best approximation is a minimization problem. Which means we need to use calculus to find critical points. But in order to do this, the function (in this case, the norm) needs to be differentiable. This is often not the case. For example, the `1`-norm and `oo`-norm are both defined with absolute values, which is not differentiable. So we restrict our attention to a certain type of norm.

Using inner-products, we can create the norm

This theorem also tells us how to actually construct the best approximation. Assume `G` has a basis `g_1,...,g_n`. According to the theorem, we want `ltf-g,hgt = 0` for all `h in G`. `h` can be written as a linear combination of `g_1,...,g_n`, so it suffices to have `ltf-g,g_igt = 0` for all `i`.

We can start by looking for `g` in the form `sum_(j=1)^nc_jg_j`. Then `ltf-sum_(j=1)^n c_jg_j,g_igt = 0` for all `i`.

This is an `nxxn` linear system of equations that we can solve to find the `c_j`'s. This system is referred to as normal equations.

If we want the best-fitting degree `le n` polynomial, then we change `G` to be the set of degree `le n` polynomials. The normal equations would be an `(n+1)xx(n+1)` linear system. It seems natural to use `{1, x, ..., x^n}` as the basis for `G`, but it turns out that the Gram matrix becomes the Hilbert matrix, which is ill conditioned.

Instead, we should choose a different basis. It would be ideal if the Gram matrix was a diagonal matrix or the identity matrix. To be a diagonal matrix, we need `ltg_i,g_jgt = 0` for all `i != j`. Then `g_1,...,g_n` would be an orthogonal basis. To be the identity matrix, we additionally need `ltg_i,g_jgt = 1` for all `i = j`. Then `g_1,...,g_n` would be an orthonormal basis.

If `g_1, ..., g_n` are orthonormal, then `c_i = ltf,g_igt` and the best approximation `g` has the form `sum_(i=1)^nltf,g_igt g_i`. If `g_1, ..., g_n` are orthogonal, then we could make them orthonormal by using `g_1/norm(g_1), ..., g_n/norm(g_n)`. So `g` would have the form `sum_(i=1)^nltf,g_i/norm(g_i)gt g_i/norm(g_i)`.

We could use Legendre polynomials to form an orthogonal basis (and therefore get a diagonal matrix). They are defined recursively:

where

Note that `p_0, ..., p_n` are all monic polynomials.

At the beginning of the section, we mentioned that we can't use calculus to study best approximation under the `oo`-norm. So if we want to find a best approximation under the `oo`-norm, then we need another approach.

Let `f` be a polynomial of degree `n`. Consider the interval `[-1,1]` and the norm `norm(f)_(oo) = max_(x in [-1,1]) abs(f(x))` for `f in C[-1,1]`. What is the best approximation among `Pi_(n-1)` under the `oo`-norm?

If `g in Pi_(n-1)` is the best approximation, then `g` minimizes `norm(f-g)_(oo)`. Well,

`f(x)-g(x)` is a polynomial of degree `n`. Recall that the monic degree `n` Chebyshev polynomial minimizes `max_(x in [-1,1]) abs(p(x))` among all monic degree `n` polynomials `p`. So we can use Chebyshev polynomials to minimize `max_(x in [-1,1])abs(f(x)-g(x))` if `f(x)-g(x)` is monic.

We can write `f` as:

Then

Since `f` has degree `n` and `g` has degree `n-1`, `a_n = b_n`.

So the monic version of `f(x)-g(x)` is `(f(x)-g(x))/b_n`.

which is minimized when `(f(x)-g(x))/b_n` is the monic degree `n` Chebyshev polynomial. Then we can solve for `g(x) = f(x)-b_n*T_n/2^(n-1)`.