Section 1: Basin Boundaries.  The Julia set.

1.1 Basin of attraction.  In the study of various dynamical systems, we attempt answer the question about the eventual fate of the system, or its asymptotic state.  As we have seen in the previous Labs, the asymptotic states can be stable (or attracting) fixed point, periodic and chaotic orbits, or infinity, which can also be considered an attracting fixed point.  Very often dynamical systems have several attracting asymptotic states, so orbits starting from different initial conditions may converge to different asymptotic states.  Therefore, in order to answer the question about the eventual fate of the system, not only do we need to find all possible asymptotic states, but also to determine which initial conditions converge to each asymptotic state.  To this end, we define the basin of attraction of an asymptotic state:

Definition: The set of all initial conditions from which orbits converge to a given asymptotic state is called the basin of attraction of that state.

1.2 Basin boundaries.  In the above example,  the two basins are separated by points x = -1, and x =1, which do not belong to either of the two basins.  Indeed the point x = 1 is itself an unstable fixed point, while x = -1 is mapped to the fixed point x = 1 on the next iterate, and so it is called an eventually fixed points.   In general, orbits starting from points exactly on the basin boundaries, always remain on the boundaries. Identifying the basin boundaries allows us to determine the eventual state of the orbit starting from any initial condition.

In the example above, the basin boundaries consist of only two points and are, therefore very simple.  However, there are many dynamical systems where the basin boundaries may have a very complicated fractal structure.

1.3 Complex maps. Very often 2-dimensional maps can be written as maps of a single complex number z = x + iy, where i is the imaginary unit and x and y are two real numbers representing real and imaginary parts of z.

Q1 Show that the complex map z_n+1 = z_n² is equivalent to the 2-dimensional map x_n+1 = x_n² - y_n² ;  y_n+1 = 2x_n y_n where z_n = x_n + iy_n.

Q2 The map z_n+1 = z_n² has a stable fixed point at z = 0. Show that the basin of attraction of this fixed point is the interior of the unit circle |z| < 1 in the complex plane. Show that any orbit starting outside the unit circle escapes to infinity.
Hint: consider the absolute value of the map iterates |z_n|.

Definition: The Julia set of a complex function h(z) is the boundary of the set of points whose orbits escape to infinity under the iteration of the map  z_n+1 = h(z_n).  In other words, the Julia set is the boundary between the escaping orbits and those that never escape.

It turns out that such a simple structure of the Julia set is exceptionally rare.  Typically, the Julia sets have much more complex structure, which is often a fractal.

1.4 The Julia set of the function z² + c. Consider the function h(z) = z² + c, where c = a + ib is a constant complex number.  Let us use Matlab to determine the Julia set of this function.  The nice thing about Matlab is that it understands complex numbers just as well as it does the real ones.  The basic idea is to consider points on a dense grid in the complex plane as starting points for orbits of the map

Program 1

c = 0.36 + i*0.1;
x = -2:0.02:2; y = -2:0.02:2;
[x0, y0] = meshgrid(x,y);
n = zeros(size(x0));
for k = 1:length(y),
  for m = 1:length(x),
    z = x0(k,m) + i*y0(k,m);
    for itr = 1:30,
      z = z.^2 + c;
      if abs(z) > 10,
        n(k,m) = itr; break,
      end,
    end,
  end,
end,
plotcolour(x0,y0,n);
pixgrid(size(x0));
 

Comments to the program:

• In the 1st line of the program we specify the value of the parameter c.
• 2nd line defines the values of x and y coordinates for the grid points.
• The function meshgrid creates arrays with coordinates of the grid points.
• 4th line defines an array n where the number of iterates required for an orbit to escape outside the circle  of radius 10 is stored.
• Next, we have two loops that run over the values of the grid coordinates in the x- and y-directions.
• The initial value of z is specified by the grid point (x0, y0).
• The third loop runs over the map iterates.
• If |z| > 10, the orbit is outside the circle of radius 10, the number of iterates is saved in the array n and the loop is terminated using the break command.  If the orbit does not escape after 30 iterates, the value of n remains equal to zero.  Such points are coloured black.

Q3   Check that the Julia set of the function with c = 0 is indeed a unit circle.   Construct plots of the Julia set for parameter values, c = 0.255, c = -1, c = 0.360284 + i*0.100376, c = -0.1 + i*0.8.  You can also modify the values of the grid coordinates to zoom in on a particular region of the set.  For example, try c = 0.4 + i*0.1, and x = 0:0.004:1; y = 0:0.004:1;
Note: No report is required for this exercise.  Simply explore various shapes of the Julia set and appreciate their complexity.

Section 2: The Newton-Raphson method.

2.1 The Newton-Raphson method for functions of a single variable.  Consider the problem of locating roots of some function f(x).  That is, we want to find solutions to the equation

There are many functions where such roots cannot be found analytically.  However, we can still look for such solutions using numerical procedures.  One such procedure was introduced by Isaac Newton over 300 years ago and successfully used long before the advent of computer age.  The idea is to start from an initial point x₀ which is reasonably close to the solution.  Then a line tangent to the function f at the point x₀ a is extended until it crosses the horizontal axis.  The crossing point x₁ is used as the next estimate for the solution.  We can easily determine x₁ from the fact that the slope of the tangent line is, on one hand, equal to the ratio f(x₀)/(x₀ - x₁) and, on the other hand, to the derivative of the function f'(x₀).  Thus

Repeating the same procedure at the point x₁, we obtain the next estimate of the solution according to the formula

It is easy to see that the procedure of the Newton-Raphson method for locating roots of the function f(x) is nothing else but a map  x_n+1 = N(x_n), where N(x) = x - f(x)/f'(x).  To understand how the Newton-Raphson method locates the roots of the function f(x), we can study the map N using methods of the dynamical systems theory.

Q4 Show that every root x^* of the function f is a fixed point of the map N.  Show also that every fixed point is superstable, that is, that N'(x^*) = 0.

Program 2

x = -3:0.03:3; y = -3:0.03:3;
[x0, y0] = meshgrid(x,y);
n = zeros(size(x0));
for k = 1:length(y),
  for m = 1:length(x),
    z = x0(k,m) + i*y0(k,m);
    for itr = 1:30,
      z = z - (z.^3 - 1)./(3*z.^2);
      if abs(z - 1) < 1e-5,
        n(k,m) = itr; break,
      end,
    end,
  end,
end,
plotcolour(x0,y0,n);
pixgrid(size(x0));

Comments to the program:

• The difference between Programs 1 and 2 is essentially in two lines: computation of the iterates of the map N(z) = z - f(z)/f'(z), and the condition for terminating the loop when the iterates converge to the solution z = 1 to within 10^-5.

Q5 How to modify Program 2 so that it plots the basin of attraction of the root z = 0 of the function f(z) = z - z⁴ ?