Math 448/548 - Examples 5.3 & 5.4. RLC Circuits (p.176-183)

Contents

Example 5.3

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))
 
 

Example 5.4

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');