%
%   NLIN_DEMO1.M  MATLAB demo that compares the linear and nonlinear 
%                    solutions for a simple nonlinear system
%  
%   A simple second-order nonlinear system is simulated using the ode45 routine to  
%   numerically integrate the nonlinear state equations.  The nonlinear system has  
%   two equilibrium points and these are treated as both the initial condition and
%   the reference point for a linear approximation to the nonlinear system. 
%
%   The system has no physical interpretation, since the equations were simply "made 
%   up" to present a simple example for illustration purposes.  It is interesting,
%   however, since one of the equilibrium points is in a region of stability and the 
%   other one is an unstable equilibrium point (where deviations from equilibrium 
%   continually move away from the equilibrium state). The linearized systems do
%   quite well in both cases near the linearization point, but the one that uses 
%   the unstable equilibrium as the reference state tends to deviate significantly as
%   the system moves away from the reference state.
%
%   File prepared by J. R. White, UMass-Lowell (last update: Feb. 2020)
%

      clear all,   close all,   nfig = 0;
%
%   anonymous function for ode45
      ftx = @(t,x) [-x(1) + 2*x(2) + x(2)*x(2);    
                    -x(1) - 3*x(2)              ];
%
%   let's focus on the xe = [0 0]^T equilibrium point first
%
%   numerical solution of nonlinear system (unit impulse input)
      timeo = 0;  timef = 3.5;  xo1 = [0 0]'; b = [1 0]';  xo1p = xo1 + b;
      options = odeset('RelTol',1.0e-6);
      [tn1,xn1] = ode45(ftx,[timeo,timef],xo1p,options);
%
%   linear solution (unit impulse input)
      Nt1 = 36;  tl1 = linspace(timeo,timef,Nt1);
      A1 = [-1 2;-1 -3];  B = [1 0]';  C = eye(2);  D = [0 0]';
      sys1 = ss(A1,B,C,D);  dx1 = impulse(sys1,tl1);
      xl1 = zeros(Nt1,2);
      for i = 1:2  % add reference point and deviation variable to get state vector
        xl1(:,i) = xo1(i) + dx1(:,i);
      end
%
%   compare linear and nonlinear solutions
      nfig = nfig+1;   figure(nfig)
      plot(tn1,xn1(:,1),'r-',tl1,xl1(:,1),'r.', ...
           tn1,xn1(:,2),'b--',tl1,xl1(:,2),'b.','LineWidth',2),grid
      title('NLin\_Demo1:  Impulse Response for a Nonlinear System (x_o = [0 0]^T)')
      ylabel('State Variables'), xlabel('Time (sec)')
      legend('NLin x_1(t)','Lin x_1(t)','NLin x_2(t)','Lin x_2(t)')
%
%   now let's focus on the xe = [15 -5]^T equilibrium point 
%
%   numerical solution of nonlinear system (unit impulse input)
      timeo = 0;  timef = 3.5;  xo2 = [15 -5]'; b = [1 0]';  xo2p = xo2 + b;
      options = odeset('RelTol',1.0e-6);
      [tn2,xn2] = ode45(ftx,[timeo,timef],xo2p,options);
%
%   linear solution (unit impulse input)
      Nt2 = 36;  tl2 = linspace(timeo,timef,Nt2);
      A2 = [-1 -8;-1 -3];  B = [1 0]';  C = eye(2);  D = [0 0]';
      sys2 = ss(A2,B,C,D);  dx2 = impulse(sys2,tl2);
      xl2 = zeros(Nt2,2);
      for i = 1:2   % add reference point and deviation variable to get state vector
        xl2(:,i) = xo2(i) + dx2(:,i);
      end
%
%   compare linear and nonlinear solutions
      nfig = nfig+1;   figure(nfig)
      plot(tn2,xn2(:,1),'r-',tl2,xl2(:,1),'r.', ...
           tn2,xn2(:,2),'b--',tl2,xl2(:,2),'b.','LineWidth',2),grid
      title('NLin\_Demo1:  Impulse Response for a Nonlinear System (x_o = [15 -5]^T)')
      ylabel('State Variables'), xlabel('Time (sec)')
      legend('NLin x_1(t)','Lin x_1(t)','NLin x_2(t)', ...
             'Lin x_2(t)','Location','NorthWest')
%
%  end of demo