Numerical Analysis: A Nonlinear Story

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

Notes provided by Professor Li-Tien Cheng (UCSD)

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 x2,x3,cos(x),ex,...

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*. This seems like an easy and obvious idea until we get to numbers like π, which has an infinite number of digits. Does this mean we need an infinite amount of memory to store π? But computers don't have infinite memory, so what happens then?

*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 π 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×10-16 and 1000000000000000 as 1×1015. Then we would only need to store the 45 and the -16 for 0.45×10-16. And the 1 and the 15 for 1×1015.

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×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×10-2 would become 0.214×10-2. As the name implies, the rest of the digits (after the 3rd) are "chopped off". Under rounding, 0.2146×10-2 would become 0.215×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 =|x-fl(x)|

• relative error =|x-fl(x)||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 |x-fl(x)| and |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 dn+1=dn+2=...=9) and making the denominator as small as possible (by letting d1=1 * and d2=d3=...=0).

*d1 cannot be 0 because of the way numbers are written in scientific notation. We wouldn't write something like 0.083×103; it would be 0.83×102.

*0.999...=1

So the relative error — which we can denote by ε — can be at most 101-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*. 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.

* 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(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 ε to be the largest possible relative error. So all of these relative errors are smaller than ε, i.e. |δ1|,|δ2|,|δ3|≤ε

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 |xx-y| and |yx-y| can be very large when x is close to y. When this happens, |δ| 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 a0,b0 such that f(a0),f(b0) have different signs (e.g. f(a0) is positive and f(b0) 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 c0 (the midpoint) and suppose that f(c0) is positive.

Now we have two choices for what our new interval [a1,b1] could be: [a0,c0] or [c0,b0]. The IVT doesn't tell us if there's a root between a0 and c0, but it does guarantee that there is a root between c0 and b0. So the new interval [a1,b1] should be [c0,b0].

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 dn and use that as our approximation. In this picture, there would be the most error if the actual solution was at bn. In general,

We usually pick as our approximation the midpoint of an and bn. 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 an or at bn.)

|exact-approximate|=|r-cn|≤|bn-an|2

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

1. Start with an interval [a0,b0] where f(a0) and f(b0) have different signs

2. Cut the interval in half at c0, which is the midpoint of [a0,b0]

3. Pick a new (smaller) interval such that either f(a0) and f(c0) have different signs or f(c0) and f(b0) 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 limn→∞cn=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 |b1-a1|? Thinking back to the steps of the algorithm, [a1,b1] represented the new interval after splitting [a0,b0] in half. So [a1,b1] is really half of [a0,b0]. 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 ∞, the right-hand side goes to 0.

This proves that what we want to say is true: limn→∞|r-cn|=0⇒limn→∞cn=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 12n+1(b0-a0), 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 x0 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 x0:

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 x0 and x1.

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 |xn+1-r|. Since xn+1=xn-f(xn)f′(xn), we need xn, f(xn), and f′(xn) to get xn+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, xn-f(xn)f′(xn)=xn+1 from the iterative formula for Newton's method, so we have:

Then, we can let en=xn-r be the error between the approximation and the root (so en+1=xn+1-r).