clear all, close all, clc syms x1 x2 C f1 = -4*x1-x1^3-x2; f2 = x1/C; [x1steady, x2steady] = solve(f1,f2)
x1steady = 0 x2steady = 0
Stability (Eigenvalue Method)
f11 = diff(f1,x1); f12 = diff(f1,x2);
f21 = diff(f2,x1); f22 = diff(f2,x2);
A = [f11 f12; f21 f22];
A = subs(A,[x1,x2],[x1steady,x2steady])
lambda = eig(A)
ezplot(lambda(1), [0,1]); hold on
ezplot(lambda(2), [0,1]);
axis([0 1 -4 1])
A = [ -4, -1] [ 1/C, 0] lambda = 1/2/C*(-4*C+2*(4*C^2-C)^(1/2)) 1/2/C*(-4*C-2*(4*C^2-C)^(1/2))
Work with the modified system
clear all, close all syms x1 x2 f1 = x1-x1^3-x2; f2 = x1;
The following code displays the direction field f1,f2 and the nullcline f1=0.
M=10; % number of samples points x1min=-2; x1max=2; % domain specification x2min=-2; x2max=2; [X1,X2]=meshgrid(x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max); dX1=X1-X1.^3-X2; % x1-component dX2=X1; % x2-component quiver(X1,X2,dX1,dX2); % matlab routine axis([x1min x1max x2min x2max]); grid on title('Direction field (the vectors may be rescaled!)'); hold on xlabel('Current I_R'); ylabel('Voltage V_c'); ezplot(f1,[x1min x1max x2min x2max]); hold on plot([0 0],[x1min x2max],'r') % plot of the vertical nulcline
We now study the linearized system around the equilibrium (0,0), by plotting a few solution curves to the linearized system.
f11 = diff(f1,x1); f12 = diff(f1,x2); f21 = diff(f2,x1); f22 = diff(f2,x2); A = [f11 f12; f21 f22]; A = subs(A,[x1,x2],[0,0])
A = [ 1, -1] [ 1, 0]
Compute the eigenvalues and eigenvectors of the matrix A
[U,D] = eig(A)
U = [ 1/2+1/2*i*3^(1/2), 1/2-1/2*i*3^(1/2)] [ 1, 1] D = [ 1/2+1/2*i*3^(1/2), 0] [ 0, 1/2-1/2*i*3^(1/2)]
Since eigenvalues are complex (conjugate), so are the eigenvectors.
X_1 = U(:,1); X_2 = U(:,2);
lambda = [D(1,1);D(2,2)]
syms t
X_1 = X_1*exp(lambda(1)*t);
X_2 = X_2*exp(lambda(2)*t);
lambda = 1/2+1/2*i*3^(1/2) 1/2-1/2*i*3^(1/2)
We first plot the direction field for the linearized systems around (0,0)
figure M=10; % number of sample points x1min=-2; x1max=2; % domain specification x2min=-2; x2max=2; [X1,X2]=meshgrid(x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max); A = double(A) % We need A in numeric format dX1=A(1,1)*X1+A(1,2)*X2; % x1-component dX2=A(2,1)*X1+A(2,2)*X2; % x2-component quiver(X1,X2,dX1,dX2); % matlab routine axis([x1min x1max x2min x2max]); grid on hold on
A = 1 -1 1 0
Plot two distinct solutions
c1=1; c2=0; X = c1*(X_1+X_2)/2+c2*(X_1-X_2)/(2*i); ezplot(X(1),X(2), [-10,1]), hold on c1=0; c2=1; X = c1*(X_1+X_2)/2+c2*(X_1-X_2)/(2*i); ezplot(X(1),X(2), [-10,1]) title('Phase Portrait the linearized system');