UP
PDF letter size
PDF legal size

How to solve basic engineering and mathematics problems using Mathematica, Matlab and Maple

Nasser M. Abbasi

May 30, 2016 page compiled on May 30, 2016 at 2:41am

This is a collection of how to examples showing the use of Mathematica and Matlab to solve basic engineering and mathematics problems. Few examples are also in Maple, Ada, and Fortran.

This was started as a cheat sheet few years ago, and I continue to update it all the time.

Few of the Matlab examples require the use of toolboxs such as signal processing toolbox and the control systems toolbox (these come free as part of the student version). Most examples require only the basic system installation.

1 Control systems, Linear systems, transfer functions, state space related problems
2 Linear algebra, Linear solvers, common operations on vectors and matrices
3 signal processing, Image processing, graphics, random numbers
4 Differential, PDE solving, integration, numerical and analytical solving of equations
5 plotting how to, ellipse, 3D

Chapter 1
Control systems, Linear systems, transfer functions, state space related problems

 1.1 Creating tf and state space and different Conversion of forms
 1.2 Obtain the step response of an LTI from its transfer function
 1.3 plot the impulse and step responses of a system from its transfer function
 1.4 Obtain the response of a transfer function for an arbitrary input
 1.5 Obtain the poles and zeros of a transfer function
 1.6 Generate Bode plot of a transfer function
 1.7 How to check that state space system x ′ = Ax + Bu  is controllable?
 1.8 Obtain partial-fraction expansion
 1.9 Obtain Laplace transform for a piecewise functions
 1.10 Obtain Inverse Laplace transform of a transfer function
 1.11 Display the response to a unit step of an under, critically, and over damped system
 1.12 View steady state error of 2nd order LTI system with changing undamped natural frequency
 1.13 Show the use of the inverse Z transform
 1.14 Find the Z transform of sequence x(n)
 1.15 Sample a continuous time system
 1.16 Find closed loop transfer function from the open loop transfer function for a unity feedback
 1.17 Compute the Jordan canonical/normal form of a matrix A
 1.18 Solve the continuous-time algebraic Riccati equation
 1.19 Solve the discrete-time algebraic Riccati equation
 1.20 Display impulse response of H(z) and the impulse response of its continuous time approximation H(s)
 1.21 Find the system type given an open loop transfer function
 1.22 Find the eigenvalues and eigenvectors of a matrix
 1.23 Find the characteristic polynomial of a matrix
 1.24 Verify the Cayley-Hamilton theorem that every matrix is zero of its characteristic polynomial
 1.25 How to check for stability of a continuous system represented as a transfer function and state space
 1.26 Check continuous system stability in the Lyapunov sense
 1.27 Given a closed loop block diagram, generate the closed loop Z transform and check its stability
 1.28 Determine the state response of a system to only initial conditions in state space
 1.29 Determine the response of a system to only initial conditions in state space
 1.30 Determine the response of a system to step input with nonzero initial conditions
 1.31 Draw the root locus from the open loop transfer function
 1.32 Find eAt  where A  is a matrix
 1.33 Draw root locus for a discrete system
 1.34 Plot the response of the inverted pendulum problem using state space
 1.35 How to build and connect a closed loop control systems and show the response?
 1.36 Compare the effect on the step response of a standard second order system as ζ  changes
 1.37 Plot the dynamic response factor Rd  of a system as a function of     -ω
r = ωn   for different damping ratios
 1.38 How to find closed loop step response to a plant with a PID controller?
 1.39 How to make Nyquist plot?

1.1 Creating tf and state space and different Conversion of forms

1.1.1 Create continuous time transfer function given the poles, zeros and gain

Problem: Find the transfer function H (s)  given zeros s = − 1,s = − 2  , poles s = 0,s = − 4,s = − 6  and gain 5.



Mathematica

Clear["Global`*"]; 
num  = (s+1)(s+2); 
den  = (s)(s+4)(s+6); 
gain = 5; 
sys  = TransferFunctionModel[ 
           gain (num/den),s]
  Out[30]= TransferFunctionModel[{
       {{5*(1 + s)*(2 + s)}},
       s*(4 + s)*(6 + s)}, s]

Matlab

clear all; 
m_zeros = [-1 -2]; 
m_poles = [0 -4 -6]; 
gain  = 5; 
sys   = zpk(m_zeros,m_poles,gain); 
[num,den] = tfdata(sys,'v'); 
printsys(num,den,'s')
  num/den =
      5 s^2 + 15 s + 10
     -------------------
     s^3 + 10 s^2 + 24 s





Maple

restart; 
alias(DS=DynamicSystems): 
zeros :=[-1,-2]: 
poles :=[0,-4,-6]: 
gain  :=5: 
sys   :=DS:-TransferFunction(zeros,poles,gain): 
#print overview of the system object 
DS:-PrintSystem(sys); 
 
exports(sys); 
#To see fields in the sys object, do the following 
# tf, inputcount, outputcount, statecount, sampletime, 
# discrete, systemname, inputvariable, outputvariable, 
# statevariable, parameters, systemtype, ModulePrint 
tf:=sys[tf]; #reads the transfer function 
Matrix(1,1,{(1,1)=(5*s^2+15*s+10)/(s^3+10*s^2+24*s)}) 
numer(tf[1,1]); #get the numerator of the tf 
            5*s^2+15*s+10 
denom(tf[1,1]); #get the denominator of the tf 
             s*(s^2+10*s+24)

     [ 5s2+15s+10]
tf =  s3+10s2+24s



1.1.2 Convert transfer function to state space representation

problem 1

Problem: Find the state space representation for the continuous time system defined by the transfer function

        -----5s-----
H (s) = s2 + 4s + 25



Matematica

Clear["Global`*"]; 
sys = TransferFunctionModel[(5 s)/(s^2+4 s+25),s]; 
ss  = MinimalStateSpaceModel[StateSpaceModel[sys]]

pict

(a=Normal[ss][[1]])//MatrixForm 
(b=Normal[ss][[2]])//MatrixForm 
(c=Normal[ss][[3]])//MatrixForm 
(d=Normal[ss][[4]])//MatrixForm

pict





Matlab

clear all; 
s      = tf('s'); 
sys_tf = (5*s) / (s^2 + 4*s + 25); 
[num,den] = tfdata(sys_tf,'v'); 
[A,B,C,D] = tf2ss(num,den)

  A =
      -4   -25
       1     0
  B =
       1
       0
  C =
       5     0
  D =
       0





Maple

restart; 
alias(DS=DynamicSystems): 
sys := DS:-TransferFunction(5*s/(s^2+4*s+25)): 
sys := DS:-StateSpace(sys); 
exports(sys); #to see what fields there are 
a := sys:-a;

[        ]
   0  1
  − 25 − 4

b:=sys:-b;

[  ]
 0
 1

c:=sys:-c;

[  ]
 05

d:=sys:-d;

[ ]
0



problem 2

Problem: Given the continuous time S transfer function defined by

H (s) =  --1-+-s---
         s2 + s + 1

convert to state space and back to transfer function.



Mathematica

Clear["Global`*"]; 
sys = TransferFunctionModel[(s + 1)/(s^2 + 1 s + 1), s]; 
ss = MinimalStateSpaceModel[StateSpaceModel[sys]]

(        )
   0 1 0
( − 1 − 11 )

   1 1 0

TransferFunctionModel[ss]

---s +-1--
s2 + s + 1



Matlab

clear all; 
s   = tf('s'); 
sys = ( s+1 )/(s^2  + s  +  1); 
[num,den]=tfdata(sys,'v'); 
%convert to state space 
[A,B,C,D] = tf2ss(num,den)

  A =
      -1    -1
       1     0
  
  B =
       1
       0
  
  C =
       1     1
  
  D =
       0

%convert from ss to tf 
[num,den] = ss2tf(A,B,C,D); 
sys=tf(num,den)

  sys =
       s + 1
    -----------
    s^2 + s + 1



1.1.3 Create a state space representation from A,B,C,D and find the step response

Problem: Find the state space representation and the step response of the continuous time system defined by the Matrices A,B,C,D as shown

  A =
       0     1     0     0
       0     0     1     0
       0     0     0     1
    -100   -80   -32    -8
  
  B =
       0
       0
       5
      60
  
  C =
       1     0     0     0
  
  D =
       0


Mathematica

Clear["Global`*"]; 
 
a = {{0,1,0,0}, 
     {0,0,1,0}, 
     {0,0,0,1}, 
     {-100,-80,-32,-8} 
    }; 
 
b = {{0} ,{0}, {5},{60}}; 
c = {{1,0,0,0}}; 
d = {{0}}; 
 
sys = StateSpaceModel[{a,b,c,d}]; 
y   = OutputResponse[sys,UnitStep[t],{t,0,10}]; 
 
p = Plot[Evaluate@y, {t, 0, 10}, 
         PlotRange -> All, 
         GridLines -> Automatic, 
         GridLinesStyle -> LightGray]

pict



Matlab

clear all; 
A = [0 1 0 0; 
     0 0 1 0; 
     0 0 0 1; 
   -100 -80 -32 -8]; 
B = [0;0;5;60]; 
C = [1 0 0 0]; 
D = [0]; 
 
sys = ss(A,B,C,D); 
 
step(sys); 
grid; 
title('unit step response');

pict



Maple

restart: 
alias(DS=DynamicSystems): 
alias(size=LinearAlgebra:-Dimensions); 
a:=Matrix([[0,1,0,0],[0,0,1,0], 
           [0,0,0,1],[-100,-90,-32,-8]]): 
b:=Matrix([[0],[0],[5],[6]]): 
c:=Matrix([[1,0,0,0]]): 
d:=Matrix([[1]]): 
sys:=DS:-StateSpace(a,b,c,d); 
 
DS:-ResponsePlot(sys, DS:-Step(), 
  duration=10,color=red,legend="step");

pict



1.1.4 Convert continuous time to discrete time transfer function, gain and phase margins

Problem: Compute the gain and phase margins of the open-loop discrete linear time system sampled from the continuous time S transfer function defined by

           1 + s
H (s) =  -2--------
         s + s + 1

Use sampling period of 0.1 seconds.



Mathematica

Clear["Global`*"]; 
tf = TransferFunctionModel[(s+1)/(s^2+1 s+1),s];

pict

ds = ToDiscreteTimeModel[tf, 0.1, z, 
                Method -> "ZeroOrderHold"]

pict

{gm, pm} = GainPhaseMargins[ds];

  gm
     Out[28]={{31.4159,19.9833}}
  pm
     Out[29]={{1.41451,1.83932},{0.,Pi}}

max = 100; 
yn  = OutputResponse[ds, 
          Table[UnitStep[n],{n, 0, max}]]; 
 
ListPlot[yn, 
      Joined -> True, 
      InterpolationOrder -> 0, 
      PlotRange -> {{0, max}, {0, 1.4}}, 
      Frame -> True, 
      FrameLabel ->{{"y[n]", None}, 
                   {"n", "unit step response"}}, 
      ImageSize -> 400, 
      AspectRatio -> 1,BaseStyle -> 14]

pict



Matlab

clear all; close all; 
Ts   = 0.1; %sec, sample period 
s    = tf('s'); 
sys  = ( s + 1 ) / (  s^2  + s  +  1); 
%convert s to z domain 
sysd = c2d(sys,Ts,'zoh'); 
tf(sysd)

     0.09984 z - 0.09033
    ----------------------
    z^2 - 1.895 z + 0.9048
  
  Sample time: 0.1 seconds
  Discrete-time transfer function.
  
  Gm =
     19.9833
  
  Pm =
    105.3851
  
  Wcg =
     31.4159
  
  Wcp =
      1.4145

step(sysd,10); 
[Gm,Pm,Wcg,Wcp] = margin(sysd)

pict



1.1.5 Convert differential equation to transfer functions and to state space

Problem: Obtain the transfer and state space representation for the differential equation

 d2y     dy
3--2-+ 2 ---+ y (t) = u(t)
 dt      dt



Mathematica

Clear[y, u, t]; 
eq =  3  y''[t] + 2 y'[t] +  y[t] == u[t]; 
ss = StateSpaceModel[ {eq}, {{y'[t], 0}, 
        {y[t], 0}}, {{u[t], 0}}, {y[t]}, t]; 
ss = MinimalStateSpaceModel[ss]

pict

tf = Simplify[TransferFunctionModel[ss]]

pict



Matlab

clear all; 
 
syms y(t) Y 
eq  = 3*diff(y,2)+2*diff(y)+y; 
lap = laplace(eq); 
lap = subs(lap,'y(0)',0); 
lap = subs(lap,'D(y)(0)',0); 
lap = subs(lap,'laplace(y(t),t,s)',Y); 
H   = solve(lap-1,Y); 
pretty(H)

          1
    --------------
       2
    3 s  + 2 s + 1

[num,den] = numden(H); 
num = sym2poly(num); 
den = sym2poly(den); 
[A,B,C,D] = tf2ss(num,den)

  A =
     -0.6667   -0.3333
      1.0000         0
  B =
       1
       0
  C =
           0    0.3333
  D =
       0



Maple

restart; 
alias(DS=DynamicSystems): 
 
ode:=3*diff(y(t),t$2)+2*diff(y(t),t)+y(t)= 
                             Heaviside(t): 
 
sys:=DS:-DiffEquation(ode, 
          'outputvariable'=[y(t)], 
          'inputvariable'=[Heaviside(t)]): 
 
sys:=DS:-TransferFunction(sys): 
sys:-tf(1,1);

     1
--2---------
3s +  2s + 1

sys:=DS:-StateSpace(sys): 
#show the state space matrices 
{sys:-a,sys:-b,sys:-c,sys:-d};

{ [         ] [ ]          }
     0   1     0  [    ] [ ]
   − 1∕3− 2∕3 , 1 , 1∕30 , 0



1.1.6 Convert from continuous transfer function to discrete time transfer function

Problem: Convert

H (s) = ------1-------
        s2 + 10s + 10

To Z  domain transfer function using sampling period of 0.01  seconds and using the zero order hold method.



Mathematica

Clear["Global`*"] 
sys = TransferFunctionModel[ 
              1/(s^2 +10 s+20),s]; 
 
sysd= ToDiscreteTimeModel[sys,0.01,z, 
        Method->"ZeroOrderHold"]

pict



Matlab

s    = tf('s'); 
sys  = 1/(s^2 + 10*s + 20); 
T    = 0.01; %seconds 
sysz = c2d(sys,T,'zoh')

  sysz =
  
    4.837e-05 z + 4.678e-05
    -----------------------
    z^2 - 1.903 z + 0.9048
  
  Sample time: 0.01 seconds
  Discrete-time transfer function.



1.1.7 Convert a Laplace transfer function to an ordinary differential equation

Problem: Give a continuous time transfer function, show how to convert it to an ordinary differential equation. This method works for non-delay systems. The transfer function must be ratio of polynomials. For additional methods see this question at stackexchange

Mathematica

tfToDiff[tf_, s_, t_, y_, u_] := Module[{rhs, lhs, n, m}, 
  rhs = Numerator[tf]; 
  lhs = Denominator[tf]; 
  rhs = rhs /. m_. s^n_. -> m D[u[t], {t, n}]; 
  lhs = lhs /. m_. s^n_. -> m D[y[t], {t, n}]; 
  lhs == rhs 
  ] 
 
tf = (5s)/(s^2+4s+25); 
tf = C0 s/(R0 C0 s + 1); 
eq=tfToDiff[tf, s, t, y, u]

 Out = 25 + 4 y'[t]+y''[t]==5 u'[t] 

1.2 Obtain the step response of an LTI from its transfer function

Problem: Find the unit step response for the continuous time system defined by the transfer function

        -----25-----
H (s) = s2 + 4s + 25



Mathematica

Clear["Global`*"]; 
 
tf=TransferFunctionModel[25/(s^2+4s+25),s]; 
y = OutputResponse[tf, UnitStep[t], t]; 
 
Plot[Evaluate@y, {t, 0, 4}, 
 PlotRange -> {{0, 4}, {0, 1.4}}, 
 Frame -> True, 
 FrameLabel -> {{"y(t)", None}, 
                {t, "step response"}}, 
 GridLines -> Automatic, 
 GridLinesStyle -> Dashed, 
 ImageSize -> {300, 300}, 
 PlotStyle -> Red, 
 AspectRatio -> 1]

pict



Matlab

clear all; 
s     = tf('s'); 
sys   = 25/(s^2+4*s+25); 
[y,t] = step(sys,4); 
 
plot(t,y,'r'); 
xlim([0 4]); 
ylim([0 1.4]); 
title('step response'); 
xlabel('t'); 
ylabel('y(t)'); 
grid 
set(gcf,'Position',[10,10,310,310]);

pict



Maple

restart: 
with(DynamicSystems): 
sys := TransferFunction(25/(s^2+4*s+25)): 
ResponsePlot(sys, Step(),duration=4);

pict



1.3 plot the impulse and step responses of a system from its transfer function

Problem: Find the impulse and step responses for the continuous time system defined by the transfer function

              1
H (s) = -2-----------
        s  + 0.2s + 1

and display these on the same plot up to some t  value.

Side note: It was easier to see the analytical form of the responses in Mathematica and Maple so it is given below the plot.



Mathematica

Clear["Global`*"]; 
SetDirectory[NotebookDirectory[]] 
sys = TransferFunctionModel 
        [1/(s^2 + 2/10 s + 1), s]; 
 
yStep = Assuming[t > 0, 
  Simplify@OutputResponse[ 
             sys, UnitStep[t], t]]

  {1 - E^(-t/10) Cos[(3 Sqrt[11] t)/10]
   - (
  E^(-t/10) Sin[(3 Sqrt[11]t)/10])/
                         (3 Sqrt[11])}

yImpulse = Simplify@ 
     OutputResponse[sys,DiracDelta[t],t]

  {(10 E^(-t/10) HeavisideTheta[t]
    Sin[(3 Sqrt[11] t)/10])/(3 Sqrt[11])}

\begin{X301} 
p = Grid[{ 
  {Plot[Evaluate@{yStep, yImpulse}, 
   {t, 0, 50}, 
  PlotRange -> {{0, 50}, {-0.8, 2.0}}, 
  Frame -> True, 
  FrameLabel -> {{"y(t)", None}, 
     {t, "step and impulse reponse"}}, 
  GridLines -> Automatic, 
  GridLinesStyle -> Dashed, 
  ImageSize -> {300, 300}, 
  PlotStyle -> {Red, Black}, 
  AspectRatio -> 1] 
  }, 
  {Row[{"step response=", Re@yStep}]}, 
  {Row[{"impulse response=",Re@yImpulse}]} 
  }, Alignment -> Left, Frame -> True]

pict



Matlab

clear all; close all; 
t   = 0:0.05:50; 
s   = tf('s'); 
sys = 1/(s^2+0.2*s+1); 
y   = step(sys,t); 
plot(t,y,'-r') 
hold on 
y  = impulse(sys,t); 
plot(t,y,'-k') 
title('Step and Impulse responses'); 
xlabel('t'); 
ylabel('y(t)'); 
xlim([0 50]); 
ylim([-0.8 2]); 
legend('step','impulse'); 
grid on; 
set(gcf,'Position',[10,10,310,310]);

pict



Maple Using Maple DynamicSystems

restart: 
alias(DS=DynamicSystems): 
sys:=DS:-TransferFunction(1/(s^2+0.2*s+1)): 
p1:=DS:-ResponsePlot(sys, DS:-Step(), 
   duration=50,color=red,legend="step"): 
p2:=DS:-ImpulseResponsePlot(sys, 
                   50, 
                   legend="impulse"): 
plots:-display([p1,p2],axes=boxed, 
   title=`step and impulse reponses`);

pict

Using Laplace transform method:

with(inttrans): 
with(plottools): 
H:=1/(s^2+0.2*s+1); 
input:=laplace(Heaviside(t),t,s): 
yStep:=invlaplace(input*H,s,t); 
input:=laplace(Dirac(t),t,s): 
yImpulse:=invlaplace(input*H,s,t); 
plot([yStep,yImpulse],t=0.001..50, 
     title=`step response`, 
     labels=[`time`,`y(t)`], 
     color=[red,black], 
     legend=["step","impulse"]);

pict



1.4 Obtain the response of a transfer function for an arbitrary input

Problem: Find the response for the continuous time system defined by the transfer function

H (s) = ------1------
        s2 + 0.2s + 1

when the input is given by

u (t) = sin(t)

and display the response and input on the same plot.

Side note: It was easier to see the analytical form of the responses in Mathematica and Maple so it is given below the plot.



Mathematica

Clear["Global`*"]; 
SetDirectory[NotebookDirectory[]] 
sys=TransferFunctionModel[1/(s^2+2/10 s+1),s]; 
u = Sin[t]; 
y = OutputResponse[sys, u, t]; 
p = Grid[{ 
   {Plot[Evaluate@{y, u}, {t, 0, 50}, 
     PlotRange -> {{0, 50}, {-5, 5}}, 
     Frame -> True, 
     FrameLabel -> {{"y(t)", None}, 
           {"t", "input and its reponse"}}, 
     GridLines -> Automatic, 
     GridLinesStyle -> Dashed, 
     ImageSize -> {300, 300}, 
     PlotStyle -> {Red, Black}, 
     AspectRatio -> 1] 
    } 
   }, Alignment -> Left, Frame -> True]

pict



Matlab

clear all; 
close all; 
t   = 0:0.05:50; 
u   = sin(t); 
s   = tf('s'); 
sys = 1/(s^2+0.2*s+1); 
[num,den] = tfdata(sys,'v'); 
y = lsim(num,den,u,t); 
plot(t,u,'k',t,y,'r'); 
title('input and its response'); 
xlabel('t'); ylabel('y(t)'); 
xlim([0 50]); 
ylim([-5 5]); 
grid on; 
set(gcf,'Position',[10,10,310,310]);

pict



Maple

restart; 
with(inttrans):H:=1/(s^2+0.2*s+1); 
input:=sin(t); 
inputS:=laplace(input,t,s): 
y:=invlaplace(inputS*H,s,t); 
plot([y,input],t=0..25, 
 title=`system response`, 
 labels=[`time`,`y(t)`], 
 color=[black,red], 
 legend=["response","input"]);

pict

Using DynamicSystem package

restart: 
alias(DS=DynamicSystems): 
sys :=DS:-TransferFunction(1/(s^2+0.2*s+1)): 
p1:=DS:-ResponsePlot(sys, sin(t), 
 duration=25,color=black,legend="response"): 
p2:=plot(sin(t),t=0..25,color=red, 
   legend="input",size=[400,"default"]): 
plots:-display([p2,p1],axes=boxed, 
  title="step and impulse reponses", 
  legendstyle= [location=right]);

pict



1.5 Obtain the poles and zeros of a transfer function

Problem: Find the zeros, poles, and gain for the continuous time system defined by the transfer function

H (s) = -----25-----
        s2 + 4s + 25



Mathematica

Clear["Global`*"]; 
sys=TransferFunctionModel[25/(s^2+4 s+25),s]; 
TransferFunctionZeros[sys]
     {{{}}}
TransferFunctionPoles[sys]//N
     {{{-2.-4.58258 I,-2.+4.58258 I}}}

Matlab

clear all; 
s = tf('s'); 
sys = 25 / (s^2 + 4*s + 25); 
[z,p,k] =zpkdata(sys,'v')
  
  z =
     Empty matrix: 0-by-1
  
  p =
    -2.0000 + 4.5826i
    -2.0000 - 4.5826i



Maple

restart; 
alias(DS=DynamicSystems): 
sys := DS:-TransferFunction(25/(s^2+4*s+25)): 
r   := DS:-ZeroPolePlot(sys,output=data): 
zeros  := r[1]; 
poles  := r[2];
  zeros := []
  poles := [-2.000000000-4.582575695*I,
            -2.+4.582575695*I]



1.6 Generate Bode plot of a transfer function

Problem: Generate a Bode plot for the continuous time system defined by the transfer function

             5s
H (s) = -2----------
        s  + 4s + 25



Mathematica

Clear["Global`*"]; 
 
tf=TransferFunctionModel[(5 s)/(s^2+4s+25),s]; 
 
BodePlot[tf, GridLines -> Automatic, 
 ImageSize -> 300, 
 FrameLabel -> {{{"magnitude (db)", None}, 
       {None, "Bode plot"}}, 
       {{"phase(deg)", None}, 
        {"Frequency (rad/sec)", None}}}, 
 ScalingFunctions -> {{"Log10", "dB"}, 
                      {"Log10", "Degree"}}, 
 PlotRange -> {{{0.1, 100}, Automatic}, 
               {{0.1, 100}, Automatic}} 
 ]

pict



Matlab

clear all; 
s   = tf('s'); 
sys = 5*s / (s^2 + 4*s + 25); 
bode(sys); 
grid; 
set(gcf,'Position',[10,10,400,400]);

pict



Maple

restart: 
alias(DS=DynamicSystems): 
sys:=DS:-TransferFunction(5*s/(s^2+4*s+25)): 
DS:-BodePlot(sys,output=dualaxis, 
             range=0.1..100);

or can plot the the two bode figures on top of each others as follows

DS:-BodePlot(sys,size=[400,300], 
 output=verticalplot,range=0.1..100);

pict



1.7 How to check that state space system  ′
x = Ax  + Bu  is controllable?

A system described by

pict

Is controllable if for any initial state x0   and any final state xf  there exist an input u  which moves the system from x0   to xf  in finite time. Only the matrix A  and B  are needed to decide on controllability. If the rank of

[B AB  A2B  ... An− 1B ]

is n  which is the number of states, then the system is controllable. Given the matrix

     (               )
        0  1   0   0
     |  0  0  − 1  0 |
A =  |(  0  0   0   1 |)

        0  0   5   0

And

    (      )
        0
B = ||   1  ||
    (   0  )
       − 2



Mathematica

A0 = {{0, 1, 0, 0}, 
      {0, 0, -1, 0}, 
      {0, 0, 0, 1}, 
      {0, 0, 5, 0}}; 
B0 = {{0}, {1}, {0}, {-2}}; 
sys = StateSpaceModel[{A0, B0}]; 
m = ControllabilityMatrix[sys]

(               )
   0  1   0   2
||  1  0   2   0 ||
(  0 − 2  0 − 10)
  − 2 0 − 10  0

ControllableModelQ[sys]

True

MatrixRank[m]

4



Matlab

A0 = [0 1 0 0; 
      0 0 -1 0; 
      0 0 0 1; 
      0 0 5 0]; 
B0 = [0 1 0 -2]'; 
sys = ss(A0,B0,[],[]); 
m   = ctrb(sys)

  m =
       0     1     0     2
       1     0     2     0
       0    -2     0   -10
      -2     0   -10     0

rank(m)

4



Maple

restart: 
alias(DS=DynamicSystems): 
A:=Matrix( [ [0,1,0,0], 
             [0,0,-1,0], 
             [0,0,0,1], 
             [0,0,5,0] 
           ] 
          ); 
B:=Matrix([[0],[1],[0],[-2]]); 
sys:=DS:-StateSpace(A,B); 
m:=DS:-ControllabilityMatrix(sys);

⌊              ⌋
  0  1   0   2
|| 1  0   2   0 ||
|              |
|⌈ 0 − 2  0 − 10|⌉

 − 2 0 − 10  0

DS:-Controllable(sys,method=rank);

true

DS:-Controllable(sys,method=staircase);

true

LinearAlgebra:-Rank(m);

4



1.8 Obtain partial-fraction expansion

Problem: Given the continuous time S transfer function defined by

        s4 +-8s3-+-16s2-+-9s-+-9
H (s) =    s3 + 6s2 + 11s + 6

obtain the partial-fractions decomposition.

Comment: Mathematica result is easier to see visually since the partial-fraction decomposition returned in a symbolic form.



Mathematica

Remove["Global`*"]; 
 
expr = (s^4+8 s^3+16 s^2+9 s+6)/ 
            (s^3+6 s^2+11 s+6); 
Apart[expr]

  2 + s +3/(1+s) -4/(2+s) -6/(3+s)



Matlab

clear all; 
s=tf('s'); 
tf_sys = (s^4+8*s^3+16*s^2+9*s+6)/... 
          (s^3+6*s^2+11*s+6); 
[num,den] = tfdata(tf_sys,'v'); 
[r,p,k] = residue(num,den)

  r =
     -6.0000
     -4.0000
      3.0000
  
  p =
     -3.0000
     -2.0000
     -1.0000
  
  k =
       1     2



Maple

p:=(s^4+8*s^3+16*s^2+9*s+9)/(s^3+6*s^2+11*s+6); 
p0:=convert(p,parfrac);

s + 2 − --7--−  --9----+ ---9---
        s + 2   2s + 6   2s + 2

[op(p0)];

[s, 2,--7--,− ---9--, --9---]
     s + 2   2s + 6  2s + 2



1.9 Obtain Laplace transform for a piecewise functions

Problem: Obtain the Laplace transform for the function defined in the following figure.

pict

Function f(t) to obtain its Laplace transform

Comment: Mathematica solution was easier than Matlab’s. In Matlab the definition of the Laplace transform is applied to each piece separately and the result added. Not finding the piecewise maple function to access from inside MATLAB did not help.



Mathematica

Remove["Global`*"]; 
f[t_] := Piecewise[{{0,t<0}, 
                    {t,t>= 0 && t< T}, 
                    {T, t>T}}] 
Simplify[LaplaceTransform[f[t],t,s] ,{T>0}]

  Out[]= (1 - E^((-s)*T))/s^2



Matlab

clear all; 
syms T t s; 
syms s positive; 
I1 = int(t*exp(-s*t),t,0,T); 
I2 = int(T*exp(-s*t),t,T,Inf); 
result = simple(I1+I2); 
pretty(result)

    1 - exp(-T s)
    -------------
          2
         s



Maple

With Maple, had to use Heaviside else Laplace will not obtain the transform of a piecewise function.

restart; 
assume(T>0): 
interface(showassumed=0): 
f:=t->piecewise(t<0,0,t>=0 and t<T,t,t>T,T): 
r:=convert(f(t),Heaviside): 
r:=inttrans[laplace](r,t,s);

    − sT
1 −-e----
   s2



1.10 Obtain Inverse Laplace transform of a transfer function

Problem: Obtain the inverse Laplace transform for the function

           4    3     2
H (s) = --s-+-5s--+-6s--+-9s-+-30-
        s4 + 6s3 + 21s2 + 46s + 30

Mathematica

Remove["Global`*"]; 
f = (s^4+5 s^3+6 s^2+9 s+30)/(s^4+6 s^3+21 s^2+46 s+30); 
InverseLaplaceTransform[f,s,t]; 
Expand[FullSimplify[%]]
       (           )
          1     i     (− 1− 3i)t(            6it               )   3e− 3t   23e −t
δ(t) +  ---- + ----  e        (73 + 326i)e  + (− 326 − 73i)  − ----- + ------
        234    234                                              26       18

Matlab

clear all; 
syms s t 
f = (s^4+5*s^3+6*s^2+9*s+30)/(s^4+6*s^3+21*s^2+46*s+30); 
pretty(f)

      4      3      2
     s  + 5 s  + 6 s  + 9 s + 30
    -----------------------------
     4      3       2
    s  + 6 s  + 21 s  + 46 s + 30
pretty(ilaplace(f))

                                                      /            399 sin(3 t) \
                                          253 exp(-t) | cos(3 t) + ------------ |
    23 exp(-t)   3 exp(-3 t)                          \                253      /
    ---------- - ----------- + dirac(t) - ---------------------------------------
        18           26                                     117

maple

restart; 
interface(showassumed=0): 
p:=(s^4+5*s^3+6*s^2+9*s+30)/(s^4+6*s^3+21*s^2+46*s+30); 
r:=inttrans[invlaplace](p,s,t);
              −3t                                      −t
Dirac (t) − 3e---+  (− 506-cos-(3t)-−-798-sin-(3t) +-299-)e-
             26                      234

1.11 Display the response to a unit step of an under, critically, and over damped system

Problem: Obtain unit step response of the second order system given by the transfer function

                 2
H  (s) = ------ω-n-------
         s2 + 2 ξωns + ω2n

in order to illustrate the response when the system is over, under, and critically damped. use ωn =  1  and change ξ  over a range of values that extends from under damped to over damped.



Mathematica

Clear["Global`*"]; 
Needs["PlotLegends`"] 
 
sys=TransferFunctionModel[w^2/(s^2+2*z*w*s+w^2),s]; 
zValues = Range[.2,1.2,.2]; 
fun = OutputResponse[sys/.{w->1,z->#}, UnitStep[t], 
       {t,0,20}]&/@zValues; 
 
Plot[Evaluate@Flatten@Table[fun[[i]], 
                         {i,1,Length[fun]}], 
   {t,0,20}, 
   Frame->True, 
   FrameLabel->{{"y(t)",None}, 
 {"t","step response for different \[Xi] values"}}, 
   PlotRange->{{0,20},{0,1.6}}, 
   GridLines->Automatic, 
   GridLinesStyle->Dashed, 
   PlotLegend->zValues, 
   LegendPosition->{0.76,-0.5},LegendSize -> 1, 
   LegendLabel -> "\[Xi] values", 
   ImageSize -> 450, 
   LabelStyle -> 14,AspectRatio -> 1 
]

pict



Matlab

clear all; close all; 
 
wn = 1; 
z  = .2:.2:1.2; 
t  = 0:0.05:20; 
y  = zeros(length(t),length(z)); 
legendStr = cell(length(z),1); 
 
for i = 1:length(z) 
    [num,den] = ord2(wn,z(i)); 
    num = num*wn^2; 
    [y(:,i),~,~] = step(num,den,t); 
    legendStr(i) = {sprintf('%1.1f',z(i))}; 
end 
plot(t,y); 
legend(legendStr); 
 
title(sprintf( 
 '2nd order system step response with changing %s', 
 '\zeta')); 
 
xlabel('time (sec)'); 
ylabel('amplitude'); 
grid on; 
set(gcf,'Position',[10,10,400,400]);

pict



Maple

restart; 
alias(DS=DynamicSystems): 
H   := (w,zeta)->w^2/(s^2+2*zeta*w*s+w^2): 
sys := (w,zeta)->DS:-TransferFunction(H(w,zeta)): 
zetaValues := [seq(i,i=0.1..1.2,.2)]: 
sol:=map(z->DS:-Simulate(sys(1,z), 
                 Heaviside(t)),zetaValues): 
c :=ColorTools:-GetPalette("Niagara"): 
 
plots:-display(seq(plots[odeplot] 
       (sol[i],t=0..20,color=c[i], 
 
 legend=typeset(zeta,"=",zetaValues[i]), 
 legendstyle=[location=right]), 
         i=1..nops(zetaValues)), 
         gridlines=true, 
 title="step response for different damping ratios", 
 labels=[typeset(t),typeset(y(t))]);

pict

Instead of using Simulate as above, another option is to use ResponsePlot which gives same plot as above.

restart; 
alias(DS=DynamicSystems): 
c:=ColorTools:-GetPalette("Niagara"): 
H:=(w,zeta)->w^2/(s^2+2*zeta*w*s+w^2): 
sys:= (w,zeta)-> 
    DS:-TransferFunction(H(w,zeta)): 
zetaValues := [seq(i,i=0.1..1.2,.2)]: 
 
sol:= map(i-> DS:-ResponsePlot(sys(1,zetaValues[i]), 
 Heaviside(t), 
 duration=14, 
 color=c[i], 
 legend=typeset(zeta,"=",zetaValues[i]), 
 legendstyle=[location=right] 
 ), 
  [seq(i,i=1..nops(zetaValues))] 
   ): 
 
plots:-display(sol,gridlines=true, 
 title="step response for different damping ratios", 
 labels=[typeset(t),typeset(y(t))] 
  );



1.12 View steady state error of 2nd order LTI system with changing undamped natural frequency

Problem: Given the transfer function

         ------ω2n-------
H  (s) = s2 + 2 ξω s + ω2
                 n     n

Display the output and input on the same plot showing how the steady state error changes as the un damped natural frequency ωn  changes. Do this for ramp and step input.

The steady state error is the difference between the input and output for large time. In other words, it the difference between the input and output at the time when the response settles down and stops changing.

Displaying the curve of the output and input on the same plot allows one to visually see steady state error.

Use maximum time of 10  seconds and ξ = 0.707  and change ωn  from 0.2  to 1.2  .

Do this for ramp input and for unit step input. It can be seen that with ramp input, the steady state do not become zero even at steady state. While with step input, the steady state error can become zero.

1.12.1 Mathematica

ramp input
Remove["Global`*"]; 
 
sys = TransferFunctionModel[w^2 /(s^2+2 z w  s+w^2 ),s]; 
z   = .707; 
nTrials = 6; 
maxTime = 10; 
 
tb = Table[Plot[ 
  Evaluate@{t,OutputResponse[sys/.w->0.2*i,t,t]}, 
     {t,0,maxTime}, 
     Frame->True, 
     PlotRange->All, 
     FrameLabel->{ 
      {"y(t)",None}, 
      {"t","Subscript[\[Omega], n]="<>ToString[0.2*i]} 
      } 
    ], 
    {i,1,nTrials}]; 
Grid[Partition[tb,2]]

pict

step input
tb = Table[Plot[ 
  Evaluate@{UnitStep[t], 
     OutputResponse[sys/.w->0.2*i,UnitStep[t],t]}, 
     {t,0,maxTime}, 
     Frame->True, 
     PlotRange->All, 
     FrameLabel->{ 
      {"y(t)",None}, 
      {"t","Subscript[\[Omega], n]="<>ToString[0.2*i]} 
      } 
    ], 
    {i,1,nTrials}]; 
 
Grid[Partition[tb,2]]

pict

1.13 Show the use of the inverse Z transform

These examples show how to use the inverse a Z transform.

1.13.1 example 1

Problem: Given

F (z) = --z---
        z − 1

find x [n ] = F −1 (z)  which is the inverse Ztransform.



Mathematica

Clear["Global`*"]; 
 
x[n_]:= InverseZTransform[z/(z-1),z,n] 
 
ListPlot[Table[{n,x[n]},{n,0,10}], 
  Frame->True, 
  FrameLabel->{{"x[n]",None}, 
      {"n","Inverse Ztransform of z/(z-1)"}}, 
  FormatType->StandardForm, 
  RotateLabel->False, 
  Filling->Axis, 
  FillingStyle->Red, 
  PlotRange->{Automatic,{0,1.3}}, 
  PlotMarkers->{Automatic,12}]

pict





Matlab

function e19() 
 
nSamples = 10; 
z = tf('z'); 
h = z/(z-1); 
[num,den] = tfdata(h,'v'); 
[delta,~] = impseq(0,0,nSamples); 
xn = filter(num,den,delta); 
 
stem(xn); 
title('Inverse Z transform of z/(z-1)'); 
xlabel('n'); ylabel('x[n]'); 
ylim([0 1.3]); 
set(gcf,'Position',[10,10,400,400]); 
 
end 
 
%function from Signal Processing by Proakis 
%corrected a little by me for newer matlab version 
function [x,n] = impseq(n0,n1,n2) 
% Generates x(n) = delta(n-n0); n1 <= n,n0 <= n2 
% ---------------------------------------------- 
% [x,n] = impseq(n0,n1,n2) 
% 
if ((n0 < n1) || (n0 > n2) || (n1 > n2)) 
    error('arguments must satisfy n1 <= n0 <= n2') 
end 
n = n1:n2; 
x = (n-n0) == 0; 
end

pict



1.13.2 example 2

Problem: Given

           5z
F (z) =  ------2-
         (z − 1 )

find          −1
x [n ] = F   (z)

In Mathematica analytical expression of the inverse Z transform can be generated as well as shown below



Mathematica

Clear["Global`*"]; 
 
x[n_]:= InverseZTransform[(5 z)/(z-1)^2,z,n]; 
Print["Inverse Z is ", x[n]]

Inverse Z is 5 n

ListPlot[Table[{n,x[n]},{n,0,10}], 
   Frame->True, 
   FrameLabel->{{"x[n]",None}, 
  {"n","Inverse Ztransform of (5 z)/(z-1)^2"}}, 
   FormatType->StandardForm, 
   RotateLabel->False, 
   Filling->Axis, 
   FillingStyle->Red, 
   PlotRange->{Automatic,{0,60}}, 
   PlotMarkers->{Automatic,12}, 
   ImageSize->350]

pict





Matlab

function e19_2() 
 
nSamples = 10; 
z = tf('z'); 
h = (5*z)/(z-1)^2; 
[num,den] = tfdata(h,'v'); 
[delta,~] = impseq(0,0,nSamples); 
xn = filter(num,den,delta); 
 
stem(xn); 
title('Inverse Z transform of 5z/(z-1)^2'); 
xlabel('n'); ylabel('x[n]'); 
ylim([0 60]); 
set(gcf,'Position',[10,10,400,400]); 
 
end

pict



1.14 Find the Z transform of sequence x(n)

1.14.1 example 1

Find the Z transform for the unit step discrete function

Given the unit step function x[n] = u[n]  defined as x =  {1,1,1,⋅⋅⋅}  for n ≥  0  , find its Z transform.



Mathematica

Remove["Global`*"]; 
ZTransform[UnitStep[n],n,z]
  Out[] = z/(-1+z)

Matlab

syms n 
pretty(ztrans(heaviside(n)))
      1
    ----- + 1/2
    z - 1



1.14.2 example 2

Find the Z transform for        ( )
x[n] =  13 n u (n ) + (0.9)n−3 u(n)



Mathematica

f[n_]:=((1/3)^n+(9/10)^(n-3))UnitStep[n]; 
ZTransform[f[n],n,z]
  Out[18]= z (3/(-1+3 z)+10000/(729 (-9+10 z)))

Matlab

syms n 
pretty(ztrans( 
 ((1/3)^n+(0.9)^(n-3)) 
       *heaviside(n)))
         100           1
    ------------- + ------- + 1729/1458
    81 (z - 9/10)   3 z - 1
  



1.15 Sample a continuous time system

Given the following system, sample the input and find and plot the plant output

pict

Use sampling frequency f = 1  Hz and show the result for up to 14  seconds. Use as input the signal u (t) = exp (− 0.3t)sin(2π(f∕3 )t)  .

Plot the input and output on separate plots, and also plot them on the same plot.

Mathematica

Clear["Global`*"]; 
 
nSamples =14; 
f =1.0; 
sampleTime =1/f; 
u[t_] := Sin[2*Pi*(f/3)*t]*Exp[-0.3*t] 
ud[n_]:= u[n*sampleTime]; 
plant = TransferFunctionModel[ 
            (3 s)/((s+1)(s+2)),s]; 
sysd  = ToDiscreteTimeModel[plant, 
         sampleTime, 
         z, 
         Method->"ZeroOrderHold"]

pict

opts = {Joined->True,PlotRange->All, 
  AspectRatio->1,InterpolationOrder->0, 
  ImageSize->200,Frame->True, 
  PlotRange->{{0,nSamples},{-0.5,0.5}}, 
  GridLines->Automatic,GridLinesStyle->Dashed}; 
 
inputPlot = ListPlot[Table[ud[k],{k,0,nSamples}], 
   Evaluate@opts, 
   FrameLabel->{{"u(t)",None}, 
             {"n","plant input u(nT)"}}, 
   PlotStyle->{Thick,Red}]; 
 
plantPlot=ListPlot[OutputResponse[sysd,Table[ud[k], 
 {k,0,nSamples}]], 
 Evaluate@opts, 
 FrameLabel->{{"y(t)",None}, 
             {"n","plant output y(nT)"}}, 
 PlotStyle->{Thick,Blue}]; 
 
Grid[{{inputPlot,plantPlot}}]

pict

Show[{inputPlot,plantPlot}, 
     PlotLabel->"input/output plot"]

pict

Matlab

clear all; close all; 
 
nSamples = 14; 
f = 1; 
T = 1/f; 
u=@(t) sin(2*pi*(f/3).*t).*exp(-0.3.*t); 
ud=@(n) u(n*T); 
U=ud(1:14);  %sampled input 
s=tf('s'); 
plant  = (3*s)/((s+1)*(s+2)); 
plantD = c2d(plant,T,'zoh')
  plantD =
  
       0.6976 z - 0.6976
    ------------------------
    z^2 - 0.5032 z + 0.04979
  
  Sample time: 1 seconds
  Discrete-time transfer function.
lsim(plantD,U,0:nSamples-1) 
grid

pict

1.16 Find closed loop transfer function from the open loop transfer function for a unity feedback

Problem: Given

       -------s------
L(s) = (s + 4)(s + 5)

as the open loop transfer function, how to find G (s)  , the closed loop transfer function for a unity feedback?

pict



Mathematica

Clear["Global`*"]; 
sys = TransferFunctionModel[ 
             s/((s+4)(s+5)),s]

pict

closedLoop = SystemsModelFeedbackConnect[sys];

pict

The system wrapper can be removed in order to obtain the rational polynomial expression as follows

First[Flatten[closedLoop[[1,1]]/ 
         closedLoop[[1,2]]]]
------s-------
s2 + 10s + 20

Matlab

clear all; close all; 
s=tf('s'); 
sys=s/((s+4)*(s+5)); 
closedloop=feedback(sys,1)
  Transfer function:
         s
  ---------------
  s^2 + 10 s + 20
  



1.17 Compute the Jordan canonical/normal form of a matrix A



Mathematica

Remove["Global`*"]; 
m={{3,-1,1,1,0,0}, 
   {1, 1,-1,-1,0,0}, 
  {0,0,2,0,1,1}, 
  {0,0,0,2,-1,-1}, 
  {0,0,0,0,1,1}, 
  {0,0,0,0,1,1}}; 
MatrixForm[m]
( 3− 1 1  1 0  0 )
|                |
| 1 1 − 1− 1 0 0 |
|| 0 0  2  0 1  1 ||
|| 0 0  0  2 − 1− 1||
( 0 0  0  0 1  1 )
  0 0  0  0 1  1

{a,b}=JordanDecomposition[m]; 
b
(        )
  000000
|| 021000 ||
|| 002000 ||
| 000210 |
|( 000021 |)

  000002

Matlab

clear all; 
A=[3 -1  1  1  0  0; 
   1  1 -1 -1  0  0; 
   0  0  2  0  1  1; 
   0  0  0  2 -1 -1; 
   0  0  0  0  1  1; 
   0  0  0  0  1  1;] 
 
jordan(A)
  ans =
   0     0     0     0     0     0
   0     2     1     0     0     0
   0     0     2     1     0     0
   0     0     0     2     0     0
   0     0     0     0     2     1
   0     0     0     0     0     2



1.18 Solve the continuous-time algebraic Riccati equation

Problem: Solve for X  in the Riccati equation

A ′X  + XA  −  XBR  −1B ′X +  C′C =  0

given

pict


Mathematica

Clear ["Global`*"]; 
a={{-3,2},{1,1}}; 
b={{0},{1}}; 
c={{1,-1}}; 
r={{3}}; 
sol=RiccatiSolve[{a,b},{Transpose[c].c,r}]; 
MatrixForm[N[sol]]
(                )
 0.5895171.82157
  1.82157 8.81884

Matlab

%needs control system 
clear all; close all; 
a = [-3 2;1 1]; 
b = [0 ; 1]; 
c = [1 -1]; 
r = 3; 
x = care(a,b,c'*c,r)
  x =
      0.5895    1.8216
      1.8216    8.8188



1.19 Solve the discrete-time algebraic Riccati equation

Problem: Given a continuous-time system represented by a transfer function

----1-----
s(s + 0.5)

convert this representation to state space and sample the system at sampling period of 1  second, and then solve the discrete-time Riccati equation.

The Riccati equation is given by

A ′X  + XA  −  XBR  −1B ′X +  C′C =  0

Let R  = [3]  .



Mathematica

Clear["Global`*"]; 
sys=TransferFunctionModel[1/(s(s+.5)),s]; 
dsys=ToDiscreteTimeModel[sys, 
        1,z,Method->"ZeroOrderHold"]

pict

ss = StateSpaceModel[dsys]

pict

a = ss[[1,1]]; 
b = ss[[1,2]]; 
c = ss[[1,3]]; 
d = ss[[1,4]]; 
r = {{3}}; 
DiscreteRiccatiSolve[{a,b}, 
               {Transpose[c].c,r}]; 
MatrixForm[%]
(                     )
   0.671414 − 0.977632
  − 0.977632  2.88699

Matlab

clear all; close all; 
s    = tf('s'); 
sys  = 1/(s*(s+0.5)); 
dsys = c2d(sys,1)
  dsys =
      0.4261 z + 0.3608
    ----------------------
    z^2 - 1.607 z + 0.6065
  
  Sample time: 1 seconds
  Discrete-time transfer function.
[A,B,C,D]=dssdata(dsys)
  A =
      1.6065   -0.6065
      1.0000         0
  B =
       1
       0
  C =
      0.4261    0.3608
  D =
       0
  
      2.8870   -0.9776
     -0.9776    0.6714
dare(A,B,C'*C,3)
  ans =
  
      2.8870   -0.9776
     -0.9776    0.6714



1.20 Display impulse response of H(z) and the impulse response of its continuous time approximation H(s)

Plot the impulse response of             2
H (z) = z∕(z −  1.4z + 0.5)  and using sampling period of T =  0.5  find continuous time approximation using zero order hold and show the impulse response of the system and compare both responses.

Mathematica

maxSimulationTime = 10; 
samplePeriod      = 0.5; 
tf = TransferFunctionModel[z/(z^2-1.4 z+0.5), 
          z, 
          SamplingPeriod->samplePeriod];

pict

Find its impulse response

discreteResponse=First@OutputResponse[tf,DiscreteDelta[k], 
                                           {k,0,maxSimulationTime}]

  {0.,1.,1.4,1.46,1.344,1.1516,0.94024,0.740536,
  0.56663,0.423015,0.308905,0.22096,0.154891,
  0.106368,0.0714694,0.0468732,0.0298878,
  0.0184063,0.0108249,0.00595175,0.00291999}

approximate to continuous time, use ZeroOrderHold

ctf = ToContinuousTimeModel[tf, s, Method -> "ZeroOrderHold"]

pict

Find the impulse response of the continuous time system

continouseTimeResponse=Chop@First@OutputResponse[ctf,DiracDelta[t],t]
          −0.693147t
− 1.25559e        (− 13.3012 θ(t)sin (0.283794t ) − 1.θ(t)cos(0.283794t))

continouseTimeResponse=Chop@First@OutputResponse[ctf,DiracDelta[t], 
                                         {t,0,maxSimulationTime}]

pict

Plot the impulse response of the discrete system

ListPlot[ 
discreteResponse, 
DataRange->{0,maxSimulationTime}, 
Filling->Axis, 
FillingStyle->{Red, 
  PlotMarkers->Graphics[{PointSize[0.03], 
  Point[{0,0}]}]},PlotRange->All, 
  ImageSize->300,ImagePadding->{{45,10},{35,10}}, 
  AspectRatio->0.9,AxesOrigin->{0,0},Frame->True,Axes->None, 
  ImageSize->300, 
  Frame->True 
  ]

pict

Plot the impulse response of the continuous system

Plot[samplePeriod *continouseTimeResponse,{t,0,maxSimulationTime}, 
     PlotRange->All,Frame->True,ImageSize->300]

pict

Plot both responses on the same plot

p = Show[ 
  ListPlot[ 
   discreteResponse, 
   Filling -> Axis, 
   FillingStyle -> {Red, 
     PlotMarkers -> Graphics[{PointSize[0.03], Point[{0, 0}]}]}, 
   PlotRange -> All, 
   DataRange -> {0, maxSimulationTime}, 
   ImageSize -> 300, 
   ImagePadding -> {{45, 10}, {35, 50}}, 
   AspectRatio -> 0.9, 
   AxesOrigin -> {0, 0}, 
   Frame -> True, 
   Axes -> None], 
  Plot[samplePeriod *continouseTimeResponse, {t, 0, maxSimulationTime}, 
   PlotRange -> All], 
  PlotRange -> All, 
  FrameLabel -> {{Style["response", 10], None}, 
     {Style["time (sec)", 10], 
    "impulse response of discrete system with its continouse time approximatio"}} 
  ]

pict

Do the same plot above, using stair case approximation for the discrete plot

Show[ 
ListPlot[ 
 discreteResponse, 
 Filling->None,PlotRange->All,DataRange->{0,maxSimulationTime}, 
 ImageSize->300,ImagePadding->{{45,10},{35,50}},AspectRatio->0.9, 
 AxesOrigin->{0,0},Frame->True,Axes->None, 
 Joined->True,InterpolationOrder->0,PlotStyle->Red], 
Plot[samplePeriod *continouseTimeResponse,{t,0,maxSimulationTime}, 
 PlotRange->All], 
 PlotRange->All,FrameLabel->{{Style["response",10],None}, 
 {Style["time (sec)",10],"impulse response"}}]

pict

Matlab

clear all; 
maxSimulationTime = 15; 
samplePeriod = 0.5; 
 
z=tf('z',samplePeriod); 
tf=z/(z^2-1.4*z+0.5) 
impulse(tf,0:samplePeriod:maxSimulationTime) 
 
ctf=d2c(tf,'zoh') 
 
hold on 
impulse(ctf,0:maxSimulationTime) 
title(['impulse response of discrete system with',... 
       ' its continuous time approximation']); 
xlabel('time (sec)'); 
ylabel('response');

pict

1.21 Find the system type given an open loop transfer function

Problem: Find the system type for the following transfer functions

  1. ss+21−s
  2. -s3+12
s −s
  3. s+1
 s5

To find the system type, the transfer function is put in the form  k∑ (s−si)
sM-∑ji(s−-sj)   , then the system type is the exponent M  . Hence it can be seen that the first system above has type one since the denominator can be written as  1
s (s − 1)  and the second system has type 2 since the denominator can be written as  2
s (s − 1)  and the third system has type 5. The following computation determines the type



Mathematica

  Clear["Global`*"];
  p=TransferFunctionPoles[TransferFunctionModel[
      (s+1)/(s^2-s),s]];
  Count[Flatten[p],0]
  
  p=TransferFunctionPoles[TransferFunctionModel[
      (s+1)/( s^3-s^2),s]];
  Count[Flatten[p],0]
  
  p=TransferFunctionPoles[
     TransferFunctionModel[(s+1)/ s^5 ,s]];
  Count[Flatten[p],0]
  
  Out[171]= 1
  Out[173]= 2
  Out[175]= 5

Matlab

  clear all;
  s=tf('s');
  [~,p,~]=zpkdata((s+1)/(s^2-s));
  length(find(p{:}==0))
  
  [~,p,~]=zpkdata((s+1)/(s^3-s^2));
  length(find(p{:}==0))
  
  
  [~,p,~]=zpkdata((s+1)/s^5);
  length(find(p{:}==0))
  
  >>
  ans =
       1
  
  ans =
       2
  
  ans =
       5



1.22 Find the eigenvalues and eigenvectors of a matrix

Problem, given the matrix

(          )
   1  2  3
(  4  5  6 )

   7  8  9

Find its eigenvalues and eigenvectors.



Mathematica

  Remove["Global`*"]
  (a = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}})
     // MatrixForm

(    )
  123
( 456)
  789

  N[Eigenvalues[a]]//MatrixForm

{ 3 (    √ --)  3 (    √ --)   }
  -- 5 +   33  ,-- 5 −   33  ,0
  2             2

  N[Eigenvectors[a]]//MatrixForm

(       √ --  (  √--)  )
  − −-15−√-33 4-6+√33- 1
||    33+7√-33  (33+7√ 33)  ||
|( − -15−-3√3-4-−6+-√33-1|)
    −33+7 33 −33+7 33
      1        − 2    1



Matlab

Matlab generated eigenvectors are such that the sum of the squares of the eigenvector elements add to one.

  clear all; close all;
  
  A=[1 2 3;4 5 6;7 8 9];
  
  [v,e]=eig(A)

  v =
     -0.2320   -0.7858    0.4082
     -0.5253   -0.0868   -0.8165
     -0.8187    0.6123    0.4082
  
  e =
     16.1168         0         0
           0   -1.1168         0
           0         0   -0.0000



1.23 Find the characteristic polynomial of a matrix

(          )
   1  2  3
(  4  5  6 )
   7  8  0



Mathematica

  a = {{1, 2, 3}, {4, 5, 6}, {7, 8, 0}};
  CharacteristicPolynomial[a,x]
  
  Out[6]= 27+72 x+6 x^2-x^3

Matlab

Note: Matlab generated characteristic polynomial coefficients are negative to what Mathematica generated.

  clear all;
  A=[1 2 3;4 5 6;7 8 0];
  
  p=poly(A)
  poly2str(p,'x')
  >>
  p =
   1.0000 -6.0000 -72.0000 -27.0000
  ans =
   x^3 - 6 x^2 - 72 x - 27



1.24 Verify the Cayley-Hamilton theorem that every matrix is zero of its characteristic polynomial

Problem, given the matrix

(      )
  1  2
  3  2

Verify that matrix is a zero of its characteristic polynomial. The Characteristic polynomial of the matrix is found, then evaluated for the matrix. The result should be the zero matrix.



Mathematica

  Remove["Global`*"]
  a = {{1,2},{3,2}};
  n = Length[a];
  p = CharacteristicPolynomial[a,x]

x2 − 3x − 4

  (-4 IdentityMatrix[n] - 3 a +
    MatrixPower[a,2])//MatrixForm
(00 )

 00

Another way to do the above is as follows

  a  = {{1,2},{3,2}};
  p  = CharacteristicPolynomial[a,x];
  cl = CoefficientList[p,x];
  Sum[MatrixPower[a,j-1] cl[[j]],
       {j,1,Length[cl]}]
(   )
 00
 00

Matlab

MATLAB has a build-in function polyvalm() to do this more easily than in Mathematica. Although the method shown in Mathematica can easily be made into a Matlab function

  clear all;
  
  A=[1 2;3 2];
  p=poly(A);
  poly2str(p,'x')
  polyvalm(p,A)
  
  ans =
     x^2 - 3 x - 4
  ans =
       0     0
       0     0



1.25 How to check for stability of a continuous system represented as a transfer function and state space

Problem: Given a system Laplace transfer function, check if it is stable, then convert to state space and check stability again. In transfer function representation, the check is that all poles of the transfer function (or the zeros of the denominator) have negative real part. In state space, the check is that the matrix A is negative definite. This is done by checking that all the eigenvalues of the matrix A have negative real part. The poles of the transfer function are the same as the eigenvalues of the A matrix. Use

           5s
sys =  -2----------
       s + 4s + 25

1.25.1 Checking stability using transfer function poles



Mathematica

  Remove["Global`*"];
  sys   = TransferFunctionModel[(5s)/(s^2+4s+25),s];
  poles = TransferFunctionPoles[sys]
  
  Out[47]= {{{-2-I Sqrt[21],-2+I Sqrt[21]}}}
  
  Re[#]&/@poles
  Out[42]= {{{-2,-2}}}
  
  Select[%,#>=0&]
  Out[44]= {}

Matlab

   clear all;
  s=tf('s');
  sys=5*s/( s^2+4*s+25);
  [z,p,k]=zpkdata(sys,'v');
  p
  >>
  p =
    -2.0000 + 4.5826i
    -2.0000 - 4.5826i
  find(real(p)>=0)
  >>
  ans =
     Empty matrix: 0-by-1



1.25.2 Checking stability using state space A matrix



Mathematica

  ss=StateSpaceModel[sys];
  a=ss[[1,1]]
  
  Out[49]= {{0,1},{-25,-4}}
  
  e=Eigenvalues[a]
  Out[50]= {-2+I Sqrt[21],-2-I Sqrt[21]}
  
  Re[#]&/@e
  Out[51]= {-2,-2}
  
  Select[%,#>=0&]
  Out[52]= {}

Matlab

  sys=ss(sys);
  [A,B,C,D]=ssdata(sys);
  A
  >>>
  A =
     -4.0000   -6.2500
      4.0000         0
  
  e=eig(A);
  e
  >>
  e =
    -2.0000 + 4.5826i
    -2.0000 - 4.5826i
  
  find(real(e)>=0)
  >>
  ans =
     Empty matrix: 0-by-1



1.26 Check continuous system stability in the Lyapunov sense

Problem: Check the stability (in Lyapunov sense) for the state coefficient matrix

     ⌊             ⌋
       0    1    0
A  = ⌈ 0    0    1 ⌉

       − 1 − 2  − 3

The Lyapunov equation is solved using lyap() function in MATLAB and LyapunovSolve[] function in Mathematica, and then the solution is checked to be positive definite (i.e. all its eigenvalues are positive).

We must transpose the matrix A  when calling these functions, since the Lyapunov equation is defined as   T
A  P + P A =  − Q  and this is not how the software above defines them. By simply transposing the A  matrix when calling them, then the result will be correct.



Mathematica

  Remove["Global`*"];
  (mat = {{0,1,0},{0,0,1},{-1,-2,-3}})
(          )
   0  1  0
(  0  0  1 )
  − 1− 2− 3

  p = LyapunovSolve[Transpose@mat,
      -IdentityMatrix[Length[mat]]];
  MatrixForm[N[p]]
(          )
  2.32.10.5
( 2.14.61.3 )
  0.51.30.6

  N[Eigenvalues[p]]

{6.18272, 1.1149,0.202375 }

Matlab

  clear all;
  
  A=[0   1   0
     0   0   1
     -1 -2  -3];
  
  p=lyap(A.',eye(length(A)))
  
  p =
    2.3   2.1  0.5
    2.1   4.6  1.3
    0.5   1.3  0.6
  
  e=eig(p)
  
  e =
       0.20238
         1.1149
         6.1827





Maple

  with(LinearAlgebra):
  A:=<<0,1,0;0,0,1;-1,-2,-3>>;
  p,s:=LyapunovSolve(A^%T,-<<1,0,0;0,1,0;0,0,1>>);
  Eigenvalues(p);

⌊                          ⌋
  6.18272045921436  +  0.0i
|| 1.11490451203192  +  0.0i ||
⌈                          ⌉
  0.202375028753723  + 0.0i



1.27 Given a closed loop block diagram, generate the closed loop Z transform and check its stability

Problem: Given the following block diagram, with sampling time T  = 0.1 sec  , generate the closed loop transfer function, and that poles of the closed loop transfer function are inside the unit circle

pict

System block diagram.



Mathematica

  Remove["Global`*"]
  plant  = TransferFunctionModel[s/(1.0+s),s];
  plantd = ToDiscreteTimeModel[
               plant,
               1,
               z,
               Method->"ZeroOrderHold"]

pict

  d = ToDiscreteTimeModel[
       TransferFunctionModel[(s+1)/(s+5.0),s],
            1,
            z,
            Method->"ZeroOrderHold"]

pict

  sys=SystemsModelSeriesConnect[d,plantd]

pict

  loopBack=SystemsModelFeedbackConnect[sys]

pict

Now generate unit step response

  ListPlot[OutputResponse[loopBack,
        Table[1, {k, 0, 20}]],
          Joined->True,
          PlotRange->All,
          InterpolationOrder->0,
          Frame->True,
          PlotStyle->{Thick,Red},
          GridLines->Automatic,
          GridLinesStyle->Dashed,
          FrameLabel->{{"y(k)",None},
            {"n","plant response"}},
          RotateLabel->False]

pict

  ss=StateSpaceModel[loopBack]

pict

  poles=TransferFunctionPoles[loopBack]

{0.543991 − 0.325556i, 0.543991  + 0.325556i}

  Abs[#]&/@poles

(                    )
 {0.633966,0.633966 }

Poles are inside the unit circle, hence stable.



Matlab

  clear all; close all
  
  s      = tf('s');
  plant  = s/(s+1);
  T      = 1; %sampeling time;
  plantd = c2d(plant,T,'zoh');
  d      = c2d( (s+1)/(s+5), T , 'zoh');
  % obtain the open loop
  sys    = series(d,plantd);
  % obtain the closed loop
  loop   = feedback(sys,1)

  loop =
  
     z^2 - 1.801 z + 0.8013
    ------------------------
    2 z^2 - 2.176 z + 0.8038
  
  Sample time: 1 seconds
  Discrete-time transfer function.

  step(loop)
  grid

pict

  [z,p,k]=zpkdata(loop,'v');
  fprintf('poles of closed loop discrete'...
    'system as inside unit circle\n');
  abs(p)

  ans =
      0.6340
      0.6340



1.28 Determine the state response of a system to only initial conditions in state space

Problem: Given a system with 2 states, represented in state space, how to determine the state change due some existing initial conditions, when there is no input forces?



Mathematica

  Remove["Global`*"];
  a = {{-0.5572, -0.7814}, {0.7814, 0}};
  c = {{1.9691, 6.4493}};
  sys=StateSpaceModel[{a,{{1},{0}},c,{{0}}}]

pict

  x0 = {1,0};
  {x1,x2} = StateResponse[{sys,x0},0,{t,0,20}];
  Plot[{x1,x2},{t,0,20},
    PlotRange->All,
    GridLines->Automatic,
    GridLinesStyle->Dashed,
    Frame->True,
    ImageSize->350,
    AspectRatio->1,
    FrameLabel->{{"y(t)",None},
      {"t","first and second state change with time"}},
      RotateLabel->False,
      PlotStyle->{Red,Blue},BaseStyle -> 12]

pict



Matlab

  clear all;
  
  A = [-0.5572   -0.7814;0.7814  0];
  B = [1;0];
  C = [1.9691  6.4493];
  x0 = [1 ; 0];
  
  sys = ss(A,B,C,[]);
  
  [y,t,x]=initial(sys,x0,20);
  plot(t,x(:,1),'r',t,x(:,2),'b');
  title('x1 and x2 state change with time');
  xlabel('t (sec)');
  ylabel('x(t)');
  grid
  set(gcf,'Position',[10,10,320,320]);

pict



1.29 Determine the response of a system to only initial conditions in state space

Problem: Given a system represented by state space, how to determine the response y(t)  due some existing initial conditions in the states. There is no input forces.



Mathematica

  Remove["Global`*"];
  a = {{-0.5572, -0.7814}, {0.7814, 0}};
  c = {{1.9691, 6.4493}};
  sys=StateSpaceModel[{a,{{1},{0}},c,{{0}}}]

pict

  x0 = {1,0};  (*initial state vector*)
  y  = OutputResponse[{sys,x0},0,{t,0,20}];
  Plot[y,{t,0,20},
    PlotRange->All,
    GridLines->Automatic,
    GridLinesStyle->Dashed,
    Frame->True,
    ImageSize->300,
    AspectRatio->1,
    FrameLabel->{{"y(t)",None},
       {"t","system response"}},
    RotateLabel->False,PlotStyle->Red]

pict

Matlab

  clear all;
  
  A = [-0.5572   -0.7814;0.7814  0];
  B = [1;0];
  C = [1.9691  6.4493];
  x0 = [1 ; 0];
  
  sys = ss(A,B,C,[]);
  
  [y,t,x]=initial(sys,x0,20);
  plot(t,y);
  title('y(t) plant response');
  xlabel('t (sec)');
  ylabel('y(t)');
  grid
  set(gcf,'Position',[10,10,320,320]);

pict



1.30 Determine the response of a system to step input with nonzero initial conditions

Problem: Given a system represented by state space, how to determine the response with nonzero initial conditions in the states and when the input is a step input?



Mathematica

  Remove["Global`*"];
  a = {{-0.5572, -0.7814}, {0.7814, 0}};
  c = {{1.9691, 6.4493}};
  sys=StateSpaceModel[{a,{{1},{0}},c,{{0}}}]

pict

  x0 = {1,0};
  y  = OutputResponse[{sys,x0},
           UnitStep[t],{t,0,20}];
  
  Plot[y,{t,0,20},PlotRange->{{0,20},{0,13}},
    GridLines->Automatic,
    GridLinesStyle->Dashed,
    Frame->True,
    ImageSize->300,
    AspectRatio->1,
    FrameLabel->{{"y(t)",None},
    {"t",
    "system response to initial conditions and step input"}},
    RotateLabel->False,
    PlotStyle->Red]

pict

Matlab

  A = [-0.5572   -0.7814;0.7814  0];
  B = [1;0];
  C = [1.9691  6.4493];
  x0 = [1 ; 0];
  
  sys = ss(A,B,C,[]);
  
  [y1,t] = initial(sys,x0,20);
  y2    = step(sys,t);
  
  plot(t,y1+y2,'r');
  
  title(['initial conditions and step'...
            'input response']);
  xlabel('t *sec)'); ylabel('y(t)');
  grid
  set(gcf,'Position',[10,10,320,320]);

pict



1.31 Draw the root locus from the open loop transfer function

Problem: Given L (s)  , the open loop transfer function, draw the root locus. Let

                  2
L (s) = ---------s-+-2s-+-4----------
        s(s + 4)(s + 6)(s2 + 1.4s + 1 )

Root locus is the locus of the closed loop dominant pole as the gain k  is varied from zero to infinity.



Mathematica

  Remove["Global`*"]
  sys = TransferFunctionModel[
   k*(s^2+2 s+4)/(s(s+4)(s+6)(s^2+1.4s+1)),s];
  RootLocusPlot[sys,{k,0,100},
      ImageSize->300,
      GridLines->Automatic,
      GridLinesStyle->Dashed,
      Frame->True,
      AspectRatio -> 1]

pict

Matlab

  clear all; close all;
  s=tf('s');
  sys=(s^2+2*s+4)/...
   (s*(s+4)*(s+6)*(s^2+1.4*s+1));
  
  rlocus(sys,0:100)
  set(gcf,'Position',[10,10,420,420]);

pict



1.32 Find eAt  where A  is a matrix

Mathematica

  ClearAll["Global`*"]
  mat={{0,1},
       {-2,-3}}
  MatrixExp[mat t];
  MatrixForm[%]
(  − e− 2t + 2e−t − e−2t + e− t)
     −2t     −t    − 2t    −t
   2e   −  2e    2e    − e

Now verify the result by solving for eAt  using the method would one would do by hand, if a computer was not around. There are a number of methods to do this by hand. The eigenvalue method, based on the Cayley Hamilton theorem will be used here. Find the eigenvalues of |A − λI |

  m = mat-lambda IdentityMatrix[Length[mat]]
(              )
  − λ     1
   − 2 − λ − 3

Det[m]

 2
λ  + 3λ + 2

  sol=Solve[%==0,lambda]
      Out[15]= {{lambda->-2},{lambda->-1}}

  eig1=lambda/.sol[[1]]
  eig2=lambda/.sol[[2]]
  
        Out[16]= -2
        Out[17]= -1

  (*setup the equations to find b0,b1*)
  eq1 = Exp[eig1 t]==b0+b1 eig1;
  eq2 = Exp[eig2 t]==b0+b1 eig2;
  sol = First@Solve[{eq1,eq2},{b0,b1}]
{          (       )            (     )}
 b0 →  e− 2t 2et − 1  ,b1 →  e−2t et − 1

  (*Now find e^At*)
  b0=b0/.sol[[1]]
 −2t(  t    )
e    2e −  1

  b1=b1/.sol[[2]]
     (     )
e −2t et − 1

  b0 IdentityMatrix[Length[mat]]+b1 mat;
  Simplify[%]
  MatrixForm[%]
(                                    )
    e−2t(− 1 + 2et)   e−2t(− 1 + et)
   − 2e−2t(− 1 + et) − e−2t(− 2 + et)

The answer is the same given by Mathematica’s command MatrixExp[]

Matlab

  syms t
  A=[0 1;-2 -3];
  expm(t*A)
  
  ans =
  
  [2/exp(t)-1/exp(2*t),1/exp(t)-1/exp(2*t)]
  [2/exp(2*t)-2/exp(t),2/exp(2*t)-1/exp(t)]
  
  pretty(ans)
  
  +-                                        -+
  | 2 exp(-t)-  exp(-2 t),exp(-t)-exp(-2 t)  |
  |                                          |
  | 2 exp(-2 t)-2 exp(-t),2 exp(-2 t)-exp(-t)|
  +-                                        -+

1.33 Draw root locus for a discrete system

Problem: Given the open loop for a continuous time system as

         5s + 1
sys = -2---------
      s  + 2s + 3

convert to discrete time using a sampling rate and display the root locus for the discrete system.



Mathematica

  Remove["Global`*"];
  sys=TransferFunctionModel[k(5s+1)/(s^2+2s+3),s];
  samplePeriod=0.5; (*sec*)
  sysd=ToDiscreteTimeModel[ sys,samplePeriod,z]
  
    Out[] = TransferFunctionModel[
           {{{0.5*k*(1. + z)*(-9.5 + 10.5*z)}},
           2.75 - 6.5*z + 6.75*z^2},
           z, SamplingPeriod -> 0.5]

pict

  p1 = RootLocusPlot[sysd,{k,0,100},
         ImageSize->300,
         GridLines->Automatic,
         GridLinesStyle->Dashed,
         Frame->True,
         AspectRatio->1,
         PlotPoints->200,
         PlotStyle->PointSize[0.01]];
  
  p2=ParametricPlot[{Sin[t],Cos[t]},{t,0,2Pi},
         AspectRatio->1,
         GridLines->Automatic,
         GridLinesStyle->Dashed];
  
  Show[p2,p1]

pict





Matlab

  clear all; close all;
  s    = tf('s');
  sys  = (5*s+1)/(s^2+2*s+3);
  Ts   = 0.5; %sample time
  sysd = c2d(sys,Ts,'zoh');
  
  rlocus(sysd)
  grid
  set(gcf,'Position',[10,10,520,520]);

pict



1.34 Plot the response of the inverted pendulum problem using state space

Problem: Given the inverted pendulum shown below, use state space using one input (the force on the cart) and 2 outputs (the cart horizontal displacement, and the pendulum angle. Plot each output separately for the same input.

pict



Mathematica

  Remove["Global`*"];
  a = {{0, 1, 0, 0},
       {0, 0, ((-m)*g)/M, 0},
       {0, 0, 0, 1},
       {0, 0, ((M + m)*g)/(M*L), 0}};

(           )
  01   0gm  0
|| 00  − M--0||
( 00   0   1)
  00 g(m+M-)0
      LM

  b = {{0}, {1/M}, {0}, {-(1/(M*L))}};
  MatrixForm[b]

(      )
    0
||  M1  ||
(   0  )
  − -1-
    LM

  c = {{1, 0, 0, 0}, {0, 0, 1, 0}};
  MatrixForm[b]

(      )
  1000
  0010

  sys=StateSpaceModel[{a,b ,c}]

pict



It is now possible to obtain the response of the system as analytical expression or an interpolatingFunction.

It is much more efficient to obtain the response as interpolatingFunction. This requires that the time domain be given.

Here is example of obtaining the analytical expression of the response



  sys=sys/.{g->9.8,m->1,M->20,L->1};
  y=Chop[OutputResponse[sys,UnitStep[t],t]]

e−6.41561t(0.0238095e6.41561tt2θ(t)+
0.000115693e3.2078tθ(t)− 0.000231385e6.41561tθ(t)
+0.000115693e9.62341tθ(t))
,e−6.41561t(− 0.00242954e3.2078tθ(t)+ 0.00485909
e6.41561tθ(t)−
0.00242954e9.62341tθ(t))

  Plot[y[[1]],{t,0,3},
    PlotLabel->"y(t) response",
    FrameLabel->{"time (sec)","y (meter)"},
    Frame->True,
    PlotRange->All,
    ImageSize->300,
    GridLines->Automatic,
    GridLinesStyle->Dashed,
    RotateLabel->False]

pict

  Plot[ y[[2]],{t,0,3},
   PlotLabel->"\[Theta](t) response",
   FrameLabel->{"time (sec)","\[Theta] (radians)"},
   Frame->True,
   PlotRange->All,
   ImageSize->300,
   GridLines->Automatic,
   GridLinesStyle->Dashed,
   RotateLabel->False]

pict





Matlab

  clear all; close all;
  
  m=1; M=20; L=1; g=9.8;
  
  A=[0 1      0           0;
     0 0   -m*g/M         0;
     0 0      0           1;
     0 0  (M+m)*g/(M*L)   0
     ];
  B=[0  1/M  0 -1/(M*L)]';
  C=[1 0 0 0;
     0 0 1 0];
  D=[0];
  
  sys=ss(A,B,C(1,:),0);
  step(sys,3);
  grid on;
  set(gcf,'Position',[10,10,320,320]);
  title('y(t) due to step input');

pict

  figure;
  sys=ss(A,B,C(2,:),0);
  step(sys,3);
  title('\theta(t) due to step input');
  grid on;
  set(gcf,'Position',[10,10,320,320]);

pict



1.35 How to build and connect a closed loop control systems and show the response?

1.35.1 example 1, single input, single output closed loop

Given the following simple closed loop system, show the step response. Let mass M  =  1kg  , damping coefficient c = 1newton -seconds per meter  and let the stiffness coefficient be k =  20newton  per meter  .

pict

Using propertional controller J (s) = kp  where kp = 400  . First connect the system and then show y(t)  for 5 seconds when the reference yr(t)  is a unit step function.



Mathematica

  m = 1; c = 1; k = 20; kp = 400;
  plant      = TransferFunctionModel[1/(m s^2 + c s + k), s];
  controller = TransferFunctionModel[kp, s];
  sys        = SystemsModelSeriesConnect[plant, controller];
  sys        = SystemsModelFeedbackConnect[sys];
  o          = OutputResponse[sys, UnitStep[t], t];
  Plot[o, {t, 0, 10}, Frame -> True, PlotRange -> All,
    FrameLabel -> {{"y(t)", None},
      {"t (sec)", "response of closed loop"}},
      BaseStyle -> 14]

pict





Matlab

  close all;
  m = 1; c = 1; k = 20;
  s          = tf('s');
  plant      = 1/(m*s^2+c*s+k);
  controller = 400;
  sys        = series(controller,plant);
  sys        = feedback(sys,1);
  [Y,T]      = step(sys,0:.01:10);
  plot(T,Y);
  xlabel('time (sec)');
  ylabel('y(t)');
  title('response of closed loop');

pict

Another way to do the above is

  m=1;
  c=1;
  k=20;
  s=tf('s');
  sys3=1/(m*s^2+c*s+k);
  sys2=400;
  sys1=1;
  sys=append(sys1,sys2,sys3);
  Q=[2 1 -3;3 2 0];
  input=[1];
  output=[2 3];
  sys=append(sys1,sys2,sys3);
  sys=connect(sys,Q,input,output);
  T=0:.01:10;
  [Y,T,X]=step(sys,T);



1.36 Compare the effect on the step response of a standard second order system as ζ  changes

Problem: Given a standard second order system defined by the transfer function

         ------ω2n-------
G (s) =  s2 + 2ζ ωns + ω2n

Change ζ  and plot the step response for each to see the effect of changing ζ  (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

Y--(s)    ------ω2n-------
U (s) =  s2 + 2 ζω s + ω2
                 n     n

Where Y (s)  and U (s)  are the Laplace transforms of the output and input respectively.

Hence the differential equation, assuming zero initial conditions becomes

 ′′            ′      2          2
y  (t) + 2ζωn  y (t) + ωn y (t) = ω n u (t)

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

[  ′]   [             ] [  ]   [   ]
  x1  =    0      1      x1  +   0  u (t)
  x′2     − ω2n  − 2ζωn    x2     ω2n

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

[  ]   [             ] [  ]   [   ]
 x′1       0      1      x1      0
 x′  =   − ω2  − 2ζω    x   +  ω2
  2         n       n    2       n

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]

pict

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];

pict



1.37 Plot the dynamic response factor Rd  of a system as a function of     -ω
r = ωn   for different damping ratios

Problem: Plot the standard curves showing how the dynamic response Rd  changes as r = ωωn   changes. Do this for different damping ratio ξ  . Also plot the phase angle.

These plots are result of analysis of the response of a second order damped system to a harmonic loading. ω  is the forcing frequency and ωn  is the natural frequency of the system.



Mathematica

  Rd[r_,z_]:=1/Sqrt[(1-r^2)^2+(2 z r)^2];
  
  phase[r_,z_]:=Module[{t},
     t=ArcTan[(2z r)/(1-r^2)];
     If[t<0,t=t+Pi];
     180/Pi t
  ];
  
  plotOneZeta[z_,f_] := Module[{r,p1,p2},
   p1 = Plot[f[r,z],{r,0,3},PlotRange->All,
         PlotStyle->Blue];
  
   p2 = Graphics[Text[z,{1.1,1.1f[1.1,z]}]];
   Show[{p1,p2}]
  ];
  
  p1 = Graphics[{Red,Line[{{1,0},{1,6}}]}];
  p2 = Map[plotOneZeta[#,Rd]&,Range[.1,1.2,.2]];
  
  Show[p2,p1,
   FrameLabel->{{"Subscript[R, d]",None},
   {"r= \[Omega]/Subscript[\[Omega], n]",
    "Dynamics Response vs. Frequency ratio for different \[Xi]"}},
   Frame->True,
   GridLines->Automatic,GridLinesStyle->Dashed,
   ImageSize -> 300,AspectRatio -> 1]

pict

  p = Map[plotOneZeta[#,phase]&,Range[.1,1.2,.2]];
  Show[p,FrameLabel->{{"Phase in degrees",None},
    {"r= \[Omega]/Subscript[\[Omega], n]",
    "Phase vs. Frequency ratio for different \[Xi]"}},
    Frame->True,
    GridLines->Automatic,GridLinesStyle->Dashed,
    ImageSize->300,AspectRatio->1]

pict



1.38 How to find closed loop step response to a plant with a PID controller?

Find and plot the step response of the plant ---1---
s2+2s+1   connected to a PID controller with P  = 10,I = 3.7,D  = 0.7  . Use negative closed loopback.



Mathematica

  plant = TransferFunctionModel[1/(s^2 + 2*s + 1), s];
  kip   = 10; ki = 3.7; kid = 0.7;
  pid   = TransferFunctionModel[
            (kip*s + ki + kid*s^2)/s, s];
  openLoop   = SystemsModelSeriesConnect[
     TransferFunctionModel[plant], pid];
  closedLoop = SystemsModelFeedbackConnect[
                    openLoop];
  input  = UnitStep[t];
  output = OutputResponse[closedLoop, input, t];
  
  Plot[{input, output}, {t, 0, 5}, PlotRange ->
   All, GridLines -> Automatic,
   GridLinesStyle -> Directive[LightGray, Dashed],
   Frame -> True,
   FrameLabel -> {{"y(t)", None},
      {"t (sec)", "Step response"}},
   BaseStyle -> 12]

pict



Matlab

  close all;
  clear all;
  s        = tf('s');
  plant    = 1/(s^2 + 2*s + 1)
  kip      = 10; ki = 3.7; kid = 0.7;
  thePID   = pid(kip,kid,kid);
  BLKSYS   = append(thePID,plant);
  Q        = [2 1; 1 -2];
  closedLoop = connect(BLKSYS,Q,1,2);
  [y,t]   = step(closedLoop);
  plot(t(1:140),y(1:140));
  xlabel('t (sec)'); ylabel('y(t)');
  title('step response');
  grid

pict



1.39 How to make Nyquist plot?

Nyquist command takes as input the open loop transfer function (not the closed loop!) and generates a plot, which was can look at to determine if the closed loop is stable or not. The closed loop is assumed to be unity feedback. For example, if the open loop is G (s)  , then we know that the closed loop transfer function is  G(s)
1+G-(s)-   . But we call Nyquist with G (s)  .

Here are two examples.

1.39.1 Example 1

Given G(s) =  1−-s0.2s  generate Nyquist plot.



Mathematica

  NyquistPlot[s/(1 - 0.2 s),
     NyquistGridLines -> Automatic]

pict



Matlab

  clear all; close all;
  nyquist([1 0],[-0.2 1])
  grid

pict



1.39.2 Example 2

Given         ---5(1−0.5s)----
G(s) =  s(1+0.1s)(1−0.25s)   generate Nyquist plot.



Mathematica

  NyquistPlot[5(1-0.5 s)/(s(1 + 0.1 s)(1 - 0.25 s))]

pict



Matlab

  clear all; close all;
  s=tf('s');
  sys=5*(1-0.5*s)/(s*(1+0.1*s)*(1-0.25*s))
  
           2.5 s - 5
    ------------------------
    0.025 s^3 + 0.15 s^2 - s
  
  Continuous-time transfer function.
  
  >> nyquist(sys)

pict



However, there is a better function to do this called nyquist1.m which I downloaded and tried. Here is its results



  clear all; close all;
  nyquist1()

pict



Chapter 2
Linear algebra, Linear solvers, common operations on vectors and matrices

 2.1 Multiply matrix with a vector
 2.2 Insert a number at specific position in a vector or list
 2.3 Insert a row into a matrix
 2.4 Insert a column into a matrix
 2.5 Build matrix from other matrices and vectors
 2.6 Generate a random 2D matrix from uniform (0 to 1) and from normal distributions
 2.7 Generate an n by m zero matrix
 2.8 Rotate a matrix by 90 degrees
 2.9 Generate a diagonal matrix with given values on the diagonal
 2.10 Sum elements in a matrix along the diagonal
 2.11 Find the product of elements in a matrix along the diagonal
 2.12 Check if a Matrix is diagonal
 2.13 Find all positions of elements in a Matrix that are larger than some value
 2.14 Replicate a matrix
 2.15 Find the location of the maximum value in a matrix
 2.16 Swap 2 columns in a matrix
 2.17 Join 2 matrices side-by-side and on top of each others
 2.18 Copy the lower triangle to the upper triangle of a matrix to make symmetric matrix
 2.19 extract values from matrix given their index
 2.20 Convert N  by M  matrix to a row of length N M
 2.21 find rows in a matrix based on values in different columns
 2.22 Select entries in one column based on a condition in another column
 2.23 Locate rows in a matrix with column being a string
 2.24 Remove set of rows and columns from a matrix at once
 2.25 Convert list of separated numerical numbers to strings
 2.26 Obtain elements that are common to two vectors
 2.27 Sort each column (on its own) in a matrix
 2.28 Sort each row (on its own) in a matrix
 2.29 Sort a matrix row-wise using first column as key
 2.30 Sort a matrix row-wise using non-first column as key
 2.31 Replace the first nonzero element in each row in a matrix by some value
 2.32 Perform outer product and outer sum between two vector
 2.33 Find the rank and the bases of the Null space for a matrix A
 2.34 Find the singular value decomposition (SVD) of a matrix
 2.35 Solve Ax  = b
 2.36 Find all nonzero elements in a matrix
 2.37 evaluate f(x) on a vector of values
 2.38 generates equally spaced N points between x1   and x2
 2.39 evaluate and plot a f(x,y)  on 2D grid of coordinates
 2.40 Find determinant of matrix
 2.41 Generate sparse matrix with n  by n  matrix repeated on its diagonal
 2.42 Generate sparse matrix for the tridiagonal representation of the classic second difference operator on 1D grid
 2.43 Generate sparse matrix for the Laplacian differential operator   2
∇  u  for 2D grid
 2.44 Generate sparse matrix for the Laplacian differential operator ∇2u  for 3D grid
 2.45 Generate the adjugate matrix for square matrix
 2.46 Multiply each column by values taken from a row
 2.47 extract submatrix from a larger matrix by removing row/column
 2.48 delete one row from a matrix
 2.49 delete one column from a matrix
 2.50 generate random matrix so that each row adds to 1
 2.51 generate random matrix so that each column adds to 1
 2.52 sum all rows in a matrix
 2.53 sum all columns in a matrix
 2.54 find in which columns values that are not zero
 2.55 How to remove values from one vector that exist in another vector
 2.56 How to find mean of equal sized segments of a vector
 2.57 find first value in column larger than some value and cut matrix from there
 2.58 make copies of each value into matrix into a larger matrix
 2.59 repeat each column of matrix number of times
 2.60 How to apply a function to each value in a matrix?
 2.61 How to sum all numbers in a list (vector)?
 2.62 How to find maximum of each row of a matrix?
 2.63 How to find maximum of each column of a matrix?
 2.64 How to add the mean of each column of a matrix from each column?
 2.65 How to add the mean of each row of a matrix from each row?
 2.66 Find the different norms of a vector
 2.67 Check if a matrix is Hermite
 2.68 Obtain the LU decomposition of a matrix
 2.69 Linear convolution of 2 sequences
 2.70 Circular convolution of two sequences
 2.71 Linear convolution of 2 sequences with origin at arbitrary position
 2.72 Visualize a 2D matrix
 2.73 Find the cross correlation between two sequences
 2.74 Find orthonormal vectors that span the range of matrix A
 2.75 Solve A x= b and display the solution
 2.76 Determine if a set of linear equations A x= b, has a solution and what type of solution
 2.77 Given a set of linear equations automatically generate the matrix A  and vector b  and solve Ax = b
 2.78 Convert a matrix to row echelon form and to reduced row echelon form
 2.79 Convert 2D matrix to show the location and values

Mathematica and Matlab allows one to do pretty much the same operations in the area of linear algebra and matrix manipulation. Two things to keep in mind is that Mathematica uses a more general way to store data.

Mathematica uses ragged arrays or a list of lists. This means rows can have different sizes. (these are the lists inside the list).

Hence each list (or row) can be a different size. In Matlab a matrix must have same size rows and same size columns.

In Matlab one can also have ragged arrays, these are called cells. In Mathematica, there is one data structure for both.

Another thing to keep in mind is that Matlab, due to its Fortran background is column major when it comes to operations on matrices. This simple example illustrate this difference. Suppose we have a matrix A  of 3 rows, and want to find the location where A (i,j) = 2  where i  is the row number and j  is the column number. Given this matrix

     (  1  2  3  4 )
     (             )
A =     2  3  1  5
        5  6  7  2

Then the result of using the find() command in Matlab is

  A=[1 2 3 4;
     2 3 1 5;
     5 6 7 2];
  [I,J]=find(A==2)
  I =
       2
       1
       3
  J =
       1
       2
       4

We see that the Matlab result gives the order of the rows it found the element at based on searching column wise since it lists the second row first in its result. Compare this to Mathematica Position[] command

  mat = {{1, 2, 3, 4},
         {2, 3, 1, 5},
         {5, 6, 7, 2}}
  Position[mat, 2]

which gives

  {{1,2},
   {2,1},
   {3,4}}

We see from the above that Mathematica searched row-wise.

Mathematica use for matrix manipulate will take more time to master compared to Matlab, since Mathematica data structure is more general and little more complex (ragged arrays) compared to Matlab since it also has to support symbolics in its commands and not just numbers

In Maple, was can use the following short cuts to enter vectors and matrices: For row vector:  v:=<1|2|3|4> and for column vector  v:=<1,2,3,4> and for matrix of say 2 rows and 3 columns  A:=<1,2|3,4|5,6> so | acts as column separator. There are other ways to do this (as typical in Maple), but I find the above the least confusing. For transpose we can use  A^%T

2.1 Multiply matrix with a vector

How to perform the following matrix/vector multiplication?

  [ 1 2 3 ]   [ 1 ]
  [ 4 5 6 ] * [ 2 ]
  [ 7 8 8 ]   [ 3 ]

In Mathematica the dot "." is used for Matrix multiplication. In Matlab the "*" is used. In Fortran, MATMUL is used.



Mathematica

  a={{1,2,3},{4,5,6},{7,8,8}};
  x={1,2,3};
  
  a.x
  Out[7]= {14,32,47}

Matlab

  A=[1 2 3;4 5 6;7 8 9];
  x=[1; 2; 3];
  c=A*x
  
      14
      32
      50



Ada

  with Ada.Text_Io; use Ada.Text_Io;
  with Ada.Numerics.Real_Arrays;
       use Ada.Numerics.Real_Arrays;
  
  procedure t1 is
    A : constant real_matrix :=
            (( 1.0,  2.0,  3.0),
            ( 4.0,  5.0,  6.0),
            ( 7.0,  8.0,  9.0));
  
    V : constant real_vector := (1.0,2.0,3.0);
    procedure Put (X : real_vector) is
      begin
          FOR e of X LOOP
              put_line(float'image(e));
          END LOOP;
    end put;
  begin
      put(A*v);
  end t1;

compile and run

  >gnatmake -gnat2012 t1.adb
  gcc -c -gnat2012 t1.adb
  gnatbind -x t1.ali
  gnatlink t1.ali
  >./t1
   1.40000E+01
   3.20000E+01
   5.00000E+01

Fortran

  PROGRAM main
  IMPLICIT NONE
  integer, parameter :: dp = kind(1.D0)
  integer, parameter :: N=3
  real(kind=dp), parameter, DIMENSION(N, N) :: &
    A=DBLE(reshape([ 1.0, 2.0, 3.0,&
                     4.0, 5.0, 6.0,&
                     7.0, 8.0, 9.0], shape(A), &
                                       order=[2,1]))
  
  real(kind=dp), parameter, DIMENSION(N, 1) :: &
             x=DBLE(reshape([1.0, 2.0, 3.0], shape(x)))
  
  real(kind=dp), DIMENSION(N, 1) :: result
  integer, DIMENSION(2)    :: siz
  
  result = MATMUL(A,x)
  siz    = SHAPE(result)
  Write(*,*) result
  
  print *, "number of rows = ",siz(1)
  print *, "number of columns = ",siz(2)
  
  end program

compile and run

  >gfortran  -fcheck=all -Wall -Wconversion
     -Wextra -Wconversion-extra -pedantic test.f90
  >./a.out
  14.000000000000000  32.000000000000000 50.000000000000000
   number of rows =            3
   number of columns =            1

Maple

  A:=Matrix([[1,2,3],[4,5,6],[7,8,9]]):
  x:=Vector([1,2,3]): #default is column vector
  c:=A.x;
⌊    ⌋
  14
⌈ 32 ⌉

  50

Python

  import numpy as np
  mat=np.array([[1,2,3],[4,5,6],[7,8,8]])
  b=np.array([1,2,3])
  b.shape=(3,1)
  r=dot(mat,b)
  
  
  array([[14],
         [32],
         [47]])
  
  

2.2 Insert a number at specific position in a vector or list

The problem is to insert a number into a vector given the index.



Mathematica

  vec={1,2,3,4};
  vec=Insert[vec,99,3];
  vec
  
  Out[11]= {1,2,99,3,4}

Matlab

  A=[1 2 3 4];
  A=[ A(1:2) 99 A(3:end) ]
  
  >>
  A =
       1     2    99     3     4



Fortran

  program f
    implicit none
    REAL(4), ALLOCATABLE :: A(:)
    A=[1,2,3,4];
    A=[A(1:2), 99.0, A(3:)]
    Print *,A
  end program f
  
  >gfortran -std=f2008 t2.f90
  >./a.out
   1.0000000 2.0000000 99.000000
            3.0000000 4.0000000

Maple

  v:=Vector[row]([1,2,3,4]);
  v:=Vector[row]([v[1..2],99,v[3..]]);

But I find using <> notation easier, so will use that from now on.

  v:=<1|2|3|4>;
  v:=<v(1..2)|99|v(3..)>;
  v :=[ 1 2 99 3 4]



Python Note that Python uses zero index.

  import numpy as np
  b=np.array([1,2,3,4])
  b=numpy.insert(b,2,99)
  b
  
     Out[86]: array([ 1,  2, 99,  3,  4])



2.3 Insert a row into a matrix

The problem is to insert a row into the second row position in a 2D matrix



Mathematica

  mat={{1,2,3},
        {4,5,6},
        {7,8,9}}
  mat = Insert[mat,{90,91,92},2]
  >>
   {{1,2,3},
    {90,91,92},
    {4,5,6},
    {7,8,9}}

Matlab

  A=[1 2 3;
     4 5 6;
     7 8 9]
  A=[ A(1,:); [90 91 92];  A(2:end,:) ]
  >>
       1     2     3
      90    91    92
       4     5     6
       7     8     9



Fortran

  program t3
    implicit none
    integer i,nRow,nCol
    integer(4), ALLOCATABLE ::A(:,:)
  
    A = reshape ([1,2,3,   &
                  4,5,6,   &
                  7,8,9], [3,3], order=[2,1])
    nRow=size(A,1)
    nCol=size(A,2)
  
    print *, 'before'
    do i=LBOUND(A, 1),UBOUND(A, 1)
        print *,A(i,:)
    end do
  
    A = reshape ([A(1,:),    &
        90,91,92,  &
        A(2:,:)], [nRow+1,nCol] , order=[2,1])
    print *, 'after'
    do i=LBOUND(A, 1),UBOUND(A, 1)
       print *,A(i,:)
    end do
  
  end program t3

Compile and run

  >gfortran -std=f2008 t3.f90
  >./a.out
   before
             1           2           3
             4           5           6
             7           8           9
   after
             1           2           3
            90          91          92
             4           7           5
             8           6           9

Maple

Using <<>> notation

  A:=< <1|2|3>,
       <4|5|6>,
       <7|8|9>>;
  
  A:=< A(1..2,..),<90|91|92>, A(3,..)>;
  
  [ 1  2   3
    4  5   6
    90 91  92
    7  8   9]

Using Matrix/Vector notation

  A:=Matrix([ [1,2,3],[4,5,6],[7,8,9]]);
  A:=Matrix([ [A[1..2,..]],
              [Vector[row]([91,92,92])],
              [A[3,..]
            ] ]);



Python

  import numpy as np
  mat=np.array([[1,2,3],[4,5,6],[7,8,9]])
  mat=numpy.insert(mat,1,[90,91,92],0)
  
  Out[98]:
  array([[ 1,  2,  3],
         [90, 91, 92],
         [ 4,  5,  6],
         [ 7,  8,  9]])



2.4 Insert a column into a matrix

The problem is to insert a column into the second column position in a 2D matrix.



Mathematica

  mat = {{1, 2, 3},
         {4, 5, 6},
         {7, 8, 9}};
  
  mat = Insert[Transpose[mat],
           {90, 91, 92}, 2];
  mat = Transpose[mat]
  >>
   {{1,90,2,3},
    {4,91,5,6},
    {7,92,8,9}}
  

Matlab

  A=[1 2 3;
     4 5 6;
     7 8 9];
  
  A=[A(:,1)  [90 91 92]'   A(:,2:end) ]
  
  >>
       1    90     2     3
       4    91     5     6
       7    92     8     9



Fortran

  program t4
    implicit none
    integer(4) ::i
    integer(4), ALLOCATABLE ::A(:,:)
  
    A = reshape ([1,2,3,   &
                  4,5,6,   &
                  7,8,9], [3,3], order=[2,1])
  
    A = reshape ([A(:,1), [90,91,92] ,
                    A(:,2:)] , [3,4])
  
    do i=LBOUND(A, 1),UBOUND(A, 1)
       print *,A(i,:)
    end do
  
  end program t4
  
  >gfortran t4.f90
  >./a.out
        1          90           2           3
        4          91           5           6
        7          92           8           9

Maple

  A:=< <1|2|3>,<4|5|6>,<7|8|9>>;
  A:=< A(..,1)|<90,91,92>|A(..,2..)>;
  
  [     1    90     2     3
        4    91     5     6
        7    92     8     9     ]
  

Using Matrix/Vector notation

  A:=Matrix([ [1,2,3],[4,5,6],[7,8,9]]);
  A:=Matrix([ [A[..,1],
              Vector[column]([91,92,92]),
              A[..,2..] ]]);



Python

  import numpy as np
  mat=np.array([[1,2,3],[4,5,6],[7,8,9]])
  mat=numpy.insert(mat,1,[90,91,92],1)
  
  Out[101]:
  array([[ 1, 90,  2,  3],
         [ 4, 91,  5,  6],
         [ 7, 92,  8,  9]])
  



2.5 Build matrix from other matrices and vectors

2.5.1 First example

Given column vectors      ⌊    ⌋
        1
v1 = ⌈  2 ⌉
        3 and      ⌊    ⌋
        4
v2 = ⌈  5 ⌉
        6 and      ⌊    ⌋
        7
v3 = ⌈  8 ⌉
        9 and      ⌊     ⌋
        10
v4 = ⌈  11 ⌉
        12 generate the matrix                    ⌊ 1   4 ⌋
                   |       |
                   | 2   5 |
     [        ]    ||       ||
        v1 v2      || 3   6 ||
m  =    v3 v4   =  | 7  10 |
                   ||       ||
                   |⌈ 8  11 |⌉

                     9  12

Matlab was the easiest of all to do these operations with. No surprise, as Matlab was designed for Matrix and vector operations. But I was surprised that Maple actually had good support for these things, using its <> notation, which makes working with matrices and vectors much easier.



Mathematica

The command ArrayFlatten is essential for this in Mathematica.

  v1 = {1, 2, 3}; v2 = {4, 5, 6};
  v3 = {7, 8, 9}; v4 = {10, 11, 12};
  m = ArrayFlatten[
        {{Transpose@{v1}, Transpose@{v2}},
        {Transpose@{v3}, Transpose@{v4}}}
       ]

(    )
  1 4
| 2 5|
||    ||
| 3 6|
|| 710||
( 811)
  912

Notice the need to use Transpose[{v}] in order to convert it to a column matrix. This is needed since in Mathematica, a list can be a row or column, depending on context.

Maple

  restart;
  v1:=<1,2,3>; #column by default
  v2:=<4,5,6>;
  v3:=<7,8,9>;
  v4:=<10,11,12>;
  m:=< <v1|v2>,
       <v3|v4>>;

⌊    ⌋
 1 4
|    |
||2 5 ||
||3 6 ||
|    |
||710 ||
||    ||
⌈811 ⌉
 912



Matlab

  v1=[1,2,3]'; v2=[4,5,6]';
  v3=[7,8,9]'; v4=[10,11,12]';
  m=[v1 v2;v3 v4]
  
  m =
       1     4
       2     5
       3     6
       7    10
       8    11
       9    12

Python

  import numpy as np
  v1=np.array((1,2,3));
  v2=np.array((4,5,6));
  v3=np.array((7,8,9));
  v4=np.array((10,11,12));
  
  r1 =np.hstack([(v1,v2)]).T
  r2 =np.hstack([(v3,v4)]).T
  mat = np.vstack((r1,r2))
  
  Out[211]:
  array([[ 1,  4],
         [ 2,  5],
         [ 3,  6],
         [ 7, 10],
         [ 8, 11],
         [ 9, 12]])

Another way

  v1=np.array([(1,2,3)]).T
  v2=np.array([(4,5,6)]).T
  v3=np.array([(7,8,9)]).T
  v4=np.array([(10,11,12)]).T
  
  mat =np.hstack(( np.vstack((v1,v3)),
                   np.vstack((v2,v4))) )



Fortran

  PROGRAM main
     IMPLICIT none
     INTEGER :: i
     INTEGER , DIMENSION(3) :: v1,v2,v3,v4
     INTEGER , DIMENSION(6,2) :: m
     v1 = [1,2,3]; v2=[4,5,6];
     v3=[7,8,9];   v4=[10,11,12];
  
     m(1:3,1) = v1; m(1:3,2)=v2;
     m(4:,1)=v3;  m(4:,2)=v4;
  
     DO i=1,size(m,1)
        PRINT *, m(i,:)
     END DO
  
  END PROGRAM main
  
  >gfortran -Wall foo.f90
  >./a.out
       </