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
	(Subscript[x, i]) = (Subscript[x, i - 1]) + (Subscript[x, i - 2]),  
      i ≥ 2
(Subscript[x, 0]) = 0
(Subscript[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 (∑_(i = 0)^10) (i^2) using Do- and While loops.
4. The Hénon map.  Make a loop that iterates the famous Hénon map
	(Subscript[x, i + 1]) = 1 + (Subscript[y, i]) - a ((Subscript[x, i])^2)
(Subscript[y, i + 1]) = b (Subscript[x, i]), a = 1.4
b = 0.3.
Try if you can store the results in a list of points ((Subscript[x, i]), (Subscript[y, i])) and plot those! (You can use the Append command to add (x, y) pairs to a list)

Solutions

1

The Do loop is somewhat similar to the Table command, where a definite range is given over which to carry out the commands. The While loop is more of a traditional loop, where a test is performed on each step.

i = 0; n = Input["Give N"]; a = {};
While[i < n, 
        a = Append[a, {i, N[Sin[i]]}]; i++];
Print[a // TableForm]

i = 0; n = Input["Give N"]; a = {};
Do[
        a = Append[a, {i, Sin[i] // N}], {i, 0, n, 
          1}];
Print[a // TableForm]

2

n = 
        Input["Calculate the first n Fibonacci numbers. Give n:"];
(Subscript[x, pp]) = 0; (Subscript[x, p]) = 1; Fib = {(Subscript[x, pp]), (Subscript[x, p])};
For[i = 0, i ≤ n, i++, (Subscript[x, 0]) = (Subscript[x, pp]) + (Subscript[x, p]); 
      (Subscript[x, pp]) = (Subscript[x, p]); (Subscript[x, p]) = (Subscript[x, 0]); 
      Fib = Append[Fib, (Subscript[x, 0])];]
Print[TableForm[Fib]]

3

S = 0;
Do[S += (i^2), {i, 0, 10, 1}]
Print[S]; S = 0; i = 0;
While[i ≤ 10, S += (i^2); i++]
Print[S];
Sum[(i^2), {i, 0, 10, 1}]

4

Note the typo... The initial values were not given in the assignment. You can try out these, although you can give pretty much anything as long as it's small enough. Also try to vary the parameters α and β.

(Subscript[x, p]) = 0.0; (Subscript[y, p]) = 0.0; 
    Henon = {{(Subscript[x, p]), (Subscript[y, p])}}; α = 1.4; β = 0.3; 
    n = 5000;
For[i = 0, i < n, i++, (Subscript[x, 0]) = 1 + (Subscript[y, p]) - α ((Subscript[x, p])^2); 
      (Subscript[y, 0]) = β (Subscript[x, p]); (Subscript[x, p]) = (Subscript[x, 0]); (Subscript[y, p]) = (Subscript[y, 0]); 
      Henon = Append[Henon, {(Subscript[x, 0]), (Subscript[y, 0])}];]
ListPlot[Henon]

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 (Subscript[x, 0]), that is, you want to compute
	\(TraditionalForm\`x\_n = \(\(\(f(f(⋯\ \(f(f(x\_0))\) ⋯))\)\+⏟\)\+\(n\ applications\ of\ f\) = \(\(\((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 (Subscript[x, i]) for its solution, Newton's method gives you an improved guess (Subscript[x, i + 1]) for the solution: 
	(Subscript[x, i + 1]) = f((Subscript[x, i])) = 
          (Subscript[x, i]) - ((g((Subscript[x, i])))/(g' ((Subscript[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)/(2  x)).
So, set a = 10, and apply the Nest command to f(x). Set the initial guess (Subscript[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 -> ∞))(Subscript[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 (Sqrt[10]).

Solutions

1

The Nest and FixedPoint commands are extra material. Especially the loops of the previous excercise are much easier to implement using Do, While, etc.

n = Input[];
Nest[{#[[1]] + #[[2]], Part[#, 1]} &, {1, 0}, 
      n]
Nest[Append[#, Last[#] + Part[#, -2]] &, {0, 1}, n]

n = Input[];
Nest[Append[#, Length[#]^2] &, {0}, n]
Sum[%[[i]], {i, 1, Length[%]}]

This is the Henon map. The Nest function performs OK with this one.

Nest[
        Append[#, {1 + Last[#][[2]] - α Last[#][[1]]^2, β Last[#][[1]]}] &, {{0, 0}}, 
        1000];
ListPlot[%]

2

\(Nest[# - \(#\^2 - 10.0\)\/\(2  #\) &, 3.0, 
      10]\)
(Sqrt[10]) // N

\(FixedPoint[# - \(#\^2 - 10.0\)\/\(2  #\) &, 3.0]\)

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 (Subscript[a, x]) = (Subscript[F, thrust,  x]) + (Subscript[F, G,  x])
m (Subscript[a, y]) = (Subscript[F, thrust,  y]) + (Subscript[F, G,  y])
Assume again that the thrust (Subscript[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
	\(TraditionalForm\`F\&_\_thrust = \(\((F\_\(thrust, \ x\), \ 
            F\_\(thrust, \ y\))\) = \(F\_thrust\)(
            v\_x\/\(\(|\)\(v\)\(|\)\), \ v\_y\/\(\(|\)\(v\)\(|\)\))\)\), 	
where
	|v| = (Sqrt[((Subscript[v, x])^2) + ((Subscript[v, y])^2)]),	({{\(v\_x = x\& . \)}, {\(v\_y = y\& . \)}}).
Earth's gravity is naturally directed towards the center of the planet, so gravity reads
	\(TraditionalForm\`F\&_\_grav = 
        m\ g\ \((R\/\(\(|\)\(r\)\( | \^2\)\))\) 
          r\&_\/\(\(|\)\(r\)\(|\)\)\),
where r = (x,  y),  
      and  | r |  = (Sqrt[(x^2) + (y^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

First, load the necessary constants.

<< Miscellaneous`PhysicalConstants`
R = EarthRadius/Meter // N;
g = 
        AccelerationDueToGravity/(Meter/(Second^2)) // N;

Then, the parameters

m0 = 1000; f0 = 28000; c = 1; p = 600;

Define the time-dependent parameters, i.e. mass and thrust force.

m[t_] := 
        If[c t < p, m0 - c t, m0 - p];
fm[t_] := If[c t < p, f0, 0];

Define the position and velocity vectors. All forces will be presented in terms of these.

Clear[r, t, x, y, f, v];
r[t_] = {x[t], y[t]};
v[t_] = {x'[t], y'[t]};
f[t_] = fm[t] v[t]/Norm[v[t]];
grav[t_] = g (R^2) r[t]/Norm[r[t]]^3;

A differential equation for both directions.

dex = 
        m[t] x''[t] == 
          Part[f[t], 1] - m[t] Part[grav[t], 1];
dey = 
        m[t] y''[t] == 
          Part[f[t], 2] - m[t] Part[grav[t], 2];

A fancy way to determine the initial conditions.

latitude = 90;
x0 = R Cos[latitude*Pi/180];
y0 = R Sin[latitude*Pi/180];
vx0 = 0.001;
vy0 = 0.1;

Now, commence with the solution.

T = 2000;
Sol = 
        Part[NDSolve[{dex, dey, x[0] == x0, x'[0] == vx0, 
              y[0] == y0, y'[0] == vy0}, {x[t], y[t]}, {t, 0, 
              T}], 1];
flight = 
        ParametricPlot[{x[t] /. Sol, y[t] /. 
Sol}, {t, 0, 
            T}, DisplayFunction -> Identity];

Draw a picture of the Earth to get the scale of the trajectory visible.

Earth = Graphics[Circle[{0, 0}, R]];

Show[flight, Earth, DisplayFunction -> $DisplayFunction, 
        PlotRange -> All, 
        Ticks -> None];
Plot[Norm[r[t] /. Sol], {t, 0, T}, Ticks -> None]

Mandelbrot set

The Mandelbrot set is defined as the set of points c for which the sequence defined by
(1)	 (Subscript[z, 0]) = 0
(Subscript[z, i + 1]) = ((Subscript[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

The Module command can be used to wrap a set of commands into a single function. The first argument is the declaration of any possible local variables.

nmax = 100; Clear[x, y, f, n, z, n0]
f[x_, y_] := 
        Module[{n}, n = 0; 
          For[z0 = 0; z = 0, Norm[z] < 2 && n < nmax, z0 = z, 
            z = (z0^2) + x + I y; n++]; 
          Return[n]];
Plot3D[f[x, y], {x, -2, 1}, {y, -1, 1}, 
      PlotPoints -> 100]