Download
<< Back
OutputCode
%% Question 2
%
% We study the ordinary differential equation:
%
% $$\begin{array}{rll}
% x_1' &= (a-b) x_1- c x_2+x_1(x_3+d(1-x_3^2)) &= f_1(x) \\
% x_2' &= c x_1+(a-b)x_2+x_2(x_3+d(1-x_3^2)) &= f_2(x) \\
% x_3' &= a x_3-(x_1^2+x_2^2+x_3^2) &= f_3(x)
% \end{array}$$
%
function q2_solution
%% a) Analysis of Steady States
%
% We first study the steady state $x^*=(0,0,a)^T \in R^3$
% analytically:
%
% We first calculate the jacobian ($A$) at $x^*$:
%
% $$A = \left(\begin{array}{ccc}
% \frac{\partial f_1}{\partial x_1}(x^*)
% & \frac{\partial f_1}{\partial x_2}(x^*)
% & \frac{\partial f_1}{\partial x_3}(x^*) \\
% \frac{\partial f_2}{\partial x_1}(x^*)
% & \frac{\partial f_2}{\partial x_2}(x^*)
% & \frac{\partial f_2}{\partial x_3}(x^*) \\
% \frac{\partial f_3}{\partial x_1}(x^*)
% & \frac{\partial f_3}{\partial x_2}(x^*)
% & \frac{\partial f_3}{\partial x_3}(x^*) \\
% \end{array}\right) = \left(\begin{array}{ccc}
% 2a-b+d-da^2 & -c & 0 \\ c & 2a-b+d-da^2 & 0 \\ 0 & 0 & -a
% \end{array}\right)$$
%
% We now calculate the spectrum $\sigma(A)$ using
%
% $$0 =\det(A-\lambda I) = ((2a-b-d-da^2-\lambda)^2+c^2)(-\lambda-a).$$
%
% From this we get that $\sigma(A) = \{-a, 2a-b+d-da^2+ci, 2a-b+d-da^2-ci\}$.
%
% The steady state is A-stable if the maximum of the real parts of
% the eigenvalue is negative; i.e.,
%
% $$\max_{\lambda\in\sigma(A)} Re(\lambda)=\max(-a,2a-b+d-da^2)<0$$
%
% We consider for $b=3$, $c=0.25$, $d=0.2$ and all values of $a$; then,
%
% $$\max_{\lambda\in\sigma(A)} Re(\lambda)=\max(-a,2a-2.8-0.2a^2)<0$$
%
% We get that $\max_{\lambda\in\sigma(A)} Re(\lambda) < 0$
% when $a\in(0,5-\sqrt{11})$ or $a>5+\sqrt{11}$.
%
% Now we study numerically for $a=1.0,1.95,2.02$, $b=3$, $c=0.25$, $d=0.2$.
b = 3;
c = 0.25;
d = 0.2;
disp('--------------------------------------------------');
steady_state(1.0,b,c,d);
disp('--------------------------------------------------');
steady_state(1.95,b,c,d);
disp('--------------------------------------------------');
steady_state(2.02,b,c,d);
disp('--------------------------------------------------');
%% b) $\omega$-limit
omega_limit(1.0,b,c,d);
omega_limit(1.95,b,c,d);
omega_limit(2.02,b,c,d);
%% Utility function for analysis of Steady State $x^*=(0,0)^T \in R^2$
function steady_state(a,b,c,d)
fprintf('(a,b,c,d) = (%f,%f,%f,%f)\n\n', a,b,c,d);
% Jacobian at steady state:
A = [2*a-b+d*(1-a^2), -c, 0; c, 2*a-b+d*(1-a^2), 0; 0, 0, -a];
disp('Jacobian:');
disp(A);
% Eigenvalues:
ev = eig(A);
disp('Eigenvalues:');
disp(ev);
% A-stable?:
fprintf('Maximum of real part of eigenvalues: %f\n\n', max(real(ev)));
if max(real(ev)) < 0
disp('Steady state is A-stable')
else
disp('Steady state is NOT A-stable')
end
end
%% Utility function to estimate $\omega$-limit
function omega_limit(a,b,c,d)
figure;
f = @(t,x)[ (a-b)*x(1)-c*x(2)+x(1)*(x(3)+d*(1-x(3)^2)); ...
c*x(1)+(a-b)*x(2)+x(2)*(x(3)+d*(1-x(3)^2)); ...
a*x(3)-(x(1)^2+x(2)^2+x(3)^2) ];
x0 = [.1;.1;.1];
xstar = [0;0;a];
plot3(x0(1),x0(2),x0(3),'r*');
hold on;
plot3(xstar(1),xstar(2),xstar(3),'b*');
[~,x]=ode23s(f, [0,2000], x0);
plot3(x(:,1),x(:,2),x(:,3),'r');
xlabel('x_1');
ylabel('x_2');
zlabel('x_3');
end
end
