📢 Updates

Lecture 9: Hermite interpolation available Assignment 3 available Lecture 8 online Lecture 7 online Schedule now available in "Admin" Problem class exercises now online Lecture 6 online Quiz 1 solutions online Lectures 2-3 exercises have solutions! Lecture 1 exercises have solutions!

1  Introduction

Lecture 1: Introduction to Math5485

Overview

1. Syllabus, Canvas, website, and admin;
   
2. Intro to notebooks;
   
3. Overview of Math5485;
   
4. Intro to Numerical Analysis;
   
5. Exercises

Warning

These lecture notes are mainly a record of what we discussed in class. They should not be used as a replacement for attending lectures or reading books. Please let me know if you find any errors or have any questions about the material.

TipHypothesis

You can use the Hypothesis commenting system available on the right side panel (you can login with your U/Google account). This allows you annotate the notes (private) and leave comments or reply to other people’s comments (public).

Note

I encourage you to download the lecture notes and play around with the Julia code.

Download (.ipynb) Ask a Question ChimeIn (opens in a new tab)

1.1 Notebooks

This website is created using Jupyter Notebooks. They allow you to include Markdown, \(\LaTeX\), and Julia (as well as python and other languages).

Tip

Follow the steps in installing-julia.html and A1 to learn about notebooks and Julia.

Notebooks are made up from “cells”.

• A adds a cell above the currently selected cell,
• B adds a cell below,
• M changes the cell from Julia to markdown,
• Y changes the cell from markdown to Julia,
• Double click or enter to see the syntax and edit cells,
• Press ctrl+enter to exit edit mode,

Exercise 1.1

Add a new markdown cell below this one and write “Hello world!”

1.2 Markdown

Text in this cell can be bold, italic, bold italic, quotes

“The only source of knowledge is experience”
– Einstein

Preformatted text,
    useful for
pseudocode

• Bullet points
  • level 2 bullet points
  • etc..
• Another point

You can also use normal html formatting if you are familiar with that…

Markdown cheat sheet: link (opens in a new tab)

1.3 \(\LaTeX\)

\(\LaTeX\) allows you to display math: e.g. \(\pi \approx \frac{22}{7}\) and \[f(x) = f(a) + f'(a) (x - a) + \frac{1}{2} f''(a) (x-a)^2 + R(x).\]

Multi-line equations are also possible:

\[ \begin{align} \int_a^b f(x) \mathrm{d}x % &\approx \int_a^b p_n(x) \mathrm{d}x \nonumber \\ % &= \sum_{j=0}^n w_j f(x_j) \end{align} \tag{1.1}\]

\(\LaTeX\) cheat sheet: link (opens in a new tab)

Draw symbol to get \(\LaTeX\) syntax: link (opens in a new tab)

1.4 Julia

1+1

2

You can use greek symbols in Julia code: type \pi and tab to get the symbol.

π 

π = 3.1415926535897...

You can add comments with a #:

# Define p
p = 22/7;

You can also use @show to print variables (which also prints the calculation):

@show abs(π-p)/π;

abs(π - p) / π = 0.0004024994347707008

You can use println() to output variables (and you can also add text):

println( "The relative error when approximating π by 22/7 is ", abs(π-p)/π)

The relative error when approximating π by 22/7 is 0.0004024994347707008

1.4.1 Arrays, vectors, and matrices

Arrays start at 1:

a = [1 2 3 4 5]

1×5 Matrix{Int64}:
 1  2  3  4  5

a[1]

1

Vectors:

b = [2, 4, 6, 8]

4-element Vector{Int64}:
 2
 4
 6
 8

Vector of objects:

c = [ a , b , "hello world" ]

3-element Vector{Any}:
 [1 2 … 4 5]
 [2, 4, 6, 8]
 "hello world"

c[3]

"hello world"

Matrices

You can either specify the columns:

X = [ [1, 4] [2, 5] ]

2×2 Matrix{Int64}:
 1  2
 4  5

Or the rows:

Y = [ 1 2 ; 4 5 ]

2×2 Matrix{Int64}:
 1  2
 4  5

Linear Algebra

In order to use a certain package, you need to use the using keyword:

using LinearAlgebra

Recall that

X

2×2 Matrix{Int64}:
 1  2
 4  5

det( X )

-3.0

tr( X )

6

(-1/3) * [ 5 -2 ; -4 1 ]

2×2 Matrix{Float64}:
 -1.66667   0.666667
  1.33333  -0.333333

inv( X )

2×2 Matrix{Float64}:
 -1.66667   0.666667
  1.33333  -0.333333

ϵ, v = eigen( X )

Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
 -0.4641016151377544
  6.464101615137754
vectors:
2×2 Matrix{Float64}:
 -0.806898  -0.343724
  0.59069   -0.939071

v1 = v[1:2,1]; # gets the first column of v
norm( ( X - ϵ[1] * Diagonal([1, 1]) ) * v1 )

2.220446049250313e-16

1.4.2 Functions

f = x -> cos( x );
f( π )

-1.0

# vectorisation = entry-wise functuions --> use dot!
f.( [-π, 0, π] )

3-element Vector{Float64}:
 -1.0
  1.0
 -1.0

# you can also specify the vectorisation at the start:
@. f( [-π, 0, π] )

3-element Vector{Float64}:
 -1.0
  1.0
 -1.0

In order to plot functions we need the Plots package and, in order to write \(\LaTeX\) in the titles and labels, we need LaTeXStrings.

using Plots, LaTeXStrings

g(x) = x/exp(x)
plot( g, -5, 1, label="hello" )

X = -5:.01:1;
Y = X .* exp.(X) ;
 
plot( X, Y; 
    label = L"f(w) = w e^w", 
    xlabel = L"w", ylabel = "not latex", 
    title = "L\"\" is a LaTeX string and \"\" a normal string" ) 

X = range(0, 10, length=100)
Y1 = @. sin(X) / X
Y2 = cos.(X)
plot( X, [Y1 Y2], 
        xlabel=L"x", 
        label=[L"\frac{\sin x}{x}" L"\cos x"], 
        linewidth=[4 2], 
        lc=[:green :blue]) 

Figure 1.1: Ploting two functions on one graph

In Figure 1.1, you can see two graphs!

Note

If you are looking at the code, don’t worry about the @fig-1 or the ::: {.callout-note} syntax - this adds the html/css elements so that the webpage looks nicer.

! is like “hold on” in Matlab (if you are familar with that!)

f3(x) = 1 / (1 + x^2);
f4(x) = 1 / (1 + sin(x)^2);
plot(f3, -π, π, 
        lw=3, 
        label = L"f_1(x) = \frac{1}{1+x^2}");
plot!(f4, -π, π, 
        lw = 3, 
        label = L"f_2(x) = \frac{1}{1+\sin^2(x)}", 
        size = (600, 250), legend = :outertopright )

Figure 1.2: Ploting functions

1.5 Math 5485

• Chapter 1: Basics of Numerical Analysis
• Chapter 2: Solving Nonlinear Equations in 1d
• Chapter 3: Polynomial Interpolation
• Chapter 4: Numerical Integration
• Chapter 5: Linear systems of equations

1.6 Chapter 1: Basics of Numerical Analysis

• Absolute and relative error
• Rate of convergence
• Using plots to show rate of convergence

Suppose we approximate \(x\) with \(\widetilde{x}\). Here, \(x\) could be e.g. a root (or fixed point) of a function, the value of an integral, the solution of an ODE or PDE evaluated at a point, …. .

Definition 1.1 (Errors) We define:

\[\begin{align} \text{Actual error:} &\quad x - \widetilde{x} \\ \text{Absolute error:} &\quad |x - \widetilde{x}| \end{align}\]

And when \(x \neq 0\), the relative error is:

\[\begin{align} \text{Relative error:} \quad \frac{|x - \widetilde{x}|}{|x|} \end{align}\]

Example 1.1 If \(x = 1.5\) and \(\widetilde{x} = 1\), then the absolute error is \(0.5\) and the relative error is \(\frac{0.5}{1.5} = \frac13\).

If \(x = 1.5 \times 10^{42}\) and \(\widetilde{x} = 10^{42}\), then the absolute error is \(0.5 \times 10^{42}\) and the relative error is \(\frac13\).

If \(x = 1.5 \times 10^{-42}\) and \(\widetilde{x} = 10^{-42}\), then the absolute error is \(0.5 \times 10^{-42}\) and the relative error is \(\frac13\).

The relative error is a better measure for small or large \(x\):

Example 1.2 Approximating Planck’s constant \(6.62607015 \times 10^{-34} \text{m}^2 \text{kg } \text{s}^{-1}\) with \(0 \text{m}^2 \text{kg } \text{s}^{-1}\) would give an absolute error of \(\approx 6 \times 10^{-34}\) (which seems small) but a relative error of \(1\). In practice \(0\) is a bad approximation to Planck’s constant!

Definition 1.2 (Number of correct significant digits) Suppose \(x \not= 0\) is approximated by \(\widetilde{x}\). Then, we say this approximation is accurate to \(n\) significant digits if

\[\left| \frac{x-\widetilde{x}}{x} \right| \leq 5 \times 10^{-n}.\]

Example 1.3 Approximating \(e := \lim_{n\to\infty} \big( 1 + \frac1n \big)^n \approx 2.7182818.....\) with \(3, 2.7,\) and \(2.718\) gives

\[\begin{align} \left| \frac{e-\widetilde{x}}{e} \right| = 1.04 \times 10^{-1}, \quad 0.673 \times 10^{-2}, \quad 0.104 \times 10^{-3}, \end{align}\]

respectively.

e = exp(1)
abs( (e-3)/e ), abs( (e-2.7)/e ), abs( (e-2.718)/e )

(0.10363832351432703, 0.006725508837105614, 0.00010367889601972718)

Remark 1.1 (Number of correct significant digits). Suppose \(x \not= 0\) is approximated by \(\widetilde{x}\). Some people say that this approximation is accurate to

\[\eta := - \log_{10} \left| \frac{x-\widetilde{x}}{x} \right|\]

significant digits.

Exercise 1.2

Compute \(\eta\) for approximating \(e\) with \(3\), \(2.7\), and \(2.718\).

We have

η( x, x_tilde ) = - log10( abs(x-x_tilde)/abs(x) )
η( exp(1), 3 ), η( exp(1), 2.7 ), η( exp(1), 2.718 )

(0.9844796209116626, 2.172274852323722, 3.984309635847452)

We now consider a sequence of approximations \(x_n\) converging to \(x\). That is, for all \(\varepsilon > 0\), there exists \(N=N(\varepsilon)\) such that \(|x - x_n| < \varepsilon\) for all \(n > N\). In numerical analysis, we are interested in (i) whether our numerical schemes converge, and (ii) the rate of convergence:

Definition 1.3 (Rate of convergence) Suppose \(\varepsilon_n \geq 0\) and \((\varepsilon_n) \to 0\). Then, \(x_n \to x\) with rate \(O(\varepsilon_n)\) if there exist \(C\) and \(N\) such that

\[ |x - x_n| \leq C\varepsilon_n \]

for all \(n > N\). We also write \(x_n = x + O(\varepsilon_n)\).

Exercise 1.3

Show that \(\frac{n + \cos(n)}{n^2} = O(\frac1n)\) and \(\frac{n^2+3}{n^{42}} = O(\frac1{n^{40}})\).

Since \(\cos(n) \leq n\) for all \(n \in \mathbb{N}\), we have

\[\begin{align} \frac{n + \cos n}{n^2} % \leq \frac{2n}{n^2} = \frac{2}{n}. \end{align}\]

That is, \(\frac{n + \cos n}{n^2}\) is bounded above by a constant multiple of \(n^{-1}\) so \(\frac{n + \cos n}{n^2} = O(n^{-1})\).

Similarly, \(3 \leq 3n^2\) so \(\frac{n^2 + 3}{n^{42}} \leq \frac{4}{n^{40}}\).

When plotting graphs, \(O(n^{-k})\) decay (that is, algebraic decay) can be seen as straight lines on log-log plots (\(-k\) is the slope in such graphs):

N = 1:100;
plot( N.^(-1) ;
    xlabel="n", xaxis=:log, yaxis=:log, 
    linestyle = :dot, lw=2, label = L"n^{-1}")
plot!( N.^(-2); 
    linestyle = :dot, lw=2, label = L"n^{-2}")
plot!( N.^(-5); 
    linestyle = :dot, lw=2, label = L"n^{-5}")

On the other hand, \(O(e^{-c N})\) (that is, exponential decay) can be seen as straight lines on a semi-log y plot (and \(-c\) is the gradient):

plot( exp.(-N);
    xlabel="n", yaxis=:log, 
    linestyle = :dot, lw=2, label = L"e^{-n}")
plot!( exp.(-2 .* N ); 
    linestyle = :dot, lw=2, label = L"e^{-2n}")
plot!( exp.(-5 .* N ); 
    linestyle = :dot, lw=2, label = L"e^{-5n}")

TipTO-DO

• Assignment 1: Due Wednesday (Install Julia),
• Exercises (below)
• (Sign-up to office hours if you like)

1.7 Exercises

Exercise 1.4 (Approximating Euler’s constant)

Let \(e\) denote Euler’s constant.

1. \(e_N := \big( 1 + \frac{1}{N} \big)^N\) is an approximation \(e\). What is the error, absolute error, and relative error in approximating \(e\) with \(e_N\) for \(N = 1,2,...,10\)? What are the number of accurate digits in this approximation?
2. What is the rate of convergence of \(e_N\) to \(e\)?

Let’s define \(e\) and the sequence of errors \(|e - e_N|\):

e = exp(1)

N = 1:100;
err = @. abs( e - (1 + 1/N)^N );

The error is equal to the absolute error in this case and is given by

err[1:10]

10-element Vector{Float64}:
 0.7182818284590451
 0.4682818284590451
 0.34791145808867485
 0.2768755784590451
 0.22996182845904567
 0.19665545671693163
 0.17178213141833298
 0.1524973145086972
 0.13710703674584712
 0.12453936835904278

Computing the relative errors give

err[1:10] ./ e

10-element Vector{Float64}:
 0.26424111765711533
 0.17227125736425472
 0.12798947277880338
 0.10185683307753335
 0.08459822894427681
 0.07234549952033957
 0.0631951145094156
 0.05610062684160521
 0.05043886005873465
 0.045815473235769066

Since \(0.0458... = 4.58... \times 10^{-2} \leq 5 \times 10^{-2}\), \(e_{10}\) is correct to \(2\) significant (decimal) digits.

One can also compute the \(\eta = - \log_{10} |\frac{e - e_N}{e}|\):

- log10.( err[1:10]/e )

10-element Vector{Float64}:
 0.5779996023832055
 0.7637871764659977
 0.8928257498997241
 0.992009830990522
 1.072638728778693
 1.1405884803607715
 1.199316494876063
 1.2510322861176486
 1.2972347372416126
 1.338987823174671

Looking at the following plot, we can see that \(| e - e_N | = O(N^{-1})\):

plot( err ;
    xaxis=:log, yaxis=:log, xlabel=L"N", title="absolute error", legend=false, lw=2, m=2 )
plot!( N.^(-1), color=:black, linestyle=:dash, label = L"n^{-1}" )

Exercise 1.5 (Approximating Euler-Mascheroni constant)

Euler-Mascheroni constant, \(\gamma\), is defined as the limit of \(\gamma_n = \sum_{k=1}^n \frac{1}{k} - \log n\) as \(n \to \infty\).

1. What is the absolute error in approximating \(\gamma\) with \(\gamma_{1000}\)?
2. What is the rate of convergence of \(\gamma_n\) to \(\gamma\)?
3. [This was just so that we can remember some calculus and also rigorously prove the rate of convergence that we have seen above using a rigorous error estimate]
   Prove that \(|\gamma - \gamma_n| \leq \frac1n\).

First, we define the Euler-Mascheroni constant:

γ = Base.MathConstants.eulergamma

γ = 0.5772156649015...

… and compute the absolute error in approximating \(\gamma\) with \(\gamma_{1000}\):

Nmax = 1000;

N = 1:Nmax;
γN = zeros( Nmax );
γN[1] = 1;
for n ∈ 2:Nmax
    γN[n] = (1/n) + γN[n-1] + log( (n-1)/n );
end

println( "γ₁₀₀₀ = ", γN[1000] )
println( "Abs. error = " , abs( γ - γN[1000] ) )

γ₁₀₀₀ = 0.5777155815682078
Abs. error = 0.0004999166666749266

Here’s a plot of the absolute errors on a log-log; the gradient matches \(O(N^{-1})\):

plot( γN .- γ , xaxis=:log, yaxis=:log , legend = false, lw = 2, m=2)   
plot!( N.^(-1) , linestyle=:dash)

Proof. First, notice that

\[\begin{align} \gamma_n - \gamma_N % &= \log \frac Nn - \sum_{k=n+1}^N \frac{1}{k}\\ % &= \sum_{k=n+1}^N \left( \log \frac{k}{k-1} - \frac{1}{k} \right). \end{align}\]

Here, \(\frac1k\) is the area under the rectangle with base \([k-1,k]\) and height \(\frac1k\) and \(\log \frac{k}{k-1}\) is the integral of \(\frac1x\) between \(x = k-1\) and \(x = k\):

k=3; x = range(k-1,k, 100);
plot( x, [x.^(-1) ones(100)./k (x .+ 1).^(-1)], 
    ylims=(0,1/(k-1)), 
    label=[L"y=\frac{1}{x}" L"\frac{1}{k}" L"y=\frac{1}{x+1}"] )

In the graph above, \(\frac1k\) is the area under the red line, \(\log \frac{k}{k-1}\) is the area under the blue line. We upper bound the absolute value of the error between these areas, by computing the area between \(y = \frac1x\) and \(y=\frac1{x+1}\):

\[\begin{align} 0 \leq \gamma_n - \gamma_N % &\leq \sum_{k=n+1}^N \left( \log \frac{k}{k-1} - \frac{1}{k} \right) \\ % &\leq \sum_{k=n+1}^N \int_{k-1}^k \left( \frac{1}{x} - \frac{1}{x + 1} \right) \mathrm{d}x % = \int_{n}^\infty \frac{1}{x(x+1)} \mathrm{d}x \\ % &\leq \int_{n}^\infty \frac{1}{x^2} \mathrm{d}x = \frac1n \end{align}\]

Here, we have shown that \((\gamma_n)\) is a Cauchy sequence and thus it converges, and we have the same estimate for \(\gamma_n - \gamma\) by taking the limit in \(N\).

Exercise 1.6 (Approximating π)

Assignment 1, Fall 2025, 30 pts The quantities

\[\begin{align} \pi_N := 4 \sum_{n=0}^N \frac{(-1)^n}{2n+1} \end{align}\]

approximate \(\pi\) by using the Taylor expansion of \(\arctan\) together with the fact that \(\pi = 4 \arctan(1)\).

using Plots, LaTeXStrings

Nmax = 100;
Ninter = 1:Nmax;

# vectors of length Nmax
p = zeros( Nmax ); q = zeros( Nmax );
p[1] = 4 * ( 1 - 1/3 ) ;
q[1] = abs( p[1] - π );
for N in 2:Nmax
    p[N] = p[N-1] + 4 * (-1)^N/(2*N + 1);
    q[N] = abs( p[N] - π );
end

plot(p, legend=false, xlabel="N", title=L"\pi_N")

Based on the following graph, what is the rate of convergence of this particular approximation to \(\pi\)? Why do we use a log-log scale for this graph?

plot( q , xaxis=:log, yaxis=:log, xlabel=L"N", title="absolute error", legend=false)
plot!( Ninter.^(-1) , ls=:dash)

We have \(|\pi - \pi_N| = O(N^{-1})\). The log-log plot makes this behaviour easy to see.

In fact, one can show that

\[\begin{align} |\pi - \pi_N| &= 4 \int_{0}^1 \frac{t^{2(N+1)}}{1 + t^2} \mathrm{d}t\nonumber\\ % &\leq \frac{4}{2N+3}. \end{align}\]