Roots of Polynomials
Lab 1B: Rootfinding
Due: Sunday, January 20th at 11:59PM
Please type in and execute all those commands below, as you read them; and read the documentation for plot. You will need to use plot frequently in this course. This part of the lab is not marked, but it’s important.
Matlab has a roots function that computes roots of polynomials, and it does a reasonable job. We’ll learn how it works in the section on numerical linear algebra after the midterm. For now, here’s an example:
n16 = 0:15; % If we use something twice, name it
c = [(-1).^n16./factorial(2*n16+1)];
cofs = c(end:-1:1); % reverse the order
r2 = roots(cofs); % roots expects highest coefficient first r = sqrt( r2 ); % We want the square roots of those
That gave us (in r2) the roots of
(2k + 1)!λ
gives the roots of xp(x2) which was our n = 15 approximation to sin(x).
p(λ) =
so after taking positive and negative square roots and adding the zero root,
15
∑ (−1)k k
(1)
k=0
x = [-r; 0; r];
Lab 1B: Rootfinding 2
Figure 1: Zeros of zp(z2) where p(λ) is given by Equation (1); zp(z2) approximates sin(z), in the complex plane for z near enough to 0. It’s a little hard to see on this plot, but the real roots are not bad approx- imations to some true roots of sin(z), namely kπ for k = ±3, ±2, ±1, and 0. By k = ±4 it’s not doing quite so well, and forget it for |k| ≥ 5. We remark that Matlab is getting these roots accurately: it’s just that the degree 15 polynomial isn’t a good approximation to sin(z) for |z| > 4π about.
Figure 1 gives a picture of those zeros. To produce this figure, use the following commands:
Make sure you are able to reproduce the above results. [That part is not for credit, though.]
Roots of other nonlinear functions (For credit)
So, polynomials are easy. What about functions that are not polynomials? Your main reference is chapter 4 in the text.
The main technique is Newton’s method: given an initial guess x0, we compute x1, x2, . . . until we stop1 by
plot(real(x), imag(x), ); axis( );
xlabel( );
ylabel( ); axis([-25,25,-25,25])
‘square’
‘Re(x)’
‘Im(x)’
‘ko’
xn+1 = xn − f(xn) f′(xn)
1The tolerance will tell the program when to stop, as discussed in class. Document your choices as to how the tolerance is used.
or
xn+1 = xn − f(xn) (xn − xn−1) f(xn) − f(xn−1)
(Secant) (the secant method requires two initial guesses x0 and x1) and these will iteratively generate hopefully
better xn.
(Newton)
Lab 1B: Rootfinding 3
Task 1
Write two functions: one for Newton’s method and one for secant method. Save the functions as newtonsMethod and secantMethod, respectively. Your function newtonsMethod should at least take the following inputs:
f: The function that you want to solve.
df: The derivative of the function that you want to solve. x0: The initial guess.
tol: The tolerance specifying the accuracy wanted in the solution and your function secantMethod should take:
f: The function that you would to solve. x0: The first initial guess.
x1: The second initial guess.
tol: The tolerance specifying the accuracy wanted in the solution
It also may be a good idea to have an input such as maxIter that can stop the loop within your functions
to limit the amount of time your program runs.
For this lab, your output for both of these functions should be a vector containing the iterations xn. Generally, your output would just be the final answer, since you wouldn’t need to see the values for every iteration, but the reason why we ask for your output to be a vector is to make Task 1 easier.
Remember to have an appropriate function definition and an appropriate header explaining what the function does.
In order for Matlab to know that f and df also accepts input (since the program needs to evaluate f and df at different values of x every iteration), we need to define them as a certain data type called function_handle, which is done using @(). As an example, the function
f(x) = (x + 1)(x − 1/2)
has roots at x = −1 and x = 1/2. The Matlab code
f = @(x) (x + 1).*(x – 1/2);
will define a function in Matlab. You can call this functions as
>> f(3) ans =
10
If you are unable to get your newtonsMethod or secantMethod functions working, obfuscated code has been provided on OWL in the newtonsMethod_solution.p and secantMethod_solution.p files. These functions are guaranteed to work correctly. You may use these functions for the problems that follow without risking losing marks if your newtonsMethod or secantMethod functions don’t work.
To use the functions download the newtonsMethod_solution.p and secantMethod_solution.p files
Lab 1B: Rootfinding 4
to the directory you are working in. You must download the files (right-click, save file as) not just copy and paste their contents.
Sample code:
% Define the function
f = @(x) (x + 1).*(x – 1/2); df = @(x) 2*x + 1/2;
x0 = -1.2;
x1 = -0.9;
tol = 1e-10; maxIter = 20;
% Find a root using the Newton’s method
x_newton = newtonsMethod_solution(f, df, x0, tol, maxIter);
% Find a root using the secant method
x_secant = secantMethod_solution(f, x0, x1, maxIter, tol);
Task 2
We begin with a polynomial example: find the zeros of y3 − 2y − 5.
When y = 1, this polynomial is −6; when y = 2, this polynomial is −1; when y = 3, this polynomial is 16.
So we expect y0 = 2 is a reasonable initial guess.
y 1 = y 0 − y 03 − 2 y 0 − 5
3 y 02 − 2
because if f(y) = y3 − 2y − 5, then f′(y) = 3y2 − 2.
Compute and print all 16 decimal digits of y0, y1, y2, y3, y4, y5, y6, y7.
Why are the last few answers the same? Compare with the answers obtained using “roots”. What happens if you start instead with an initial guess y0 = 1i.
Task 3 : Kepler’s Equation M = E − e sin E
Here M is the so-called “mean anomaly” and E is called the “eccentric anomaly” of an elliptic orbit. The
√a2 − b2
number e is the eccentricity of the orbit (here e = a if the semi-major axis is a and the semi-minor
axis is b; here e is “not necessarily the base of the natural logarithm”. Wikipedia isn’t very friendly about these things, instead flooding the reader with formulae. For background and a derivation of Kepler’s equation, you can see Boaz Katz’ video, if you like (it’s not necessary to do this assignment).
All we need: e is constant for a given orbiting body (e.g. Mercury has e = 0.205630) and the mean anomaly M is very simple: M = 2π(t−t0)+M0 where T is the orbital period. What we want to do is find E, given
T
M and e.
Write three programs to solve Kepler’s equation:
Lab 1B: Rootfinding 5
(a) by Newton’s method: take E0 = M and
Ek+1 = Ek − f(Ek)
f′(Ek)
for k = 0, 1, 2, . . . until convergence
(b) by Halley’s method: again, take E0 = M and
Ek+1 = Ek − f(Ek)
f′(Ek) − 1 f(Ek)f′′(Ek)
2 f′(Ek) (c) by the secant method: take E0 = M and E1 = M +esinM and
Test your programs with Mercury’s
Ek+1 = Ek − f(Ek)(Ek − Ek−1) . f(Ek) − f(Ek−1)
M = 174.796π , e = 0.205630 . 180
Once you have confirmed that your programs for Newton’s method, Halley’s method and secant method is correct, you will now plot the Mercury’s orbit around the sun:
1. Using the equation
M = 2π (t − t0) + M0 , T
where the orbital period T = 87.9691 days, t0 = 0, the mean anomaly at epoch M0 = 174.796π, 180
and eccentricity is 0.205630, create a vector that stores different values of M, calculated using 100 equally-spaced points t, so that you will plot an complete orbit.
2. For each value of M in your vector, use your Newton’s method function from part (a) to solve for E. Store E as a vector as well.
3. To get the coordinates of Mercury as it is orbiting around the Sun, we use
x = a(cos E − e) y = b sin E ,
where a = 0.387098 AU (you will be responsible to solve for b).
4. Plot x and y in Matlabstyling it as black points.
5. On the same plot, using the ’y*’ style, plot where the Sun would be.
Submission requirements
1. All files should be submitted on OWL on the Lab 1B submission. Only .m and .pdf files will be accepted for grading. One person from each group should submit:
• newtonsMethod.m
Lab 1B: Rootfinding 6
• secantMethod.m
• A script file (or files) that is able to reproduce the results you found in the problems • A pdf copy of your report
Your report should be written in a way such that someone with no programming experience can understand what the problem was, how you solved it and what the results were.
2. A printed copy of your report MUST be submitted no later than Friday, Jan. 18th in class.