Excercise sheet 8 -- solutions

Loops

1. What is the difference between the While and Do loops that were presented in the lecture notes?
2. Write the Mathematica commands to calculates the n first Fibonacci numbers, defined according to the relation
    x_i = x_ (i - 1) + x_ (i - 2), i≥2  x_0 = 0  x_1 = 1
ie, each Fibonacci number is the sum of two previous ones. Ask the number n from the user. To print out the numbers you can use the Print-command.
3. Calculate the sum Underoverscript[∑, i = 0, arg3]i^2 using Do- and While loops.
4. The Hénon map.  Make a loop that iterates the famous Hénon map
    x_ (i + 1) = 1 + y_i - a x_i^2  y_ (i + 1) = b x_i, a = 1.4  b = 0.3.
Try if you can store the results in a list of points (x_i,y_i) and plot those! (You can use the Append command to add (x,y) pairs to a list)

Solutions

Nest and FixedPoint

The Nest and FixedPoint commands are useful when one needs to iterate a function. Say you have a function f, and you need to apply it n times to a number x_0, that is, you want to compute
    x_n=Underscript[Underscript[f(f(⋯ f(f(x_0)) ⋯)), ︸], n applications of f]=Underscript[Underscript[(f ◦ f ◦ ⋯ ◦ f), ︸], n](x_0).
This is called nesting and to evaluate these in Mathematica, you can use the Nest command:
    Nest[f, x0, n],
where f is your function (must be a pure function, not a Mathematica expression!), x0 the argument for the innermost f, and n is the number of iterations.

1. Implement the loops of previous excersise using the Nest command.
2. Find the square root of 10 using the Newton's method and the Nest-command. If you've an equation g(x)=0, and a guess x_i for its solution, Newton's method gives you an improved guess x_ (i + 1) for the solution:
    x_ (i + 1)=f(x_i)=x_i-g(x_i)/(g ' (x_i)).
When f is applied repeatedly, you end up with better and better approximations for the solution (we're ignoring the important matter of convergence here, of course).
To find the square root of a, the function g(x)=x^2-a, and
(1)    f(x)=x-(x^2 - a)/(2x).
So, set a=10, and apply the Nest command to f(x). Set the initial guess x_0 to your liking and try different values of n. Compare the results to Sqrt[10].

The command FixedPoint works the same way as Nest, but you don't give the number of iterations as an argument. Instead, Mathematica applies the function f so many times that the result does not change, i.e. it tries to find the limit
    x=lim_ (n∞)x_n,
This x is called a fixed point, since it satisfies f(x)=x, hence the name of the command. Try the FixedPoint command with the function f from Eq. (1). Compare the result to 10^(1/2).

Solutions

1

2

Rocket science

Make the rocket example two dimensional, ie. take into account that the rocket can move in two directions. With some simplifications we can take the equations of motion to be
    m a_x = F_ (thrust, x) + F_ (G, x)  m a_y = F_ (thrust, y) + F_ (G, y)
Assume again that the thrust F_thrust is constant. Let us take the direction of this force to be the same as the velocity of the rocket. In this case, the force vector reads
    Overscript[F, _] _thrust=(F_ (thrust, x), F_ (thrust, y))=F_thrust(v_x/(| v |), v_y/(| v |)),     
where
    |v|=(v_x^2 + v_y^2)^(1/2),    

v_x = Overscript[x, .]
v_y = Overscript[y, .]
.
Earth's gravity is naturally directed towards the center of the planet, so gravity reads
    Overscript[F, _] _grav=m g (R/(| r |^2))Overscript[r, _]/(| r |),
where r=(x, y), and |r|=(x^2 + y^2)^(1/2). Note that we've taken the positions to be relative to the center of the mass, not to the surface of the planet.

Solution

Mandelbrot set

The Mandelbrot set is defined as the set of points c for which the sequence defined by
(1)     z_0 = 0  z_ (i + 1) = z_i^2 + c
does not approach infinity. Construct a function that counts the number of iterations of (1) that are required for the absolute value of z to become larger than 2. You also need to limit maximum number of iterations allowed. Then make a 3D plot of the function by taking the real and imaginary parts of c to be the x and y coordinates of the plot. The result should be a rather familiar looking object.

Solution

Created by Mathematica  (April 10, 2007) Valid XHTML 1.1!