1.36 Compare the effect on the step response of a standard second order system as \(\zeta \) changes

Problem: Given a standard second order system defined by the transfer function \[ G\left ( s\right ) =\frac {\omega _{n}^{2}}{s^{2}+2\zeta \omega _{n}s+\omega _{n}^{2}}\]

Change \(\zeta \) and plot the step response for each to see the effect of changing \(\zeta \) (the damping coefficient).

It is easy to solve this using the step command in Matlab, and similarly in Mathematica and Maple. But here it is solved directly from the differential equation.

The transfer function is written as \[ \frac {Y\left ( s\right ) }{U\left ( s\right ) }=\frac {\omega _{n}^{2}}{s^{2}+2\zeta \omega _{n}s+\omega _{n}^{2}}\] Where \(Y\left ( s\right ) \) and \(U\left ( s\right ) \) are the Laplace transforms of the output and input respectively.

Hence the differential equation, assuming zero initial conditions becomes\[ y^{\prime \prime }\left ( t\right ) +2\zeta \omega _{n}\ y^{\prime }\left ( t\right ) +\omega _{n}^{2}\ y\left ( t\right ) =\omega _{n}^{2}\ u\left ( t\right ) \] To solve the above in Matlab using ode45, the differential equation is converted to 2 first order ODE’s as was done before resulting in\[\begin {bmatrix} x_{1}^{\prime }\\ x_{2}^{\prime }\end {bmatrix} =\begin {bmatrix} 0 & 1\\ -\omega _{n}^{2} & -2\zeta \omega _{n}\end {bmatrix}\begin {bmatrix} x_{1}\\ x_{2}\end {bmatrix} +\begin {bmatrix} 0\\ \omega _{n}^{2}\end {bmatrix} u\left ( t\right ) \]

For a step input, \(u\left ( t\right ) =1\), we obtain

\[\begin {bmatrix} x_{1}^{\prime }\\ x_{2}^{\prime }\end {bmatrix} =\begin {bmatrix} 0 & 1\\ -\omega _{n}^{2} & -2\zeta \omega _{n}\end {bmatrix}\begin {bmatrix} x_{1}\\ x_{2}\end {bmatrix} +\begin {bmatrix} 0\\ \omega _{n}^{2}\end {bmatrix} \]

Now ODE45 can be used. In Mathematica the differential equation is solved directly using DSolve and no conversion is needed.

Mathematica

solve[z_]:=Module[{sol,eq,y,t,w=1}, 
  eq=y''[t]+2 z w y'[t]+ 
     w^2 y[t]==w^2 UnitStep[t]; 
  sol=First@NDSolve[{ 
   eq,y[0]==0,y'[0]==0},y[t],{t,0,15} 
]; 
 
Plot[y[t]/.sol,{t,0,15}, 
   PlotPoints->30, 
   PlotRange->All, 
   Frame->True, 
   FrameLabel->{{"y(t)",None},{"time (sec)", 
   "step response,Subscript[\[Omega], n]=1" 
      <>", different \[Xi]"}}, 
    ImageSize->400, 
    PlotStyle->RGBColor[2z,z/2,z/3], 
    LabelStyle->Directive[14]] 
  ]; 
 
zeta = Range[0,1.4,.2]; 
r    = Map[solve,zeta]; 
Show[r,GridLines->Automatic, 
     GridLinesStyle->Dashed]
 

Matlab

function e73 
close all; clear all; 
 
t     = 0:0.1:15; 
ic    = [0 0]'; 
zeta  = 0:.2:1.4; 
omega = 1; 
sol   = zeros(length(t),length(zeta)); 
 
for i = 1:length(zeta) 
    [t,y] = ode45(@f,t,ic,[],... 
              zeta(i),omega); 
    sol(:,i) = y(:,1); 
end 
 
plot(t,sol); 
legend('0','.2','.4',... 
 '.6','.8','1','1.2','1.4'); 
title('step response, $\omega=1$ rad/sec,'... 
 'different $\zeta$ values',... 
 'interpreter','latex',... 
 'fontsize',12); 
 
ylabel('y(t)'); 
xlabel('time sec'); 
grid on; 
set(gcf,'Position',[10,10,600,600]); 
 
function xdot=f(t,x,zeta,omega) 
xdot=[x(2) 
 -omega^2*x(1)-2*zeta*omega*x(2)]+... 
 [0 
 omega^2];