Description

𝑝(𝑎,𝑥) = 𝑎0 + 𝑎1𝑥 + ⋯ + 𝑎𝑛𝑥𝑛.
1. Write a function poly_eval(a, x) to evaluate 𝑝(𝑎,𝑥) in the naive way by just summing the terms as written. This should accept the coefficients 𝑎0,…,𝑎𝑛 as a vector.
2. The Horner method is an alternative evaluation algorithm in which no powers of 𝑥 are (explicitly) calculated.
For example, for a quadratic 𝑝(𝑥) = 𝑎2𝑥2+𝑎1𝑥+𝑎0 we have 𝑝(𝑥) = 𝑎0 + 𝑥(𝑎1 + 𝑎2𝑥)
Generalise this to polynomials of degree 3 and then degree 𝑛
3. Implement this as a function horner(a, x).
4. Test horner to make sure it gives the same result as poly_eval.
5. Benchmark the two methods to evaluate the polynomial 𝑝(𝑎,𝑥) = 1 + 2𝑥 + 3𝑥2 + ⋯ + 10𝑥10. Which algorithm is more efficient and by how much?
6. Count (approximately) the number of operations required by each of the two algorithms. Does this agree with your benchmarks? (Feel free to perform more benchmarks to check this, e.g. using polynomials with random coefficients.)
√
Exercise 2: Calculating√ 𝑦 The Babylonian algorithm for calculating 𝑥 = 𝑦 is the following iteration:
1 𝑦
𝑥𝑛+1 = 2 (𝑥𝑛 + 𝑥 𝑛).
1. Suppose that 𝑥𝑛 converges to some limiting value√ 𝑥∗ as 𝑛 → ∞. Show that then we must have 𝑥∗ = 𝑦.
2. To show that it does, in fact, converge, suppose (without loss of generality)√ that 𝑥0 < 𝑦.
Show that we then have 𝑥𝑛 < 𝑥𝑛+1 < 𝑦𝑛+1 < 𝑦𝑛 for all 𝑛. Hence the 𝑥𝑛 form an increasing sequence that is bounded above, and hence they converge. (This is a theorem that we will just assume for this course.) 3. Write a function babylon that implements this algorithm. It takes a tolerance 𝜖 and iterates until the residual is < 𝜖. It should return the sequence
(𝑥𝑛).
4. Test your function to make sure it gives the correct result.
5. Use rational initial values to convince yourself that calculating with rationals is a bad idea and that floating point is a good compromise.
√
6. Let 𝛿𝑛 ∶= 𝑥𝑛 − 𝑥∗ be the distance of 𝑥𝑛 from the limiting value 𝑦.
7. Plot the (absolute value of) 𝛿𝑛 as a function of 𝑛 on a suitable combination of linear and log scales. How fast does it seem to be converging?
8. Compare this on the same plot to the bisection algorithm from class. Which is better?
9. With 𝛿𝑛 defined as in [6.], show that 𝛿𝑛+1 ≃ 𝛿𝑛2 (≃ means “is approximately equal to”).
Exercise 3: Collision of two discs [In this question we will finally find an actual use case for solving a quadratic equation!]
Suppose we have two discs located at positions x1 and x2, with velocities v1 and v2, respectively. The discs each have radius 𝑟.
3. Write a function collision that takes all the data and calculates the collision time by using the quadratic formula. Make sure you take account of the answer to [2.]
4. Check that your code works by setting up some combinations of discs where you can do the calculation by hand (e.g. both discs moving at a 45-degree angle).
Exercise 4: Evaluating elementary functions
1. Implement the degree-𝑁 Taylor polynomial approximation exp𝑁(𝑥) for exp(𝑥) around 𝑥 = 0. Make sure that you do not explicitly calculate factorials in your code.
2. Make an interactive visualization using the Interact.jl package, showing exp𝑁 and exp as 𝑁 varies.
3. Make another visualization showing the truncation error, i.e. exp(𝑥) − exp𝑁(𝑥). How does it behave?
4. Calculate a bound for the truncation error using the Lagrange remainder. Use Stirling’s approximation to estimate how many terms you need to take for the error to be of a certain size 𝜖.
5. Implement range reduction for the exponential function: instead of blindly applying a Taylor expansion (which will require many terms for large 𝑥), calculate exp(𝑥) by reducing to the calculation of exp(𝑟) for 𝑟 ∈ [−0.5,0.5] using the relation
exp(2𝑥) = exp(𝑥)2.
6. Make an interactive visualization of the Taylor polynomial approximation to log(1+𝑥), showing visually that it fails outside a certain interval. Which interval is that, and why?
Exercise 5: Exactly representing irrationals using symbolic computing
Here we will represent certain real numbers in the computer exactly. Effectively we will use a type of√ symbolic computation, in which we explicitly keep the symbol 2, rather than approximating it by a floating-point number.
(Although this is not really the subject of the course, it’s useful to remember that there is a whole world of symbolic computation that can be applied to certain problems, and that with a bit of work you often do not need an expensive commercial tool to do this!)
√
Consider the subset of the real numbers given by√ 𝑆 ∶= {𝑎+𝑏 2 ∶ 𝑎,𝑏 ∈ ℚ},

i.e. the set of numbers of the form 𝑎 + 𝑏 2 where 𝑎 and 𝑏 are rational.
√
1. Write down formulae for the sum, difference and product of√ 𝑎1 + 𝑏1 2

and 𝑎2 + 𝑏2 2, showing that the results are in 𝑆.
√ √
2. Show that 1/(𝑎+ 𝑏) is also in 𝑆 by supposing that it equals 𝑐 +𝑑 2 and finding explicit equations for 𝑐 and 𝑑. What type of equations are they? Solve them to find explicit values for 𝑐 and 𝑑 in terms of 𝑎 and 𝑏. [This shows that we can also do division (by non-zero elements) and remain within the set, i.e. that 𝑆 is a field, namely an extension field of ℚ.
]
3. Define a type FieldExtension to represent these number pairs, and the corresponding operations. Also define a show method to print them nicely using a √ symbol (typed as sqrt<TAB>).
√ √
4. Find formulae to represent 1/(1 + 2 + 3) as elements of the corresponding set.
Some Julia tips
Package installation
using Pkg
Pkg.add(“BenchmarkTools”)
Load the package in each Julia session with
using BenchmarkTools
Benchmarking
“`julia
@time f(x)
“`
However, if the time taken is too short then this is not accurate. Instead, use the BenchmarkTools.jl package, with the syntax
“`julia using BenchmarkTools
@btime f(1, $x)
“`
Note that any variables you pass in must be given $ signs like this.
Plotting
The Plots.jl package is used as follows after loading it with using Plots:
• plot(x, y): plots the data with given 𝑥 coordinates and 𝑦 coordinates, joining those points with lines. These should be vectors.
• plot(y): specifying only one argument plots the data against the numbers
1,2,…
• plot!: adding ! adds a new plot to an existing one.
• scatter: plots points instead of lines
• Add yscale=:log10 inside the plotting command to use a logarithmic scale on the 𝑦 axis.
Note that there are small delays when first loading the package and for the first plot. Later plots will be quick.
Interactivity
The Interact.jl package enables us to generate simple interactive visualizations using a slider (and some other “widgets”).
To generate a single slider, wrap a for loop in the @manipulate macro, e.g.
@manipulate for i in 1:10
i^2 end
To manipulate a plot, put a plot as the result at the end of the for loop:
@manipulate for i in 1:10
plot(-5:0.01:5, x -> sin(i * x))
end
You can generate multiple sliders by using a joint for loop:
@manipulate for i in 1:10, j in 0.1:0.1:0.9
i + j
end