Numerical Analysis Notes

Math 342 - Spring 2026

Week 1 Notes

Mon, Jan 12

We talked about how computers store floating point numbers. Most modern programming languages store floating point numbers using the IEEE 754 standard.

In the IEEE 754 standard, a 64-bit floating point number has the form

$$\mathbf{value} = (-1)^\mathbf{sign} \, \times \, 2^{(\mathbf{exponent} - 1023)} \, \times \, (1 + \mathbf{fraction})$$

where

• $\mathbf{sign}$ is a single bit: 0 for positive, 1 for negative.
• $\mathbf{fraction}$ is a 52-bit fraction (in binary) with $0 \le \mathbf{fraction} < 1$.
• $\mathbf{exponent}$ is the 11-bit binary integer which ranges from 0 to 2047. Only 1 to 2046 are used for regular floating point numbers, $\mathbf{exponent}=0$ is reserved for zero and subnormal numbers, and $\mathbf{exponent}=2047$ is reserved for infinity and NaN (“not a number”).

We talked about how to convert between binary numbers and decimal numbers. This system works very well, but it leads to weird outcomes like 0.1 + 0.1 + 0.1 = 0.30000000000000004.

When you convert $\tfrac{1}{10} = 0.1$ to binary, you get a infinitely repeating binary decimal: $$\tfrac{1}{10} = (0.000110011001100\ldots)_2.$$ So any finite binary representation of $\tfrac{1}{10}$ will have rounding errors.

We did the following examples in class:

1. Convert (10110)₂ to decimal. (https://youtu.be/a2FpnU9Mm3E)

2. Convert 35 to binary. (https://youtu.be/gGiEu7QTi68)

3. What is the 64-bit string that represents the number 35 in the IEEE standard?

4. What are the largest and smallest 64-bit floating point numbers that can be stored?

5. In Python, compute 2.0**1024 and 2**1024. Why do you get different results?

6. In Python, compare 2.0**1024 with 2.0**(-1024) and 2.0**(-1070). What do you notice?

Wed, Jan 14

Today we talked about significant digits. Here is a quick video on how these work.

Next we defined absolute and relative error:

The base-10 logarithm of the relative error is approximately the number of significant digits, so you can think of significant digits as a measure of relative error.

Intuitively, addition & subtraction “play nice” with absolute error while multiplication and division “play nice” with relative error. This can lead to problems:

1. Catastrophic cancellation. When you subtract two numbers of roughly the same size, the relative error can get much worse. For example, both 53.76 and 53.74 have 4 significant digits, but $$53.76 - 53.74 = 0.02$$ only has 1 significant digit.

2. Useless precision. If you add two numbers with very different magnitudes, then having a very low relative error in the smaller one will not be useful.

We did these examples in class.

1. $\pi = 3.141592...$. What is the absolute and relative error if your round $\pi$ to $3.14$?

You can sometimes re-write algorithms on a computer to avoid issues with floating point numbers such as overflow/underflow and catastrophic cancellation.

2. Consider the function $f(x) = \dfrac{1 - \cos x}{\sin x}$. Use Python to compute $f(10^{-7})$.

3. The exact answer to previous question is $0.00000005 = 5 \times 10^{-8}$ (accurate to 22 decimal places). Use this to find the relative error in your previous calculation.

4. A better way to compute $f(x)$ is to use a trick to avoid the catastrophic cancellation: $$f(x) = \dfrac{1-\cos x}{\sin x} = \dfrac{1 - \cos x}{\sin x} \cdot \left( \frac{1+ \cos x}{1+\cos x} \right) = \dfrac{\sin x}{1 + \cos x}.$$ Use this new formula to compute $f(10^{-7})$. What is the relative error now?

Stirling’s formula is a famous approximation for the factorial function. $$n! \approx \sqrt{2 \pi n} \frac{n^n}{e^n}.$$ We approximated Stirling’s formula using the following code.

import math
n = 100
print(float(math.factorial(n)))
f = lambda n: math.sqrt(2 * math.pi * n) * n ** n / math.exp(n)
print(f(n))

Our formula worked well until $n=143$, then we got an overflow error. The problem was that $n^n$ got too big to convert to a floating point number. But you can prevent the overflow error by adjusting the formula slightly to.

n = 143
f = lambda n: math.sqrt(2 * math.pi * n) * (n / math.e) ** n
print(f(n))

Fri, Jan 26

Today we reviewed Taylor series. We recalled the following important Maclaurin series (which are Taylor series with center $c = 0$):

• $e^x = \sum_{n=0}^\infty \dfrac{x^n}{n!}$
• $\sin(x) = \sum_{n=0}^\infty \dfrac{(-1)^n \, x^{2n+1}}{(2n+1)!}$
• $\cos(x) = \sum_{n=0}^\infty \dfrac{(-1)^n \, x^{2n}}{(2n)!}$
• $\dfrac{1}{1-x} = \sum_{n = 0}^\infty x^n$
• $\ln(1+x) = \sum_{n = 0}^\infty \dfrac{(-1)^n \, x^{n+1}}{n+1}$

We graphed Maclaurin polynomials for $\cos x$ and $\dfrac{1}{1-x}$ on Desmos to see how they converge with different radii of convergence.

We also use the Maclaurin series for $\dfrac{\sin x}{x}$ to approximate the integral

$$\displaystyle \int_{-\pi}^\pi \frac{\sin x}{x} \, dx.$$

Then we did the following workshop in class.

• Workshop: Taylor series

Week 2 Notes

Wed, Jan 21

Today we reviewed some theorems that we will need throughout the course. The first is probably the most important theorem in numerical analysis since it lets us estimate error when using Taylor series approximations.

A special case of Taylor’s theorem is when $n = 0$. Then you get the Mean Value Theorem (MVT):

We did this example:

1. Use Taylor’s theorem to estimate the error in using the 0th and 2nd degree Maclaurin polynomials to estimate $\cos(0.03)$ and $\cos(0.6)$.

Then we started this workshop

• Workshop: Error bounds

Fri, Jan 23

Last time we started this workshop about using Taylor’s remainder formula and the triangle inequality to find upper bounds for functions. Today we revisited that workshop, but first we talked about the following.

This error formula gives a way to estimate the worst case (absolute) error when you use a Taylor polynomial approximation.

1. The 3rd degree Taylor polynomial for $\cos x$ centered at $c = \pi$ is $P_3(x) = -1 + \tfrac{1}{2}(x-\pi)^2$ (the coefficient on the 3rd degree term is zero). What is the worst case absolute error using this polynomial to estimate $\cos 3$?

After that we talked about the triangle inequality.

We talked about how you can use the triangle inequality to find upper bounds for functions. We also talked about tight upper bounds versus upper bounds that are not tight. We did this example.

2. Use the triangle inequality to find an upper bound for $|x^2 + 3x \cos x|$.

Week 3 Notes

Wed, Jan 28

We talked about how to find the roots of a function. Recall that a root (AKA a zero) of a function $f(x)$ is an $x$-value where the function hits the $x$-axis. We introduced an algorithm called the bisection method for finding roots of a continuous function that changes sign from positive to negative or negative to positive on an interval $[a, b]$. We wrote the following code to implement this algorithm.

def bisection(f, a, b, N):
    """
    Applies the bisection method recursively up to N times to estimate a root 
    of a continuous function f on an interval [a,b]. 
    """
    m = (a + b) / 2
    if N == 0 or f(m) == 0:
        return m
    if (f(a) > 0 and f(m) > 0) or (f(a) < 0 and f(m) < 0):
        return bisection(f, m, b, N - 1)
    else:
        return bisection(f, a, m, N - 1)

We tested this algorithm on the function $f(x) = \tan x - 1$ which has a root at $\tfrac{\pi}{4}$.

One feature of the bisection method is that we can easily find the worst case absolute error in our approximation of a root. That is because every time we repeat the algorithm and cut the interval in half, the error reduces by a factor of 2, so that $$\text{Absolute error} \le \frac{(b-a)}{2^{N+1}}.$$ We saw that it takes about 10 iterations to increase the accuracy by 3 decimal places (because $2^{10} \approx 10^3$).

Fri, Jan 30

Today we covered Newton’s method. This is probably the most important method for finding roots of differentiable functions. The formula is $$x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)}.$$ This formula comes from the idea which is to start with a guess $x_0$ for a root and then repeatedly improve your guess by following the tangent line at $x_n$ until it hits the $x$-axis.

1. Use Newton’s method to find roots of $\tan x - 1$.

2. How can you use Newton’s method to find $e$? Hint: use $f(x) = \ln x -1$.

Proof. Start with the first degree Taylor polynomial (centered at $x_n$) for $f(r)$ including the remainder term and the Newton’s method iteration formula:

$$f(r) = f(x_n) + f'(x_n)(r-x_n) + \frac{1}{2} f''(z)(r-x_n)^2 = 0,$$ and $$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \Rightarrow f'(x_n)(x_{n+1} - x_n) + f(x_n)=0.$$

Subtract these two formulas to get a formula that relates $(r-x_{n+1})$ with $(r-x_n)$.

$$f'(x_n)(r-x_{n+1}) + \frac{1}{2} f''(z)(r-x_n)^2 = 0.$$

Use this to get an upper bound on $|r-x_{n+1}|$. □

This corollary explains why, if you start with a good guess in Newton’s method, the number of correct decimal places tends to double with each iteration!

Week 4 Notes

Mon, Feb 2

The secant method is a variation of Newton’s method that uses secant lines instead of tangent lines. The advantage of the secant method is that it doesn’t require calculating a derivative. The disadvantage is that it is a little slower to converge than Newton’s method, but it is still much faster than the bisection method. Here is the formula:

$$x_{n+1} = x_n - \frac{f(x_n) \, (x_n - x_{n-1})}{f(x_n) - f(x_{n-1})}.$$

1. Solve the equation $2^x = 10$ using the secant method. What would make good initial guesses $x_0$ and $x_1$?

    # This code will do one step of the secant method.
    f = lambda x: 2 ** x - 10
    a, b = 3, 3.5
    a, b = b, b - f(b) * (b - a) / (f(b) - f(a)); print(b)

Convergence of order $\alpha = 1$ is called linear convergence and convergence of order $\alpha = 2$ is called quadratic convergence. Newton’s method converges quadratically. It turns out that the secant method converges with order $\alpha = \frac{1+\sqrt{5}}{2} \approx 1.618$ which is the golden ratio!

• Workshop: Root finding methods

Wed, Feb 4

Newton’s method is a special case of a method known as fixed point iteration. A fixed point of a function $f(x)$ is a number $p$ such that $f(p) = p$. A real-valued function $f(x)$ has a fixed point on an interval $[a, b]$ if and only if it intersects the graph $y = x$ on that interval.

1. Show that $\cos x$ has a fixed point in $[0,\tfrac{\pi}{2}]$.

2. Explain why $f(x) = e^x$ has no fixed points.

3. Does $\sin x$ have any fixed points?

A fixed point $p$ is attracting if the recursive sequence defined by $x_{n+1} = f(x_n)$ converges to $p$ for all $x_0$ sufficiently close to $p$. It is repelling if points close to $p$ get pushed away from $p$ when you apply the function $f$. You can draw a picture of these fixed point iterates by drawing a cobweb diagram.

To make a cobweb diagram, repeat these steps:

• Move vertically from your current x-value until you hit the graph $y = f(x)$.
• Now move horizontally to the graph $y = x$.

4. Show that the fixed point of $\cos x$ is attracting by repeatedly iterating.

5. Show that $g(x) = 1 - 2x -x^5$ has a fixed point, but it is not attracting.

Proof. Write down the first degree Taylor approximation for $f$ centered at $p$. Use it with $x = x_n$, $f(x_n) = x_{n+1}$ to show that if $M$ is an upper bound for $f''(z)$ (near $p$), then $$|x_{n+1} - p| \le |x_n - p| \left(|f'(p)| + \frac{M}{2}|x_n-p| \right).$$ Then as long as $x_n$ is close enough to $p$, the terms inside the parentheses are less than 1 which means that $|x_{n+1} - p| < |x_n - p|$, i.e, $x_{n+1}$ is closer to $p$ than $x_n$ for every $n$. □

5. For any function $f \in C^2[a,b]$ with a root $r$ in $(a,b)$, let $g(x) = x - \dfrac{f(x)}{f'(x)}$. Show that $r$ is a fixed point of $g$ and show that $g'(r) = 0$.

Fri, Feb 6

In any root finding method, we have two ways of measuring the error when we compare an approximation $x_n$ with the actual root. The forward error is the distance between the root $r$ and the approximation $x_n$ on the $x$-axis. Since we usually don’t know the true value of $r$, it is hard to estimate the forward error. The backward error is just $|f(x_n)|$, which is usually easy to calculate.

1. If $x_n$ is really close to $r$, then how could you use $f'(r)$ to estimate the forward error if you know the backward error?

We used the idea of forward and backward error to automate Newton’s method. The idea is to include an optional tolerance argument. We stop iterating when the backward error is smaller than the tolerance.

def newton(f, Df, x, tolerance = 10 ** (-15)):
    """
    Applies Newton's method to a function f with derivative Df and initial guess x.  
    Stops when |f(x)| < tol. 
    """
    step = 0
    while not abs(f(x)) < tolerance:
        x = x - f(x) / Df(x)
        step = step + 1
    print("Converged in", step, "steps.")
    return x

We finished with a cool fact about Newton’s method. It also works for to find complex number roots if you use complex numbers. We did a quick review of complex numbers.

We used the cmath library in Python to do the following experiments with Newton’s method.

2. Cube roots of unity. Use Newton’s method to find all 3 roots of $x^3 - 1$.

3. Euler’s identity. Use Newton’s method to solve $e^x + 1 = 0$.

We finished by talking about Newton fractals. When you use Newton’s method on a function with more than one root in the complex plane, the set of points that converge to each root (called the basins of attraction) form fractal patterns.

Week 5 Notes

Mon, Feb 9

Today we talked about how to solve systems of nonlinear equations with Newton’s method. As a motivating example, we talked about how this could be used in navigation, for example in the LORAN system.

To solve a vector-valued system $\mathbf{F}(\mathbf{x}) = \mathbf{0}$, you can iterate the formula $$\mathbf{x}_{n+1} = \mathbf{x}_n - \mathbf{J}(\mathbf{x}_n)^{-1} \mathbf{F}(\mathbf{x_n})$$ where $\mathbf{J}(\mathbf{x})$ is the Jacobian matrix $$\mathbf{J}(\mathbf{x}) = \begin{bmatrix} \frac{\partial F_1}{\partial x_1} & \ldots & \frac{\partial F_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial F_n}{\partial x_1} & \ldots & \frac{\partial F_n}{\partial x_n} \end{bmatrix}.$$

To program this formula in Python, we’ll need to load the numpy library. We wrote the following code to solve this system:

$$\begin{align*} xy + y^2 &= 1 \\ x^2 - y^2 &= 1 \end{align*}$$

import numpy as np

F = lambda x: np.array([
    x[0] * x[1] + x[1]**2 - 1, 
    x[0]**2 - x[1]**2 - 1
])
J = lambda x: np.array([
    [   x[1], x[0] + 2*x[1] ], 
    [ 2*x[0],       -2*x[1] ]
])

x = np.array([1,1])
for i in range(10): 
    x = x - np.linalg.inv(J(x)) @ F(x) # Use @ not * for matrix multiplication.
    print(x)

• Workshop: Nonlinear systems

Wed, Feb 11

Today we talked about systems of linear equations and linear algebra.

1. Suppose you have a jar full of pennies, nickles, dimes, and quarters. There are 80 coins in the jar, and the total value of the coins is $10.00. If there are twice as many dimes as quarters, then how many of each type of coin are in the jar?

We can represent this question as a system of linear equations. $$p+n+d+q = 80$$ $$p+5n+10d+25q = 1000$$ $$d = 2q$$ where $p,n,d,q$ are the numbers of pennies, nickles, dimes, and quarters respectively. It is convenient to use matrices to simplify these equations: $$\begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & 5 & 10 & 25 \\ 0 & 0 & 1 & -2 \end{pmatrix} \, \begin{pmatrix} p \\ n \\ d \\ q \end{pmatrix} = \begin{pmatrix} 80 \\ 1000 \\ 0 \end{pmatrix}.$$ Here we have a matrix equation of the form $A\mathbf{x} = \mathbf{b}$ where $A \in \mathbb{R}^{3 \times 4}$, $\mathbf{x} \in \mathbb{R}^4$ is the unknown vector, and $\mathbf{b} \in \mathbb{R}^3$. Then you can solve the problem by row-reducing the augmented matrix

$$\left( \begin{array}{cccc|c} 1 & 1 & 1 & 1 & 80 \\ 1 & 5 & 10 & 25 & 1000 \\ 0 & 0 & 1 & -2 & 0\end{array}\right)$$

which can be put into echelon form

$$\left( \begin{array}{cccc|c} 1 & 1 & 1 & 1 & 80 \\ 0 & 4 & 9 & 24 & 920 \\ 0 & 0 & 1 & -2 & 0\end{array}\right).$$

The variables $p, n$, and $d$ are pivot variables, and the last variable $q$ is a free variable. Once you pick values for the free variable(s), you can solve for the pivot variables one at a time using back substitution. We did this using Python:

q = 20 
d = 2*q
n = (920 - 9*d - 24*q) / 4
p = 80 - n - d - q
print(p, n, q)

# Only q = 20 makes sense, otherwise either p or n will be negative.

After that we reviewed some more concepts and terminology from linear algebra. The most important thing to understand is that you can think of matrix-vector multiplication as a linear combination of the columns of the matrix:

$$A \mathbf{x} = \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & 5 & 10 & 25 \\ 0 & 0 & 1 & -2 \end{pmatrix} \, \begin{pmatrix} p \\ n \\ d \\ q \end{pmatrix} = \underbrace{\begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix}p + \begin{pmatrix} 1 \\ 5 \\ 0 \end{pmatrix}n + \begin{pmatrix} 1 \\ 10 \\ 1 \end{pmatrix}d + \begin{pmatrix} 1 \\ 25 \\ -2 \end{pmatrix}q}_{A \mathbf{x} \text{ is a linear combination of the columns of } A}.$$

For any matrix $A \in \mathbb{R}^{m \times n}$:

• The column space of $A$ is the set of all possible linear combinations of the columns. It is also the range of $A$ as a linear transformation from $\mathbb{R}^n$ into $\mathbb{R}^m$.

• The rank of $A$ is the dimension of the column space. It is also the number of pivots, since the pivot columns are a basis for the column space.

• The null space of $A$ is the set $\{\mathbf{x} \in \mathbb{R}^n \, : \, A \mathbf{x} = \mathbf{0}\}$.

• The nullity of $A$ is the number of free variables which is the same as the dimension of the null space of $A$.

A matrix equation $A\mathbf{x} = \mathbf{b}$ has a solution if and only if $\mathbf{b}$ is in the column space of $A$. If $\mathbf{b}$ is in the column space, then there will be either one unique solution if there are no free variables (i.e., the nullity of $A$ is zero) or there will be infinitely many solutions if there are free variables.

If $A \in \mathbb{R}^{n \times n}$ (i.e., $A$ is a square matrix) and the rank of $A$ is $n$, then $A$ is invertible which means that there is a matrix $A^{-1}$ such that $A A^{-1} = A^{-1} A = I$ where $I$ is the identity matrix which is $$I = \begin{pmatrix} 1 & 0 & \ldots & 0 \\ 0 & 1 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & 1 \end{pmatrix}.$$
You can use row-reduction to find the inverse of an invertible matrix by row reducing the augmented matrix $\left[ \begin{array}{c|c} A & I \end{array} \right]$ until you get $\left[ \begin{array}{c|c} I & A^{-1} \end{array} \right]$. We did this example in class:

2. Find the inverse of $A = \begin{pmatrix} 1 & 2 & 0 \\ 0 & 1 & 3 \\ 0 & 0 & 1 \end{pmatrix}$.

Here is another example that we did not do in class:

3. Use row-reduction to find the inverse of $A = \begin{pmatrix} 1 & 3 \\ 2 & 5 \end{pmatrix}$. (https://youtu.be/cJg2AuSFdjw)

Fri, Feb 13

Today we talked about LU decomposition. We defined the LU decomposition as follows. The LU decomposition of a matrix $A \in \mathbb{R}^{m \times n}$ is a pair of matrices $L \in \mathbb{R}^{m \times m}$ and $U \in \mathbb{R}^{m \times n}$ such that $A = LU$ and $U$ is in echelon form and $L$ is a lower triangular matrix with 1’s on the main diagonal, 0’s above the main diagonal, and entries $L_{ij}$ in row $i$, column $j$ that are equal to the multiple of row $i$ that you subtracted from row $j$ as you row reduced $A$ to $U$.

1. Compute the LU decomposition of $A = \begin{pmatrix} 1 & 1 & 1 & 1 \\ 2 & 2 & 5 & 3 \\ -1 & -1 & 14 & 4 \end{pmatrix}$.

2. Use the LU decomposition to solve $Ax = \begin{pmatrix} 2 \\ 6 \\ 8 \end{pmatrix}$.

We also did this workshop.

• Workshop: LU decomposition

We finished with one more example.

3. For what real numbers $a$ and $b$ does the matrix $\begin{pmatrix} 1 & 0 & 1 \\ a & a & a \\ b & b & a \end{pmatrix}$ have an LU decomposition? (https://youtu.be/-eA2D_rIcNA)

Here’s another LU decomposition example if you want more practice.

4. Decompose $A = \begin{pmatrix} 2 & 4 & 3 & 5 \\ -4 & -7 & -5 & -8 \\ 6 & 8 & 2 & 9 \\ 4 & 9 & -2 & 14 \end{pmatrix}$. (https://youtu.be/BFYFkn-eOQk)

Week 6 Notes

Week 7 Notes

Week 8 Notes

Week 9 Notes

Week 10 Notes

Week 11 Notes

Week 12 Notes

Week 13 Notes

Week 14 Notes