### Introduction To MATLAB Programming

## Fractals and Chaos

A fractal is a geometric figure that can be subdivided into parts that are mathematically similar to the whole. You will be asked to plot the Mandelbrot fractal, and effectively practice constructing while loops, which terminate based on a known and specified condition. You will also learn how to use commands that help you terminate the loop prematurely and otherwise modify the execution of the loop.

In this unit, you will also use comparative and logical operations to evaluate expressions and create a truth matrix composed of zeroes and ones. The truth matrix will be useful for logical indexing, which can be used in situations that require that certain elements of a list be separated from the rest of the elements based on a certain characteristic.

## Logistic Equation

**Homework 4.*** For r varying between 0 and 4, find out the possible “limit cycles” ^{§} of the iterative map:*

*This converges to a single value for some values of* r *but for others results in an “orbit”, which can be very long. For every* 0<r<4,* “find” this orbit and plot the orbits together.*

*Use the “matrix-at-a-time” notation we learned in the last iteration example:*

*Start with a vector of*r−*values, and a vector of*x−*values (both row vectors and the same size)*.*Perform many (how many?) iterations on the whole vector of*x−*values, so that each place in the vector is updated according to its appropriate*r.*Plot the resulting*x*-values vs. the*r*values*.*Continue the iteration and plot several more iterations (how many?).**Observe the nice pattern that arises, and explore its self-similarity properties.*

Hint 1: (Am I getting the right answer?) The result should look something like this:

Graphing an iterative function.

Hint 2: (Code takes forever) If your code is running very slowly, you should consider updating all the orbits (one for each r value) at once. This implies holding a vector that corresponds to the r values you are considering, and a vector corresponding to the various orbits, and with one line of code you can update all the values in the orbit. Do this in a loop to find “late” elements of the orbit.

Hint 3: (Getting a similar plot, but not quite) You should only plot the late elements, so perhaps iterate without plotting for some time (maybe 1000 iterations?) and then plot successive elements of the orbit (say 100?).

Hint 4: (I don’t have so many points in my plot) Be sure to use `hold on`

so that each plot doesn’t erase the previous ones.

More ideas:

- Find how to make the plot have small dots as markers.
- Can you allow the user to “zoom in” on your plot? Once asked to see a region smaller than (0,4)×(0,1) you should probably increase the “density” of your r measurements, and confine the plotting of the points so that only the requested x‘s are plotted.

^{§}a limit cycle is an orbit of an iterative map that the dynamics of the problem converges to, regardless of the initial condition.

## More on Loops

The `for`

loop is a very useful tool for doing things over and over again for a certain number of times already known *in advance*. There are two possibilities that we would like to consider:

- What if we do not know in advance how many iterations we will need?
- What if we would like to stop a loop before it is due to end?

An example for the first kind would be a Newton iteration that should run until the value of f(x) is “small” enough, for example 10−12. Before actually performing the iterations we do not know how many steps it will take, so a `for`

loop is not exactly the right type of loop. We could get around this limitation if we introduce a maximum number of allowed iterations and then use the (as-of-yet unknown) mechanism for terminating a loop prematurely once we find a good enough approximate root.

A `while`

loop tells MATLAB® to continue iterating as long as a certain condition (which you specify) is satisfied. The syntax is:

`while`

`<condition> <statements>`

`end`

MATLAB evaluates the `<condition>`

and if it is true (or a non-zero number) it performs the `<statements>`

, if not, it continues after the `end`

. After each time it evaluates `<statements>`

MATLAB goes back and evaluates `<condition>`

again, etc. Note that `<condition>`

does *not*get evaluated in the middle of evaluating `<statements>`

but, rather, only before evaluating them. Here’s a simple way of adding two positive integers (very silly):

`x=5; y=6; while y>0 x=x+1; y=y-1;`

`end`

Of course, this fails miserably if `y`

is not a positive integer (doesn’t do anything, do you understand why?)

**Exercise 16 . **

*Solve the following problems using a*

`while`

*loop:*

*Show the numbers from 1 to 10**Show the numbers from 10 to -10**Find out how many divisors 28 has (*`mod`

*or*`rem`

*will be useful here)**Find out if a number is prime**Use an external*`while`

*and an internal*`for`

*loop to find the first 100 prime numbers.**A*perfect number*is a number*n*whose divisors (including 1 but excluding itself) add up to*n*itself. For example, 6 is a perfect number. Check if a number is perfect.**Use two nested*`while`

*loops to find the first 3 perfect numbers.*

**Homework 5.** *Consider the following sequence defined completely by the first element* S1^{¶}:

*A still*^{||} *open question in mathematics is whether all such sequences always arrive at 1 for large enough* n *(the alternatives being that some sequences may rise indefinitely, or that there may be a closed orbit that does not include 1). Compute the number of iterations it takes to arrive at* 1 *given a starting value* s *using a while loop. Since we do not know how long it will take to arrive at 1 (though you can assume that it will happen eventually) we might want to construct this sequence using a while-loop. What starting number smaller than 10,000 has the longest trajectory? What’s the largest number on that trajectory?*

^{§}This is the subject of the Collatz Conjecture.

^{||}Despite a recent “near” solution.

## Terminating a Loop Prematurely: Break and Continue

As you may recall, a `while`

loop will evaluate all its statements without checking the condition. Similarly a `for`

loop will run through all of its iterations. The `break`

keyword tells MATLAB® to exit the loop immediately. It will only terminate one loop (in the case of nested loop, the innermost one it is in) and will normally be protected by an `if`

statement (otherwise the loop is silly). Here is an example that computes the “trajectory” of 6 but stops if it finds a 17 in it:

`s=6; % initialize s to 6 while s~=1 % as long as s is not equal to 1 stay in the loop if s==17 % if s equals 17 sprintf('Found 17 in the loop!!') break; end if mod(s,2) % the actual "brains" of the iteration s=s/2; else s=3*s+1;`

`end end`

The keyword `continue`

is similar but different. It avoids the rest of the statements of the inner most loop, but continues in the loop (does not stop like `break`

).

Here’s example code for a `while`

loop that uses both break and continue to find the first 100 primes (not very efficiently, but it’s only an example):

```
n=1;
m=0;
while 1 % this means that unless we use "break", the loop will continue "forever"
n=n+1; % increase n
flag=0; % reset flag
for i=2:ceil(sqrt(n)) % no need to check numbers greater than the square-root of n
if mod(n,i)==0 % means that i divides n exactly
flag = 1 % to know that we found a divisor
break; % no need to remain in the for loop
end
end
if flag
continue % to avoid the next line. It could have also been done
% differently with an "if" statement, but this can be more elegant
end
sprintf('%d is prime!\n',n) % this is quite an interesting command...
% take some time to learn about it
m=m+1; % increment primes count
if m>=100 % if we have enough
break; % stop looking for primes
end
end
```

**Homework 6.** *The keywords* `break`

*and* `continue`

*are not “needed”* per se, *but they can make the code more elegant and readable. Rewrite the above code for the first 100 primes without using neither* `continue`

*nor* `break`

.

**Homework 7.** `for`

*loops and* `while`

*loops are not inherently different:*

*The “input” of a*`for`

*loop is a variable and a vector of values. Recreate the functionality of a*`for`

*loop using a*`while`

*loop*.*The “input” of a*`while`

*loop is the condition statement. Recreate the functionality of a*`while`

*loop using a for loop. (Hint: when using the notation*`for i=1:n`

*MATLAB does not*actually*create the vector*`1:n`

.*Internally it simply iterates over the values in that vector by incrementing*`i`

*until it reaches*`n`

.*This means that if you write*`for i=1:281474976710655`

*you’ll get a loop that, on its own, will “never” terminate. Explanation: 281474976710655 is the largest*integer*that MATLAB can represent internally. It is such a large number that even if every pass through the loop only takes 1 millisecond getting through the loop will take about 10000 years.)*

## Truth Statements and Logical Indexing

The Mandelbrot Set is the set of points z∈ℂ for which the sequence

remains bounded as n→∞. Color can be added for the unbounded elements by specifying how “fast” they diverge, for example, how many iterations it takes to reach some large absolute value. (Since once |xn| is large enough the sequence will be growing indefinitely.)

To decide on the color of a particular point on the screen (x,y), we define a complex number z=x+iy and iterate according to the equation mentioned above.

**Exercise 17.** *With three nested* `for`

*loops, iterate over* `x=-1.5:.01:0.5`

*and* `y=-1:0.01:1`

. *For each point, iterate, say, 100 times and*`break`

*if* `abs(z_n)`

*is ever too large, say larger than 10. There’s no need to keep the sequence of* z‘*s, but be sure not to overwrite the* z *variable that comes from* x+iy. *If the the sequence grew to be large, leave white, otherwise place a black dot at (the original) *(x,y).

The result should look like this:

Graphing the Mandelbrot set.

Notice that it takes quite a long time to run, especially if you modify the code to give you a plot with better resolution. To fix this we want to use matrix-at-once calculation. We can create matrices X and Y using `meshgrid`

:

`x=-1.5:.1:1; y=-1.5:.1:1.5; [X,Y]=meshgrid(x,y);`

`% make sure you understand what this function does!`

Each element in these matrices holds the appropriate x- or y-value for that position in the matrix. We can now calculate the iterations very simply:

`Z=X+1i*Y; % use 1i since you can overwrite i but 1i will always be sqrt(-1) Zn=0*Z; % Start with zeros everywhere Zn=Zn.^2+Z;`

`% One step of the iteration, notice the point-wise power operator.`

Now that we know how to iterate we need to figure out how to color the different points. Here’s where the name of this subsection comes in.

Just as we can perform point-wise arithmetic with a matrix A, we can also perform point-wise *comparative* and *logical* operations. The comparative operators are: `>`

(greater than), `<`

(less than) `==`

(equal to**), `~=`

(not equal to), `>=`

(greater than or equal to) `<=`

(less than or equal to), while the logical operators are: `~`

(not), `&`

(and), `|`

(or). For more help on operators type `help ops`

.

Letting `A=magic(4)`

, understand the following constructions:

`A>5`

`A<=3 | A>=10`

`A==10 | A==5`

`A>5 & A<10`

`~(A>5 & A<10)`

Given that `x = [1 5 2 8 9 0 1]`

and `y = [5 2 2 6 0 0 2]`

, execute and explain the results of the following commands:

`x > y`

`y < x`

`x == y`

`x <= y`

`y >= x`

`x | y`

`x & y`

`x & (~y)`

`(x > y) | (y < x)`

`(x > y) & (y < x)`

Together with the functions `any`

and `all`

we can ask questions like “Is any/all element of A equal to 4?”, the functions look down the first “non-singleton” dimension by default, which means that they work as you expect for column and row vectors, but for matrices they only “reduce” the first dimension:

```
>> any(magic(4)==4) % by default reduces the first dimension (rows):
ans =
1 0 0 0
>> any(magic(4)==4,2) % we can tell it to reduce a different one:
ans =
0
0
0
1
>> any(any(magic(4)==4)) % To reduce a matrix to a single number, we must apply
% the functions twice:
ans =
1
```

In addition we can use a truth matrix to *access* elements of a matrix. If we have a matrix A and a truth matrix (a matrix found by comparative and logical operations) B *of the same size as A*, we may write `A(B)`

and the result will be the list (not matrix!) of elements of A for which the corresponding elements of B are “true.” This is slightly confusing, so an example is useful:

```
>> A=magic(4);
>> B=A>10;
>> A(B)
ans =
16
11
14
15
13
12
```

**Exercise 18.** *Explain the order of the numbers in the previous example.*

In addition to getting the list of elements that correspond to some truth matrix, we can also use the truth matrix to change specific elements of a matrix:

```
>> A
A =
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
>> A>5
ans =
1 0 0 1
0 1 1 1
1 1 1 1
0 1 1 0
>> A(A>5)=0
A =
0 2 3 0
5 0 0 0
0 0 0 0
4 0 0 1
```

**Exercise 19.** *The exercises here show the techniques of logical indexing. Given* `x = 1:10`

and `y = [3 1 5 6 8 2 9 4 7 0]`

, *execute and interpret the results of the following commands:*

1(a) |
` (x > 3) & (x < 8)` |

1(b) |
` x(x > 5)` |

1(c) |
` y(x <= 4)` |

1(d) |
` x( (x < 2) | (x >= 8) )` |

1(e) |
` y( (x < 2) | (x >= 8) )` |

1(f) |
` x(y < 0)` |

*Given* `x = [3 15 9 12 -1 0 -12 9 6 1]`

, *provide the command(s) that will*

2(a) |
… set the positive elements of x to zero. |

2(b) |
… set values of x that are multiples of 3 to 3 (`rem` will help here). |

2(c) |
… multiply the even elements of x by 5. |

2(d) |
… extract the values of x that are greater than 10 into a vector called y. |

2(e) |
… set the values in x that are less than the mean value of x to zero. |

2(f) |
… set the values in x that are above the mean to their difference from the mean. |

**Homework 8.** *Going back to the Mandelbrot set calculation, we can keep a matrix that informs us of the step at which each element has become too large. Let’s call it* `Iter`

*and agree that a zero value corresponds to* *“not too large yet” and a non-zero value will be the iteration number at which the element has become big (absolute value greater than 10). At each step we update only these elements of* `Zn`

*that are still small (those that correspond to* `~Iter`

*being true), then find out which are the new big elements (*`abs(Zn)>10`

*and* `~Iter`

*) and set the corresponding values of* `Iter`

*to the iteration number (from the for loop). Notice that using* `~`

*on a matrix of numbers returns* `true`

f*or the elements that are zero and* `false`

*for those that are not zero. After 100 iterations you need to plot the results. Since for the locations that did not converge we have the number of iterations it needs to reach size 10, we can make a nicer plot than the black-and-white plot from before. For this we can use:*

```
p=pcolor(X,Y,Iter);
set(p,'edgecolor','none')
axis equal
```

*We only need to plot once, after all the iterations have happened. Notice how much faster this is than the nested loops from before. Your result should look like this:*

Graphing the Mandelbrot set with logical indexing.

*Some ideas for “extra credit”:*

*Try zooming into places of interest; see how the fractal is “self similar”?**Change the number of iterations and see how things change.*

**Make sure you understand the difference between = and ==; they mean entirely different things, and writing one instead of the other will most likely result in a strange error.