1. Loops

1.1. What is the difference between the While and Do loops that were presented in the lecture notes?

1.2. Using the For, While and Do loops, write the Mathematica commands that  calculate 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.

1.3. Calculate the sum (∑_(i = 0)^10) (i^2) using Do- and While loops.

1.4. Calculate the values (Subscript[y, i])for i = 0, …, 10 when y is given by the recurrence
	(Subscript[y, i]) = 2 (Subscript[y, i - 1]) + 4  (Subscript[y, i - 2])
(Subscript[y, 0]) = 1
(Subscript[y, 1]) = 4.
Write the numbers (Subscript[y, i]) into a file called arvot.dat.
	
1.5. Advanced: 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.1

Aside from the obvious differences in the syntax of the commands, the examples in the lecture notes differ in the number of terms that are accumulated into the sum (the following commands are copied from the lecture notes):

For[a = 0; i = 0, i < N[Pi], i += N[Pi/5], a += Sin[i]; 
      Print[a]]

The above command calculates the sum S = 
        sin(0) + sin(π/5) + sin(2  π/5) + … + 
          sin(4  π/5). The term sin(π) is excluded, because of the loop condition i < π.

a = 0
Do[a += Sin[i]; Print[a], {i, 0, N[Pi], N[Pi/5]}]

This do-loop calculates the sum S = 
        sin(0) + sin(π/5) + sin(2  π/5) + … + 
          sin(π). This time the last term missing from the for-loop is included.

a = 0
i = 0
While[i < N[Pi], i += N[Pi/5]; a += Sin[i]; Print[a]]

Above command calculates the sum S = 
        sin(π/5) + sin(2  π/5) + … + sin(π). The term sin(0) is missing, because iis incremented before sin(i) is added to a. For that same reason, the term sin(π) is included: i gets incremented before it's used, and therefore even though the loop condition is i < π, the term sin(π) comes into the sum.

1.2.

Here are a couple of alternatives: Firstly, all three commands using just three variables. In all examples, the variable i holds the index to the Fibonacci number that's being calculated on that particular iteration. Since Fibonacci numbers (Subscript[x, 0]) and (Subscript[x, 1]) are given, we start with i = 2. When an iteration of the loop starts, the variables xm holds the previous Fibonacci number ((Subscript[x, i - 1])) and xmm the one before that ((Subscript[x, i - 2])). The next Fibonacci number ((Subscript[x, i])) will be stored in variable x. Finally, x is printed. Loop conditions are set so that the last value printed will be the nth Fibonacci number, (Subscript[x, n]).

n = 10

For[xmm = 0; xm = 1; i = 2, i ≤ n, i++, x = xmm + xm; 
      xmm = xm; xm = x; Print[x]]

xmm = 0;
xm = 1;
Do[x = xmm + xm; xmm = xm; xm = x; Print[x], {i, 2, n}]

xmm = 0;
xm = 1;
i = 2;
While[i ≤ n, x = xmm + xm; xmm = xm; xm = x; 
      Print[x]; i++]

Here's a more straightforward and intuitive way of doing the same thing:

Clear[
      f]

    For[f[0] = 0; f[1] = 1; i = 2, i ≤ n, i++, 
      f[i] = f[i - 1] + f[i - 2]; 
      Print[f[i]]]

  

    Clear[f]

    f[0] = 0;

    f[1] = 1;

    Do[f[i] = f[i - 1] + f[i - 2]; 
      Print[f[i]], {i, 2, n}]

  

    Clear[f]

    f[0] = 0;

    f[1] = 1;

    i = 2;

    While[i ≤ n, f[i] = f[i - 1] + f[i - 2]; 
      Print[f[i]]; i++]

(Advanced: The above commands actually define  functions. Usually, when we define a function, we use the notation f[x_] = expression (note the underscore _ in the argument). However, it is also possible to define functions by assigning it a value for some specific value of the argument, for example the commands f[1] = 1; f[2] = 3 would define the function f so that f(1) = 1 and f(2) = 3. )

1.3

Straight-forward:

a = 0;
Do[a += i^2, {i, 0, 10}]
a
a = 0;
i = 0;
While[i ≤ 10, 
        a += i^2; i++];
a

1.4.

See the Fibonacci numbers example above for comparison.

chan = OpenWrite["arvot.dat"];
For[ymm = 1; ym = 4; i = 0, i ≤ 10, i++, 
      y = 2  ym + 4  ymm; ymm = ym; ym = y; 
      Write[chan, y]]
Close[chan]

1.5.

The following piece of code iterates the Hénon map 5000 times. The values of x and y from the previous iteration are kept in variables xm and ym, and on each iteration, the following values are computed from the given expression. Then, the new pair (x, y) is appended to the list.

xm = 0. ;
ym = 0. ;
a = 1.4;
b = 0.3;
n = 5000;
hlist = {{xm, ym}};
For[
i = 1,    
i ≤ n,    
i++,    
x = 1 + ym - a xm^2; 
      y = b xm;    
xm = 
        x; ym = y;   
hlist = Append[hlist, {x, y}]   
]
ListPlot[hlist, PlotStyle -> PointSize[0.002]]

2. Schrödinger equation

2.1 Try to find the first excited state of the hydrogen atom and the corresponding wavefunction.

2.2. Advanced - do this only if you have time. When solving the eigenvalue E we do not necessarily need to write a loop ourselves. We are only looking for a root of the function Δ(E), so, if we can write Δ as a Mathematica function, we can use the FindRoot command to locate the energies E. Try to do this. You can use the Module command to wrap the commands for calculating Δ into a function. Module command takes two arguments, first a list of local variables (that is, variables that you only use inside the function) and second, the commands separated by semicolons (;). Inside the Module command, you can use the Return -command to return the value for delta. Your function Δ in Mathematica could look something like this:

Delta[e_] := Module[
	{rata, ratl, …},
	rata=NDSolve[{yhtalo, …];
	…
	ero = logalk[rvali] - loglop[rvali];
	Return[ero]];
	
You can now plot the function Delta, and use FindRoot to get the energy eigenvalues!

Solutions

2.1. This is a part of the last homework assignment; you'll just have to do this one yourself!

2.2. I'll give only the definition for the function Delta.

Clear[Delta];
Delta[e_] := 
      Module[{}, 
Clear[y, 
            r]; 
yhtalo = y''[r] - (-2/r  - e) 
                  y[r] == 0; 
c1 = 
            1; 
c2 = 
            10^-4; 
yal[r_] = 
            c1*r; 
dyal[r_] = 
            D[yal[r], r]; 
ylo[r_] = 
            c2*E^(-Sqrt[-e]*r); 
dylo[r_] = 
            D[ylo[r], r]; 
tark = 
            10^(-12) ; 
ralku =  
            10^(-8) ; 
rloppu  = -Log[
                      eps]/(Sqrt[-e]) /. {eps -> 
                  1*^-10};  (* <- 
                  this guarantees that the tail of the function is decaying at rloppu *) 
rvali = 
            rloppu/15; 
rata = 
            NDSolve[{yhtalo, 
                y[ralku] == yal[ralku], y'[ralku] == 
                  dyal[ralku]}, y[r], {r, ralku, rvali}, 
              AccuracyGoal -> 15]; 
alku[r_] = 
            y[r] /. rata[[1]]; 
logalk[r_] = 
            D[alku[r], r]/alku[r]; 
ratl = 
            NDSolve[{yhtalo, 
                y[rloppu] == ylo[rloppu], y'[rloppu] == 
                  dylo[rloppu]}, y[r], {r, rloppu, rvali}, 
              AccuracyGoal -> 15]; 
loppu[r_] = 
            y[r] /. ratl[[1]]; 
loglop[r_] = 
            D[loppu[r], r]/loppu[r]; 
ero = 
            logalk[rvali] - 
              loglop[rvali]; 
ero
] /; 
        NumericQ[e]
 (* 
      This last command tells Mathematica NOT to evaluate the function before it has a numerical argument *) 

    
Plot[Delta[e], {e, -1.1, -0.9}]
FindRoot[
      Delta[e] == 
        0, {e, -1.1, -0.9}]