Description

√
𝑓(𝛼,𝑥) = 𝑥3 − 𝛼𝑥 + 2
by using fixed-point iterations.
1. Define a Julia function for 𝑓 and draw it; make an interactive visualization as 𝛼 changes. [Make sure to fix the plot limits so as to be able to see what’s going on.]
2. When 𝛼 varies, the number of real roots can change. For which approximate value of 𝛼 does the number of real roots change? How many real roots are there?
3. When 𝛼 = 2.5, approximately where are the roots? From now on fix 𝛼 to this value.
The simplest thing we can try is 𝑔(𝑥) = 𝑥 + 𝑓(𝛼,𝑥), since then a fixed point of 𝑔 is a root of 𝑓.
4. Plot the function 𝑔(𝑥) and the function 𝑦 = 𝑥. Taking into account the results stated in the slides, which root do you expect to be able to calculate by iterating 𝑔? Fix an initial condition and show that it does converge there.
5. What is the rate of convergence to that root?
6. What happens with the initial condition 𝑥 = 1.1? Why?
7. Draw a cobweb diagram to illustrate this behaviour.
8. By using algebraic transformations, find fixed-point iterations 𝑔 that converge to the other two roots.
There are two alternative approaches you can take here. The first is to find other iteration schemes, 𝑥 = ℎ(𝑥) by algebraically rearranging 𝑓(𝛼,𝑥) = 0 to isolate an 𝑥 on one side. The second is to introduce the generalized function 𝑔(𝑐,𝑥) ∶= 𝑥+𝑐𝑓(𝑥). Make an interactive plot of 𝑔(𝑐,𝑥) and 𝑦 = 𝑥 as 𝑐 varies. What do you notice about the slope of 𝑔(𝑐,𝑥) at the fixed points as 𝑐 changes? Can you use this to change the stabilities of the fixed points in the iteration scheme?
9. (Extra credit) Call 𝛼𝑐 the critical value of 𝛼 at which the number of roots changes. [This is called a bifurcation point.] Find 𝛼𝑐 numerically.
Exercise 2: Defining a type for Polynomials Although we have already used polynomials, we haven’t had a way to say “this object is a polynomial”. For this we need to define a new type!
1. Define a type Polynomial to represent a polynomial. It should have fields degree and coefficients.
2. Write a constructor function with the same name (Polynomial) that accepts a vector of coefficients and builds a Polynomial object whose degree it automatically calculates.
3. Write a show method to display the polynomial nicely.
4. Write a function to evaluate the polynomial at a point x, as follows:
function (p::Polynomial)(x)
# fill in end
Then if you have an object p of type polynomial, writing e.g. p(3) will evaluate that polynomial at the value 3.
5. Write a function derivative that takes a polynomial and returns a new polynomial that is its derivative.
6. Julia contains a module Test (in the standard library – no installation required) for testing that code is correct.
Write a few tests of the functionality you have defined using tests of the form
@test a == b
E.g. to test the sum of two Polynomials you can write
@test Polynomial([1, 2]) + Polynomial([3, 4]) == Polynomial([4, 6]) When you run these tests, you should see the message Test passed.
To create the tests, do the calculations by hand.
Exercise 3: Newton method
1. Write a function newton to implement the Newton method. It should take a function 𝑓, its derivative 𝑓′ and an initial guess 𝑥0. It should terminate when the residual is less than some tolerance, or when the number of iterations is too large.
2. Write a specialised method of newton for polynomials that accepts an object of type Polynomial, written p::Polynomial and calculates its derivative automatically.
3. Taking the same polynomial as in exercise 1, confirm that the Newton method has quadratic convergence by verifying numerically the defining limit.
4. Can you find different roots by taking different initial conditions 𝑥0? How will you know if/when you have found all of the roots like this?
The Newton method works fantastically well when the starting point is close enough to a root, but can also behave very badly, as follows.
5. How many roots does the function 𝑓(𝑧) = 𝑧3 − 1 have in the complex plane ℂ? Where are they?
6. Calculate those roots using the Newton method with complex starting points 𝑎 + 𝑖𝑏 forming a square grid in the complex plane. Store the imaginary part of the resulting root (or final value, if no root is reached) in a matrix.
7. Plot the matrix with the heatmap function and plot the true roots as points. What kind of object are you seeing? What does this imply for the Newton method? What happens close to the roots?
Exercise 4: Aberth method Newton’s method only finds one root at a time. The Aberth method calculates them all at once! Given estimates (𝑧𝑘) for all roots of a polynomial it calculates a new estimate via 𝑧𝑘 = 𝑧𝑘 − 𝑤𝑘, where
𝑘
𝑤𝑘 =.
′ 𝑘 𝑗≠𝑘 𝑧𝑘−𝑧𝑗 In order to guarantee convergence the initial guess must be chosen in a special way. Cauchy’s bound tells us that all the roots of a polynomial 𝑝(𝑥) = 𝑎0 + 𝑎1𝑥 + ⋯ + 𝑎𝑛𝑥𝑛 must have absolute value less than
𝑎𝑛−1∣,∣𝑎𝑛−2∣,…,∣𝑎0∣)
𝑈 = 1 + max (∣
𝑎𝑛 𝑎𝑛 𝑎𝑛
and a lower bound of
𝐿 =
The starting points should be chosen so that their absolute value lies in between these bounds.
1. Implement a function poly_bounds that takes in a Polynomial and returns 𝐿 and 𝑈.
2. Use poly_bounds to write a function initialize_points that finds a suitable starting group of points. Do this by choosing 𝜌s sampled from a uniform distribution over [𝐿,𝑈] and 𝜃s sampled from a uniform distribution over [0,2𝜋]. Then calculate 𝑧 = 𝜌𝑒𝑖𝜃.
[Hint: The rand() function generates uniform random numbers between
0 and 1. Use rand to write a function uniform(a, b) to generate uniform numbers between 𝑎 and 𝑏.]
3. Implement the Aberth method for a polynomial 𝑝 which takes in values from initialize_points and terminates either after a certain number of iterations or when a certain tolerance is met. Test your function for the polynomials we have considered so far.
4. Make an interactive visualization of the progress over time on some random polynomials of degree 10 or 20.
5. Find numerically the order of convergence of the method.
6. Try polynomials with multiple roots, such as (𝑥2 + 1)3. What happens?