Page 187 - Numerical methods for chemical engineering

Linear ODE systems and dynamic stability 173

which clearly has a steady state at (1, 1). The Jacobian matrix
 
∂ f 1 ∂ f 1

 ∂x 1 ∂x 2  (−4x 1 + 2) (1)
J(x) =   = (4.108)
∂ f 2 ∂ f 2
  (−1) (−6x 2 + 2)
∂x 1 ∂x 2
evaluated at the steady state is
(−4(1) + 2) (1) −2 1

J(x s ) = = (4.109)
(−1) (−6(1) + 2) −1 −4
The eigenvalues of the Jacobian are

Jac = [-2 1; -1 -4];
eig(Jac)’,
ans=-3-3
Both eigenvalues have negative real parts and thus the steady state (1, 1) is stable. As the
eigenvalues are real, we do not expect the response to be oscillatory. Using the routine,
function f = stable calc f(t,x);
f = zeros(2,1);
*
*
f(1) = -2 x(1)ˆ2+2 x(1) + x(2) -1;
*
*
f(2) = -x(1) - 3 x(2)ˆ2+2 x(2) + 2;
return;
the following code (ode45 is explained in further detail below) plots the system response
for a small perturbation away from the steady state,
% set initial state with small perturbation
*
x 0 = [1;1] + 0.1 randn(2,1);
% compute response
[t traj,x traj] = ode45(‘stable calc f’,[0 2],x 0);
% make plot
figure; subplot(2,1,1); plot(t traj,x traj(:,1));
xlabel(‘t’); ylabel(‘x 1(t)’); title(‘Response to small perturbation’);
subplot(2,1,2); plot(t traj,x traj(:,2));
xlabel(‘t’); ylabel(‘x 2(t)’);

A sample response is shown in Figure 4.4.
Let us next consider another system with a steady state at x s = (1, 1), but that has a
different function f 1 ,
2 2
˙ x 1 = (x 1 − 1) + 3(x 1 − 1) − (x 2 − 1) = x + x 1 − x 2 − 1 = f 1
1
2 2 (4.110)
˙ x 2 =−(x 1 − 1) − 3(x 2 − 1) − 4(x 2 − 1) =−x 1 − 3x + 2x 2 + 2 = f 2
2
The Jacobian matrix at the steady state x s is now

(2x s1 + 1) (−1) 3 −1
J(x s ) = = (4.111)
(−1) (−6x s2 + 2) −1 −4