Download Report

Math 365
Grady Wright
Homework 6
1. (Newton’s method) The almost universally used algorithm to compute
the recursion
1
a
xn+1 =
xn +
,
2
xn
Due Nov 21, 2014
√
a, where a > 0, is
(1)
easily obtained by means of Newton’s method for the function f (x) = x2 − a. One potential
problem with this method is that it requires a floating point division, which not all computer
processors support, or which may too expensive for a particular application. For these reasons, it is advantageous to devise a method for computing the square root that only uses
addition, subtraction, multiplication, and division by 2 (which can be easily done by shifting
the binary representation one bit to the right).
The trick for doing this is to use Newton’s
√
method to compute √1a , and then obtain a by multipling by a. Write down your recursion formula for computing √1a in a manner similar to (1). This formula should only involve
addition/subtraction, multiplication and division by 2.
√
Try you algorithm on the problem of computing 5. As an initial guess use x0 = 0.5. Report
the values of ax0 , ax1 ,. . . , ax5 in a nice table√and verify that your algorithm is working by
comparing these numbers to the true value of 5.
To see where this sort of software assisted acceleration is used in gaming, see the course
webpage for a link to the article : Origin of Quake3’s Fast InvSqrt().
2. (Julia sets and fractals) You may have seen the beautiful pictures that are produced by
fractals. One type of fractal, the Julia sets, are computed using Newton’s Method.
The equation f (z) = z 3 + 1 has three roots, one real root, and two complex conjugate roots.
They are given by
π
π
π
π
iπ
iπ
r1 = −1, r2 = exp
= cos
+i sin
, r3 = exp −
= cos
−i sin
,
3
3
3
3
3
3
√
where i = −1 and exp(iπ) = −1. It turns out, we can apply Newton’s method to the complex
valued function f (z) = 0 to obtain one of these roots.
If we start with an arbitrary complex number z = x + iy as a starting value for Newton’s
Method, the iterates that are produced will either diverge, or they will converge to one of the
three roots of f (z). If we color each starting value in the complex plane according to which
root it converges to, we get a surprising fractal pattern (see Figure 1). The boundary between
regions of different colors form the Julia set (named after the French mathematician Gaston
Julia).
(a) Download the code julia.m from the course website, and add the Newton step where
indicated.
(b) Convert the code julia.m to a function that you can call with several different sets of
axis limits. Your function should look like
function julia(xlimits,ylimits)
Figure 1: Julia set produced by finding the roots of f (z) = z 3 + 1. The black dots are the three
roots in the complex plane, and the color regions correspond to initial starting values for Newton’s
Method which converge to a given root.
where xlimits and ylimits are 1 × 2 arrays with the limits for each axis. These values
determine the viewing window. You will have to modify the code slightly to allow for
variable axes limits. You might also want to include resolution as an input argument to
create images with higher resolution.
(c) Use your julia function to produce a plot with a view going from −0.3 ≤ x ≤ 0.3 and
−0.25 ≤ y ≤ 0.25. Also make sure the resolution is high enough so that the resulting
plot does not look “pixelated”. Include this plot in your write-up. Put the julia.m file
in your homework 6 dropbox folder, do not include it in your write-up.
3. (Using fzero) Consider the Colebrook equation for the friction factor in a fully developed
pipe flow
1
/D
2.51
√ = −2 log10
√
+
3.7
f
ReD f
where f is the Darcy friction factor, /D is the relative roughness of the pipe material, and
Re is the Reynolds number based on a pipe diameter D. Use fzero to find the friction factor
f corresponding to parameter values /D = 0.0001 and ReD = 3 × 105 . Use a tolerance 10−8
with fzero. In your homework write-up give the value of f as well as the code you used to
solve the problem. Do not include the fzero code just how you called it.
Hint : Use the function optimset to set up the tolerance TolX.
4. (Using fzero) David Peters (SIAM Review, 1997) obtains the following equation for the
optimum damping ratio of a spring-mass-damper system designed to minimize the transmitted
force when an impact is applied to the mass:
h p
i
cos 4ζ 1 − ζ 2 = −1 + 8ζ 2 − 8ζ 4
Use fzero to find the ζ ∈ [0, 0.5] that satisfies this equation with a tolerance 10−12 . In your
homework write-up give the value of ζ as well as the code you used to solve the problem. Do
not include the fzero code just how you called it.
2
5. (Using fminbnd) The orbits of Mercury and Earth can be (ideally) parameterized with respect
to time as follows:
Mercury
xm (t) = −11.9084 + 57.9117 cos(2πt/87.97)
ym (t) = 56.6741 sin(2πt/87.97)
Earth
xe (t) = −2.4987 + 149.6041 cos(2πt/365.25)
ye (t) = 149.5832 sin(2πt/365.25)
The units on the position are in 106 km and the units on time are days. The coordinate system
has been arranged so that the sun is at the center.
(a) Use the function fminbnd to determine a time between 0 and 1000 for which the distance
between the Earth and Mercury is minimal. Report this time and plot the two orbits
together with the positions when the distance is at a minimum clearly marked on the
curves. Is the time you found the global minimum over the interval 0 ≤ t ≤ 1000, or just
a local minimum? Explain how you determined this with, for example, some kind of plot.
(b) Repeat part (a), but now find a time for which the distance is maximal. Report this value
and produce a similar plot to part (a) with the positions marked.
6. (Least Squares) NCM, Chapter 5, Problem 5.8
3