Excercise 8

1. Programming

1.1. Write a program that defines a function f(x) = sin(10 cos(x)) and plots it.

1.2. Make a program that solves the equation of a damped harmonic oscillator
	\(TraditionalForm\`m\ \(x\& .. \) + γ\ x\& .  + k\ x = 0, \ 
      x(0) = s, \ \(x\& . \)(0) = 0\).
Write the program into a separate file. The program should set values for m, γ, k and s, then solve the equation and finally plot the solution.

Solutions

All you have to do is to write the following stuff into a file and then run it using the << filename syntax.

Clear[f, 
      x]

    f[x_] = Sin[10 Cos[x]]

    Plot[f[x], {x, 0, 2  Pi}]

    Clear[f]

  

    m = 1

    k = 1

    g = 1

    s = 1

    sol = 
      DSolve[{m x''[t] + g x'[t] + k x[t] == 0, 
          x[0] == s, x'[0] == 0}, x[t], t]

    Plot[x[t] /. sol, {t, 0, 10}]

    Clear[m, k, g, s, sol]

2. The Rocket revisited

2.1. Mess with the parameters of the rocket example.

3. Piecewise functions, conditional evaluation

3.1. Define a function f as
	f(x)= {({{0,  x ≤ 0}, {1,  x > 0}})
Use both the /; method and the If command.

3.2. Define a piecewise function that gets as its value (x^2) when x < 3, x + 5 when 3 ≤ x < 7 and 14 when x ≥ 7. Plot the function f. How could you redefine the function so that it would be continuous?

Solutions

3.1 Here are the two ways, (Subscript[f, 1])and (Subscript[f, 2]):

Clear[f1, f2]
f1[x_] := 0 /; x ≤ 0
f1[x_] := 1 /; x > 0
f2[x_] := If[x ≤ 0, 0, 1]

3.2

Clear[f, g]
f[x_] := x^2 /; x < 3
f[x_] := x + 5 /; 3 ≤ x < 7
f[x_] := 14 /; x ≥ 7
g[x_] := x^2 /; x < 3
g[x_] := x + 5 + 1 /; 3 ≤ x < 7
g[x_] := 14 - 1 /; x ≥ 7
Plot[{f[x], g[x]}, {x, 0, 10}, 
      PlotStyle -> {Hue[0], Hue[1/3]}]

Advanced

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.

A.1. Implement the loops of previous excersise using the Nest command.

A.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]).

Rocket science

A.3. 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.

Über-advanced: 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.