Nonlinear Dynamics and Bifurcation Theory with Applications using MATLAB
Dynamical Systems with Applications using MATLAB
Dynamical systems are mathematical models that describe how the state of a system changes over time. They are widely used to study phenomena in science, engineering, economics, biology, and many other fields. Dynamical systems can be classified into two types: discrete and continuous. Discrete dynamical systems have state variables that change at discrete time steps, while continuous dynamical systems have state variables that change continuously with respect to time. In this article, we will explore both types of dynamical systems using MATLAB, a powerful software for numerical computation and visualization. We will also see how to use Simulink, an add-on tool for MATLAB that allows us to simulate dynamical processes in a graphical way.
dynamical systems with applications using matlab
Linear Discrete Dynamical Systems
A linear discrete dynamical system is one that can be written in the form x(n+1) = Ax(n), where x(n) is a vector of state variables at time step n, A is a matrix of coefficients that determines how the state variables interact with each other, and x(n+1) is the vector of state variables at the next time step. A linear discrete dynamical system can be easily implemented in MATLAB using the following commands:
% Define the matrix A A = [0.8 0.2; 0.1 0.9]; % Define the initial state vector x(0) x0 = [10; 5]; % Define the number of time steps to simulate N = 20; % Initialize an array to store the state vectors X = zeros(2,N+1); % Assign the initial state vector to the first column of X X(:,1) = x0; % Loop over the time steps for n = 1:N % Update the state vector using the matrix A X(:,n+1) = A*X(:,n); end % Plot the state variables over time plot(0:N,X(1,:),'b',0:N,X(2,:),'r') xlabel('Time step') ylabel('State variable') legend('x_1','x_2')
The plot shows how the state variables x_1 and x_2 change over time according to the matrix A. We can see that both state variables converge to a steady value as n increases. This steady value is called an equilibrium point or a fixed point of the system. It is a special state where x(n+1) = x(n), meaning that the system does not change anymore. We can find the equilibrium point by solving the equation Ax = x, which is equivalent to finding the eigenvector of A corresponding to the eigenvalue 1. In MATLAB, we can use the command eig to find the eigenvalues and eigenvectors of a matrix:
% Find the eigenvalues and eigenvectors of A [V,D] = eig(A); % Display the eigenvalues disp(D) % Display the eigenvectors disp(V)
The output shows that the matrix A has two eigenvalues: 0.7 and 1. The eigenvector corresponding to the eigenvalue 1 is [0.8944; 0.4472]. This means that the equilibrium point of the system is x* = [0.8944; 0.4472]. We can verify this by plugging it into the equation Ax = x:
% Define the equilibrium point xstar = [0.8944; 0.4472]; % Check if Ax = x disp(A*xstar - xstar)
The output shows that the difference between Ax and x is very small, indicating that x* is indeed an equilibrium point. We can also see that the equilibrium point is proportional to the eigenvector, meaning that any multiple of the eigenvector is also an equilibrium point. For example, x* = [8.944; 4.472] is also an equilibrium point.
Nonlinear Discrete Dynamical Systems
A nonlinear discrete dynamical system is one that cannot be written in the form x(n+1) = Ax(n), but rather has some nonlinear function that relates the state variables at different time steps. For example, a nonlinear discrete dynamical system can be written in the form x(n+1) = f(x(n)), where f is a nonlinear function. A nonlinear discrete dynamical system can be implemented in MATLAB using a similar approach as before, but with a different function for updating the state vector. For example, consider the following nonlinear discrete dynamical system:
x(n+1) = 4x(n)(1-x(n))
This system has only one state variable, x(n), which represents the population size of a species in a logistic growth model. The parameter 4 represents the growth rate of the population, and the term (1-x(n)) represents the carrying capacity of the environment. The function f(x) = 4x(1-x) is called a logistic map, and it exhibits complex behavior for different values of x(0). We can implement this system in MATLAB using the following commands:
% Define the function f f = @(x) 4*x.*(1-x); % Define the initial state x(0) x0 = 0.2; % Define the number of time steps to simulate N = 50; % Initialize an array to store the state values X = zeros(1,N+1); % Assign the initial state to the first element of X X(1) = x0; % Loop over the time steps for n = 1:N % Update the state using the function f X(n+1) = f(X(n)); end % Plot the state over time plot(0:N,X,'b') xlabel('Time step') ylabel('State variable') legend('x')
e>f^k for a given x and k:
% Define the function f f = @(x) 4*x.*(1-x); % Define the value of x x = 0.2; % Define the value of k k = 2; % Initialize a variable to store f^k(x) y = x; % Loop over k times for i = 1:k % Apply the function f to y y = f(y); end % Display the value of f^k(x) disp(y)
The output shows that f^2(0.2) = 0.4096. We can check if this is a periodic point by plugging it into the equation f^2(x) = x:
% Check if f^2(x) = x disp(f(f(y)) - y)
The output shows that the difference between f^2(x) and x is very small, indicating that y is indeed a periodic point. We can also see that the periodic point is a fixed point of the function f^2, meaning that f^2(y) = y. This implies that any fixed point of f^k is also a periodic point of f with period k. For example, x* = 0.75 is a fixed point of f, and also a periodic point of f with period 1.
Complex Iterative Maps
A complex iterative map is a nonlinear discrete dynamical system that involves complex numbers. A complex number is a number that can be written in the form a + bi, where a and b are real numbers, and i is the imaginary unit, defined by i^2 = -1. Complex numbers can be represented as points in a plane, called the complex plane, where the horizontal axis is the real part and the vertical axis is the imaginary part. A complex iterative map can be written in the form z(n+1) = f(z(n)), where z(n) is a complex number at time step n, and f is a complex function. A complex iterative map can be implemented in MATLAB using the same approach as before, but with complex numbers and functions. For example, consider the following complex iterative map:
z(n+1) = z(n)^2 + c
e>c and z(0). We can implement this map in MATLAB using the following commands:
% Define the function f f = @(z,c) z.^2 + c; % Define the parameter c c = 0.25 + 0.4i; % Define the initial state z(0) z0 = 0; % Define the number of time steps to simulate N = 50; % Initialize an array to store the state values Z = zeros(1,N+1); % Assign the initial state to the first element of Z Z(1) = z0; % Loop over the time steps for n = 1:N % Update the state using the function f Z(n+1) = f(Z(n),c); end % Plot the state in the complex plane plot(Z,'b') xlabel('Real part') ylabel('Imaginary part') legend('z')
The plot shows how the state variable z moves in the complex plane according to the function f. We can see that the state variable spirals outward as n increases. This means that the system does not have any equilibrium points or periodic points, but rather diverges to infinity. This behavior is called chaos, and it occurs when the system is sensitive to initial conditions, meaning that a small change in z(0) or c can lead to a large change in z(n). We can illustrate this by changing c slightly and plotting the new trajectory:
% Define a new parameter c c = 0.25 + 0.41i; % Update the state using the new parameter c for n = 1:N Z(n+1) = f(Z(n),c); end % Plot the new trajectory in red hold on plot(Z,'r') hold off legend('z','z with new c')
The plot shows that the new trajectory is very different from the original one, even though we only changed c by a small amount. This shows that the system is chaotic and unpredictable.
Electromagnetic Waves and Optical Resonators
An electromagnetic wave is a wave of electric and magnetic fields that propagates through space. An optical resonator is a device that traps electromagnetic waves inside a cavity, creating standing waves that resonate at certain frequencies. An optical resonator can be modeled by a nonlinear discrete dynamical system that describes how the electric field inside the cavity changes over time. For example, consider the following optical resonator model:
E(n+1) = gE(n)exp(i(phi(n)+theta))
e>exp(i(x)) represents the complex exponential function, which converts a real number x into a complex number with magnitude 1 and phase x. This model can be implemented in MATLAB using the following commands:
% Define the function phi phi = @(E) abs(E).^2; % Define the parameters g and theta g = 1.01; theta = 0.1; % Define the initial state E(0) E0 = 1; % Define the number of time steps to simulate N = 100; % Initialize an array to store the state values E = zeros(1,N+1); % Assign the initial state to the first element of E E(1) = E0; % Loop over the time steps for n = 1:N % Update the state using the model equation E(n+1) = g*E(n)*exp(1i*(phi(E(n))+theta)); end % Plot the magnitude and phase of the state over time subplot(2,1,1) plot(0:N,abs(E),'b') xlabel('Time step') ylabel('Magnitude') legend('E') subplot(2,1,2) plot(0:N,angle(E),'r') xlabel('Time step') ylabel('Phase') legend('arg(E)')
The plots show how the magnitude and phase of the electric field inside the cavity change over time according to the model equation. We can see that the magnitude oscillates between two values as n increases. These two values are called bistable states of the system. They are special states where E(n+1) = E(n), meaning that the system does not change anymore. We can find the bistable states by solving the equation gEexp(i(phi(E)+theta)) = E, which is equivalent to finding the fixed points of the function f(E) = gEexp(i(phi(E)+theta)). In MATLAB, we can use the command fzero to find a fixed point of a function near a given initial guess:
% Define the function f f = @(E) g*E*exp(1i*(phi(E)+theta)) - E; % Find a fixed point near E = 0.5 Estar1 = fzero(f,0.5); % Display the fixed point disp(Estar1)
The output shows that a fixed point near E = 0.5 is E* = 0.4999 + 0.0003i. We can check if this is a bistable state by plugging it into the equation f(E) = 0:
% Check if f(E) = 0 disp(f(Estar1))
e>0 is very small, indicating that E* is indeed a bistable state. We can also see that the bistable state has a small imaginary part, meaning that it is close to a real number. We can find another bistable state by using a different initial guess:
% Find a fixed point near E = 1.5 Estar2 = fzero(f,1.5); % Display the fixed point disp(Estar2)
The output shows that a fixed point near E = 1.5 is E* = 1.4999 - 0.0003i. We can check if this is a bistable state by plugging it into the equation f(E) = 0:
% Check if f(E) = 0 disp(f(Estar2))
The output shows that the difference between f(E) and 0 is very small, indicating that E* is also a bistable state. We can also see that this bistable state has a small imaginary part, meaning that it is close to a real number. These two bistable states are symmetric with respect to the real axis, meaning that they have the same magnitude but opposite phases.
Fractals and Multifractals
A fractal is a geometric object that has self-similarity, meaning that it looks the same at different scales. A fractal can be generated by applying an iterative map to a set of points in a plane or in a space. A multifractal is a generalization of a fractal that has different scaling properties in different regions. A multifractal can be generated by applying an iterative map to a measure or a distribution on a set of points. A multifractal can be used to model phenomena that have irregular and heterogeneous patterns, such as turbulence, clouds, earthquakes, and stock markets. In this section, we will see how to generate some examples of fractals and multifractals using MATLAB.
The Mandelbrot Set
The Mandelbrot set is one of the most famous examples of a fractal. It is defined by the following complex iterative map:
z(n+1) = z(n)^2 + c
e>z(n) diverges to infinity. The points that belong to the Mandelbrot set are colored black, and the points that do not belong to the Mandelbrot set are colored according to how many iterations it takes for z(n) to exceed a certain threshold. The Mandelbrot set can be generated in MATLAB using the following commands:
% Define the function f f = @(z,c) z.^2 + c; % Define the threshold for divergence T = 2; % Define the maximum number of iterations M = 100; % Define the resolution of the image N = 1000; % Define the range of the complex plane to plot x = linspace(-2,1,N); y = linspace(-1.5,1.5,N); % Initialize an array to store the colors C = zeros(N,N); % Loop over the x and y coordinates for i = 1:N for j = 1:N % Define the complex number c c = x(i) + 1i*y(j); % Initialize the state z(0) z = 0; % Initialize a variable to store the number of iterations k = 0; % Loop until z exceeds the threshold or reaches the maximum number of iterations while abs(z)
The image shows the Mandelbrot set in the complex plane. We can see that the Mandelbrot set has a complex and intricate shape, with many smaller copies of itself embedded in its boundary. The Mandelbrot set is an example of a fractal, because it has self-similarity at different scales. We can zoom in on any part of the Mandelbrot set and see more details and patterns that resemble the whole set. For example, we can zoom in on the region near c = -0.75 + 0.1i by changing the range of x and y:
% Define the new range of x and y x = linspace(-0.8,-0.7,N); y = linspace(0.05,0.15,N); % Repeat the same steps as before to generate and plot the image
The image shows a zoomed-in view of the Mandelbrot set near c = -0.75 + 0.1i. We can see that there are more smaller copies of the Mandelbrot set in this region, and some of them have different shapes and orientations. The Mandelbrot set is an example of a complex and beautiful mathematical object that can be generated by a simple iterative map.
The Henon Map
The Henon map is another example of a fractal that can be generated by an iterative map. It is defined by the following nonlinear discrete dynamical system:
x(n+1) = y(n) + 1 - ax(n)^2 y(n+1) = bx(n)
e>a and b. The Henon map can be implemented in MATLAB using the following commands:
% Define the parameters a and b a = 1.4; b = 0.3; % Define the initial state (x(0),y(0)) x0 = 0; y0 = 0; % Define the number of time steps to simulate N = 10000; % Initialize arrays to store the state values X = zeros(1,N+1); Y = zeros(1,N+1); % Assign the initial state to the first elements of X and Y X(1) = x0; Y(1) = y0; % Loop over the time steps for n = 1:N % Update the state using the system equations X(n+1) = Y(n) + 1 - a*X(n)^2; Y(n+1) = b*X(n); end % Plot the points in the plane plot(X,Y,'b.') xlabel('x') ylabel('y') title('Henon map')
The plot shows the points (x(n),y(n)) in the plane for a = 1.4 and b = 0.3. We can see that the points form a curved shape that looks like a butterfly. This shape is called a strange attractor, because it attracts nearby points and has a fractal structure. The strange attractor has self-similarity, meaning that it looks the same at different scales. We can zoom in on any part of the strange attractor and see more details and patterns that resemble the whole shape. For example, we can zoom in on the region near (x,y) = (0.5,-0.2) by changing the axis limits:
% Zoom in on the region near (0.5,-0.2) axis([0.4 0.6 -0.25 -0.15])
The plot shows a zoomed-in view of the strange attractor near (x,y) = (0.5,-0.2). We can see that there are more smaller copies of the strange attractor in this region, and some of them have different shapes and orientations. The Henon map is an example of a simple iterative map that can generate a complex and chaotic fractal.
The Binomial Measure
The binomial measure is an example of a multifractal that can be generated by an iterative map. It is defined by the following measure or distribution on a set of points:
mu(x) = (1/2)^n if x belongs to I_n mu(x) = 0 otherwise
e>I_1 = [0,1/2], and for n = 2, we have I_2 = [0,1/4] U [1/2,3/4]. The binomial measure can be visualized by coloring each point x in the interval [0,1] according to its weight or probability. The binomial measure can be generated in MATLAB using the following commands:
% Define the number of iterations n = 10; % Define the resolution of the image N = 1000; % Define the range of the interval to plot x = linspace(0,1,N); % Initialize an array to store the weights W = zeros(1,N); % Loop over the x coordinates for i = 1:N % Define the point x x(i) = i/N; % Initialize a variable to store the weight w = 1; % Loop over the iterations for j = 1:n % Divide the interval into two equal parts x1 = floor(x(i)*2^j)/2^j; x2 = x1 + 1/2^j; % Check if x belongs to the left or right part if x(i) >= x1 && x(i) < x2 % If x belongs to the left part, do nothing else % If x belongs to the right part, flip the weight w = 1 - w; end end % Assign the weight to the array W