1.
UP
2.
This document in PDF

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

Nasser M. Abbasi

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

Contents

1 Control systems, Linear systems, transfer functions, state space related problems
 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
  1.1.2 Convert transfer function to state space representation
  1.1.3 Create a state space representation from A,B,C,D and find the step response
  1.1.4 Convert continuous time to discrete time transfer function, compute the gain and phase margins, and show the discrete system step response
  1.1.5 Convert differential equation to transfer functions and to state space
  1.1.6 Convert from continuous transfer function to discrete time transfer function
  1.1.7 Convert a Laplace transfer function to an ordinary differential equation
 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 x0 = 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 Display the impulse response to under, critically, and over damped system
 1.13 View steady state error of 2nd order LTI system with changing undamped natural frequency
  1.13.1 Mathematica
 1.14 Show the use of the inverse Z transform
  1.14.1 example 1
  1.14.2 example 2
 1.15 Find the Z transform of sequence x(n)
  1.15.1 example 1
  1.15.2 example 2
 1.16 Sample a continuous time system
 1.17 Find closed loop transfer function from the open loop transfer function for a unity feedback
 1.18 Compute the Jordan canonical/normal form of a matrix A
 1.19 Solve the continuous-time algebraic Riccati equation
 1.20 Solve the discrete-time algebraic Riccati equation
 1.21 Display impulse response of discrete time system H}z~  along with the impulse response of its continuous time approximation H}s~
 1.22 Find the system type given an open loop transfer function
 1.23 Find the eigenvalues and eigenvectors of a matrix
 1.24 Find the characteristic polynomial of a matrix
 1.25 Verify the Cayley-Hamilton theorem that every matrix is zero of its characteristic polynomial
 1.26 How to check for stability of a continuous system represented as a transfer function and state space
  1.26.1 Checking stability using transfer function poles
  1.26.2 Checking stability using state space A matrix
 1.27 Check continuous system stability in the Lyapunov sense
 1.28 Given a closed loop block diagram, generate the closed loop Z transform and check its stability
 1.29 Determine the state response of a system to only initial conditions in state space
 1.30 Determine the response of a system to only initial conditions in state space
 1.31 Determine the response of a system to step input with nonzero initial conditions
 1.32 Draw the root locus from the open loop transfer function
 1.33 Find eAt   where A  is a matrix
 1.34 Draw root locus for a discrete system
 1.35 Plot the response of the inverted pendulum problem using state space
 1.36 How to build and connect a closed loop control systems and show the response?
  1.36.1 example 1, single input, single output closed loop
 1.37 Compare the effect on the step response of a standard second order system as   changes
 1.38 Plot the dynamic response factor Rd   of a system as a function of r = !!--
     n   for different damping ratios
 1.39 How to find closed loop step response to a plant with a PID controller?
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 Generate a random 2D matrix from uniform (0 to 1) and from normal distributions
 2.6 Generate an n  by m  zero matrix
 2.7 Rotate a matrix by 90 degrees
 2.8 Generate a diagonal matrix with given values on the diagonal
 2.9 Sum elements in a matrix along the diagonal
 2.10 Find the product of elements in a matrix along the diagonal
 2.11 Check if a Matrix is diagonal
 2.12 Find all positions of elements in a Matrix that are larger than some value
 2.13 Replicate a matrix
 2.14 Find the location of the maximum value in a matrix
 2.15 Swap 2 columns in a matrix
 2.16 Join 2 matrices side-by-side and on top of each others
 2.17 Copy the lower triangle to the upper triangle of a matrix to make symmetric matrix
 2.18 extract values from matrix given their index
 2.19 Convert N  by M  matrix to a row of length N M
 2.20 find rows in a matrix based on values in different columns
 2.21 Select entries in one column based on a condition in another column
 2.22 Locate rows in a matrix with column being a string
 2.23 Remove set of rows and columns from a matrix at once
 2.24 Convert list of separated numerical numbers to strings
 2.25 Obtain elements that are common to two vectors
 2.26 Sort each column (on its own) in a matrix
 2.27 Sort each row (on its own) in a matrix
 2.28 Sort a matrix row-wise using first column as key
 2.29 Sort a matrix row-wise using non-first column as key
 2.30 Replace the first nonzero element in each row in a matrix by some value
 2.31 Perform outer product and outer sum between two vector
 2.32 Find the rank and the bases of the Null space for a matrix A
 2.33 Find the singular value decomposition (SVD) of a matrix
 2.34 Solve Ax  = b
 2.35 Find all nonzero elements in a matrix
 2.36 evaluate f }x~  on a vector of values
 2.37 generates equally spaced N points between x1    and x2
 2.38 evaluate and plot a f}x;y~  on 2D grid of coordinates
 2.39 Find determinant of matrix
 2.40 Generate sparse matrix with n  by n  matrix repeated on its diagonal
 2.41 Generate sparse matrix for the tridiagonal representation of the classic second difference operator on 1D grid
 2.42 Generate sparse matrix for the Laplacian differential operator r2u  for 2D grid
 2.43 Generate sparse matrix for the Laplacian differential operator r2u  for 3D grid
 2.44 Generate the adjugate matrix for square matrix
 2.45 Multiply each column by values taken from a row
 2.46 extract submatrix from a larger matrix by removing row/column
 2.47 delete one row from a matrix
 2.48 delete one column from a matrix
 2.49 generate random matrix so that each row adds to 1
 2.50 generate random matrix so that each column adds to 1
 2.51 sum all rows in a matrix
 2.52 sum all columns in a matrix
 2.53 find in which columns values that are not zero
 2.54 How to remove values from one vector that exist in another vector
 2.55 How to find mean of equal sized segments of a vector
 2.56 find first value in column larger than some value and cut matrix from there
 2.57 make copies of each value into matrix into a into a larger matrix
 2.58 How to apply a function to each value in a matrix?
 2.59 How to sum all numbers in a list (vector)?
 2.60 How to find maximum of each row of a matrix?
 2.61 How to find maximum of each column of a matrix?
 2.62 How to add the mean of each column of a matrix from each column?
 2.63 How to add the mean of each row of a matrix from each row?
 2.64 Find the different norms of a vector
 2.65 Check if a matrix is Hermite
 2.66 Obtain the LU decomposition of a matrix
 2.67 Linear convolution of 2 sequences
 2.68 Circular convolution of two sequences
  2.68.1 example 1. The sequences are of equal length
  2.68.2 example 2. The sequences are of unequal length
 2.69 Find the linear convolution of 2 sequences where the origin is located at an arbitrary position in the sequence
 2.70 Visualize a 2D matrix
 2.71 Find the cross correlation between two sequences
 2.72 Find orthonormal vectors that span the range of matrix A
 2.73 Solve Ax  = b  and display the solution
 2.74 Determine if a set of linear equations Ax = b  has a solution and what type of solution
 2.75 Given a set of linear equations automatically generate the matrix A  and vector b  and solve Ax  = b
 2.76 Convert a matrix to row echelon form and to reduced row echelon form
3 Image processing, graphics, signal processing, random numbers
 3.1 Generate and plot one pulse signal of different width and amplitude
 3.2 Generate and plot train of pulses of different width and amplitude
 3.3 Find the discrete Fourier transform of a real sequence of numbers
 3.4 Generate uniform distributed random numbers
  3.4.1 How to generate 5 uniform distributed random numbers from 0 to 1?
  3.4.2 How to generate 5 random numbers from a to b?
  3.4.3 How to generate MATRIX of random numbers from a to b?
 3.5 Determine and plot the Fourier Transform for continuous time function
 3.6 Use FFT to display the power spectrum of the content of an audio wav file
 3.7 Plot the power spectral density of a signal
 3.8 Display spectrum of 2D image
 3.9 Obtain Fourier Series coefficients for a periodic function
 3.10 Obtain the statistical maximum likelihood estimates (MLE) of probability distributions
 3.11 Make a histogram of data sampled from some probability distribution
 3.12 apply a filter on 1D numerical data (a vector)
 3.13 apply an averaging Laplacian filter on 2D numerical data (a matrix)
 3.14 How to compute  10
X   --1---
    k}k  1~
k=1
 3.15 How to find the moving average of a 1D sequence?
 3.16 How to select N random values from a set of numbers?
4 Differential, PDE solving, integration, numerical and analytical solving of equations
 4.1 Generate direction field plot of a first order differential equation
 4.2 Solve the 2-D Laplace PDE  2
r T  x;y  = 0   for a rectangular plate with Dirichlet boundary conditions
 4.3 Solve homogeneous 1st order linear differential equation with constant coefficients and initial conditions
 4.4 Solve homogeneous 2nd order linear differential equation with constant coefficients and initial conditions
 4.5 Solve non-homogeneous 2nd order linear differential equation with constant coefficients and initial conditions
 4.6 Solve homogeneous 2nd order linear differential equation with constant coefficients with boundary values (BVP)
 4.7 Solve the 1-D heat partial differential equation (PDE)
 4.8 Show the effect of different boundary/initial conditions on the solution of the 1-D heat partial differential equation (PDE)
 4.9 Find the particular and homogenous solution to undetermined system of equations
 4.10 Plot the constant energy levels for a nonlinear pendulum
 4.11 Solve numerically the ODE u0000   u = f  using the point collocation method
 4.12 How to numerically solve a set of non-linear equations?
 4.13 Solve 2nd order ODE (Van Der Pol) and generate phase plot
 4.14 How to numerically solve Poisson PDE on 2D using Jacobi iteration method?
 4.15 How to solve BVP second order ODE using finite elements with linear shape functions and using weak form formulation?
 4.16 How to solve Poisson PDE in 2D using finite elements methods using rectangular element?
  4.16.1 Integrating the force vector
  4.16.2 f }x;y~ = 3  cos}4   x~   sin}3  y~
  4.16.3 f }x;y~ = 1
  4.16.4 f }x;y~ = xy
  4.16.5 f }x;y~ = 3  cos}4   x~   sin}3  y~
  4.16.6 Solving for reactangle grid
  4.16.7 f }x;y~ = 1
  4.16.8 f }x;y~ = xy
  4.16.9 f }x;y~ = 3  cos}4   x~   sin}3  y~
 4.17 How to solve Poisson PDE in 2D using finite elements methods using triangle element?
  4.17.1 
f  x;y  = 1
  4.17.2 Program output
  4.17.3 f }x;y~ = 1
 4.18 How to solve wave equation using leapfrog method?
 4.19 Numerically integrate f }x~  on the real line
 4.20 Numerically integrate f }x;y~  in 2D
5 plotting how to, ellipse, 3D
 5.1 Generate a plot using x  and y  data
 5.2 Plot the surface described by f  y;x
 5.3 Find the point of intersection of 3 surfaces
 5.4 How to draw an ellipse?
  5.4.1 draw an ellipse centered at origin and given its semi-major/minor }a;b~
  5.4.2 draw an ellipse centered at }x;y  ) and given its semi-major/minor }a;b~
  5.4.3 draw an ellipse centered at origin and given its semi-major/minor }a;b~  and rotated at angle
 5.5 How to plot a function on 2D using coordinates of grid locations?
 5.6 How to make table of x;y  values and plot them

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.

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.1.1 Create continuous time transfer function given the poles, zeros and gain
  1.1.2 Convert transfer function to state space representation
  1.1.3 Create a state space representation from A,B,C,D and find the step response
  1.1.4 Convert continuous time to discrete time transfer function, compute the gain and phase margins, and show the discrete system step response
  1.1.5 Convert differential equation to transfer functions and to state space
  1.1.6 Convert from continuous transfer function to discrete time transfer function
  1.1.7 Convert a Laplace transfer function to an ordinary differential equation
 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 x0 = 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 Display the impulse response to under, critically, and over damped system
 1.13 View steady state error of 2nd order LTI system with changing undamped natural frequency
  1.13.1 Mathematica
 1.14 Show the use of the inverse Z transform
  1.14.1 example 1
  1.14.2 example 2
 1.15 Find the Z transform of sequence x(n)
  1.15.1 example 1
  1.15.2 example 2
 1.16 Sample a continuous time system
 1.17 Find closed loop transfer function from the open loop transfer function for a unity feedback
 1.18 Compute the Jordan canonical/normal form of a matrix A
 1.19 Solve the continuous-time algebraic Riccati equation
 1.20 Solve the discrete-time algebraic Riccati equation
 1.21 Display impulse response of discrete time system H}z~  along with the impulse response of its continuous time approximation H}s~
 1.22 Find the system type given an open loop transfer function
 1.23 Find the eigenvalues and eigenvectors of a matrix
 1.24 Find the characteristic polynomial of a matrix
 1.25 Verify the Cayley-Hamilton theorem that every matrix is zero of its characteristic polynomial
 1.26 How to check for stability of a continuous system represented as a transfer function and state space
  1.26.1 Checking stability using transfer function poles
  1.26.2 Checking stability using state space A matrix
 1.27 Check continuous system stability in the Lyapunov sense
 1.28 Given a closed loop block diagram, generate the closed loop Z transform and check its stability
 1.29 Determine the state response of a system to only initial conditions in state space
 1.30 Determine the response of a system to only initial conditions in state space
 1.31 Determine the response of a system to step input with nonzero initial conditions
 1.32 Draw the root locus from the open loop transfer function
 1.33 Find eAt   where A  is a matrix
 1.34 Draw root locus for a discrete system
 1.35 Plot the response of the inverted pendulum problem using state space
 1.36 How to build and connect a closed loop control systems and show the response?
  1.36.1 example 1, single input, single output closed loop
 1.37 Compare the effect on the step response of a standard second order system as   changes
 1.38 Plot the dynamic response factor Rd   of a system as a function of r = !--
    !n   for different damping ratios
 1.39 How to find closed loop step response to a plant with a PID controller?

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;
  zeros = [-1 -2];
  poles = [0 -4 -6];
  gain  = 5;
  sys   = zpk(zeros,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 short overview of the system object
  DS:-PrintSystem(sys);

pict

  #see what fields are there in sys
  exports(sys);
   tf, inputcount, outputcount, statecount, sampletime,
   discrete, systemname, inputvariable, outputvariable,
   statevariable, parameters, systemtype, ModulePrint
  
  #reads the transfer function
  tf:=sys[tf];
  
     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 transfer function
              5*s^2+15*s+10
  
  denom(tf[1,1]); #get the denominator of the transfer function
               s*(s^2+10*s+24)

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



Mathematica

  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;
2          3
66  0    1  77
664    25    4 775

  b:=sys:-b;
26 0 37
66   77
64 1 75

  c:=sys:-c;
h      i
  0  5

  d:=sys:-d;
h   i
  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]]

pict



TransferFunctionModel[ss]

pict





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

The plot itself is not shown in these examples.



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}];
  
  Plot[Evaluate@y,{t,0,10},PlotRange->All]

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



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);
  
  size(a);
  size(b);
  size(c);
  size(d);
                                4, 4
                                4, 1
                                1, 4
                                1, 1
  DS:-ResponsePlot(sys, DS:-Step(),duration=10,color=red,legend="step");

pict

1.1.4 Convert continuous time to discrete time transfer function, compute the gain and phase margins, and show the discrete system step response

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~ = s2    s   1

Use sampling period of 0.1 seconds.



Mathematica

  Clear["Global*"];
  tf = TransferFunctionModel[(s + 1)/(s^2 + 1 s + 1), s];
  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]
  

pict





Matlab

  clear all; close all;
  Ts   = 0.1; %sec, sample period
  s    = tf('s');
  sys  = ( s + 1 ) / (  s^2  + s  +  1);
  sysd = c2d(sys,Ts,'zoh'); %convert s to z domain
  tf(sysd)
  step(sysd,10);
  [Gm,Pm,Wcg,Wcp] = margin(sysd)
  set(gcf,'Position',[10,10,410,410]);
  -------
  
     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

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

   2
3d--y    2dy-   y }t~ = u}t~
  dt2    dt



Mathematica

  Remove["Global*"];
  Clear["Global*"];
  
  eq = 3  y''[t]+2 y'[t]+  y[t]==UnitStep[t];
  ss = StateSpaceModel[ {eq},
          {{y'[t],0},{y[t],0}},
          {{UnitStep[t],0}},
          {y[t]},
          t]

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
------------
3s2   2s    1

  sys:=DS:-StateSpace(sys):
  {sys:-a,sys:-b,sys:-c,sys:-d}; #show the state space matrices
8 2             3 2   3                  9
>><  66  0      1   77 66 0 77 h         i h   i>>=
>>  666             777;666   777 ; 1=3   0  ;  0  >> 
: 4   1=3     2=3 5 4 1 5                  ;

1.1.6 Convert from continuous transfer function to discrete time transfer function

Problem: Convert

              1
H }s~ = -------------
        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

  sys = (5s)/(s^2+4s+25);
  
  rhs = Numerator[tf];
  lhs = Denominator[tf];
  rhs = rhs/.m_. s^n_.->  m Derivative[n][u[t]];
  lhs = lhs/.m_. s^n_.->  m Derivative[n][y[t]];
  
  eq= lhs==rhs

  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~ = -------------
        s2    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]];
  yImpulse = Simplify@OutputResponse[sys,DiracDelta[t],t];
  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]
  
  yStep
  
  {1 - E^(-t/10) Cos[(3 Sqrt[11] t)/10]
    - (
    E^(-t/10) Sin[(3 Sqrt[11] t)/10])/(3 Sqrt[11])}
  
  yImpulse
  
  {(10 E^(-t/10) HeavisideTheta[t]
    Sin[(3 Sqrt[11] t)/10])/(3 Sqrt[11])}

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"]);
  #exportplot(FileTools:-JoinPath([currentdir(),"maple_e2.png"]), %);

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

             1
H}s~  = -------------
        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*"];
  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]
  Export["mma_e3.eps", p]

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

       -----25-----
H}s~ = s2    4s   25

Mathematica

  Clear["Global*"];
  sys = TransferFunctionModel[25/(s^2+4 s+25),s];
  TransferFunctionZeros[sys]
  
     {{{}}}
  
  In[228]:= 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
  p
  
  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

H}s~ = -----5s-----
       s2    4s   25



Mathematica

  Clear["Global*"];
  
  tf = TransferFunctionModel[(5 s)/(s^2 + 4 s + 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);

pict

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

pict



1.7 How to check that state space system x0 = 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 x
 0    to x
 f   in finite time. Only the matrix A  and B  are needed to decide on controllability. If the rank of

        2        n  1
[B AB  A B  ::: A   B]

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
         1
B =
          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); #generate controllability matrix
26  0   1    0     2  37
66                    77
66                    77
66  1   0    2     0  77
66                    77
66  0     2   0      10 77
66                    77
64    2  0      10   0  75

  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]
  
  Out[6]= 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');
  
  %this gives
  %r1/(s-p1) + r2/(s-p2) +..... + k
  
  [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);
         7       9        9
s   2    -----    ------   ------
       s    2   2s   6   2s    2

  [op(p0)];
     --7--   --9--- --9---
[s;2;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);
1    e  sT
----2---
   s



1.10 Obtain Inverse Laplace transform of a transfer function

Problem: Obtain the inverse Laplace transform for the function

         s4 --5s3----6s2---9s----30--
H}s~ =  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[%]]

pict



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);
           3e  3t   }  506 cos}3t~    798sin}3t~    299~e  t
Dirac }t~   -----    -------------------------------------
            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

        -------!2n-------
H  }s~ = s2   2  !  s   !2
                 n     n

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

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

pict



1.12 Display the impulse response to under, critically, and over damped system

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

               !2n
H  }s~ =-2-------------2
        s     2  !ns    ! n

in order to illustrate the impulse 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

  Remove["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->#},DiracDelta[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","impulse response for different \[Xi] values"}},
       PlotRange->{{0,12},{-0.8,0.9}},
       GridLines->Automatic,
       GridLinesStyle->Dashed,
       PlotLegend->zValues,
       LegendPosition -> {0.95, -0.5},
       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:10;
  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),~,~] = impulse(num,den,t);
      legendStr(i) = {sprintf('%1.1f',z(i))};
  end
  plot(t,y);
  legend(legendStr);
  title(sprintf
    ('2nd order system impulse 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):
  c:=ColorTools:-GetPalette("Niagara"):
  sys:= (w,zeta)->DS:-TransferFunction(H(w,zeta)):
  zetaValues := [seq(i,i=0.1..1.2,.2)]:
  sol:= map(z->DS:-ImpulseResponse(sys(1,z)),zetaValues):
  seq(plot(sol[i][1,1],t=0..10,
      gridlines=true,color=c[i],
      legend=typeset(zeta,"=",zetaValues[i]),
      legendstyle=[location=right]),
      i=1..nops(sol)):
  plots:-display([%],
    title="impulse response for different damping ratios",
    labels=[typeset(t),typeset(y(t))]);

pict



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

Problem: Given the transfer function

                 2
H  }s~ =-------!n-------
        s2    2  !ns    !2n

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.13.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.14 Show the use of the inverse Z transform

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

1.14.1 example 1

Problem: Given

          z
F }z~ = -----
        z   1

find            1
x [n ] = F  }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.14.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.15 Find the Z transform of sequence x(n)

1.15.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 = f1;1;1;      g  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.15.2 example 2

Find the Z transform for          1  n            n  3
x[n] =  3   u}n~    }0:9~    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.16 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
  set(gcf,'Position',[10,10,320,320]);

pict



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

Problem: Given

L}s~ = ------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?



Mathematica

  Clear["Global*"];
  sys = TransferFunctionModel[s/((s+4)(s+5)),s];
  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]]]]

pict

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.18 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]

pict

  {a,b}=JordanDecomposition[m];
  MatrixForm[b]

pict

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.19 Solve the continuous-time algebraic Riccati equation

Problem: Solve for X  in the Riccati equation

A0X    XA     XBR    1B0X    C0C  = 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]]

pict





Matlab

  %requires the control system toolbox
  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.20 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

A0X    XA     XBR    1B0X    C0C  = 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[%]

pict





Matlab

  clear all; close all;
  s    = tf('s');
  sys  = 1/(s*(s+0.5));
  dsys = c2d(sys,1)
  [A,B,C,D]=dssdata(dsys)
  dare(A,B,C'*C,3)

  dsys =
      0.4261 z + 0.3608
    ----------------------
    z^2 - 1.607 z + 0.6065
  
  Sample time: 1 seconds
  Discrete-time transfer function.
  
  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



1.21 Display impulse response of discrete time system H}z~  along with the impulse response of its continuous time approximation H}s~

Plot the impulse response of H}z~ = z=}z2    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];
  
  (*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"]

  TransferFunctionModel[{
   {{5.6099211033954965 + 1.255589621432896*s}},
   0.5609921103395501 + 1.3862943611198901*s + s^2},
   s,SamplingPeriod -> None,SystemsModelLabels -> None]



Find the impulse response of the continuous time system

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

  (19.767574172152145*HeavisideTheta[t]*
    Sin[0.2837941092083292*t])/E^(0.6931471805599452*t) +
    (1.255589621432896*(0.9999999999999996*Cos[0.2837941092083292*t]*
     HeavisideTheta[t] - 2.4424297688685117*HeavisideTheta[t]*
    Sin[0.2837941092083292*t]))/E^(0.6931471805599452*t)



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*)
  
  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"}}
  ]

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.22 Find the system type given an open loop transfer function

Problem: Find the system type for the following transfer functions

1.
s--1-
s2  s
2.
 s--1-
s3  s2
3.
s-s51

To find the system type, the transfer function is put in the form   P
 k--i}Ps--si~-
sM  j}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 s1}s    1~    and the second system has type 2 since the denominator can be written as s2}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.23 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

pict



  N[Eigenvalues[a]]//MatrixForm

pict



  N[Eigenvectors[a]]//MatrixForm

pict





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.24 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.25 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]
  Out[15]= -4-3 x+x^2
  (-4 IdentityMatrix[n] - 3 a +
      MatrixPower[a,2])//MatrixForm

pict

Another way to do the above is as follows

  Remove["Global*"]
  a  = {{1,2},{3,2}};
  p  = CharacteristicPolynomial[a,x]
  Out[35]= -4-3 x+x^2
  cl = CoefficientList[p,x];
  Sum[MatrixPower[a,j-1] cl[[j]],
       {j,1,Length[cl]}]
  Out[37]= {{0,0},{0,0}}

Matlab MATLAB has a build-in function polyvalm() to do this more easily than in Mathematica. Although the above shown method in Mathematica can be easily made into a 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.26 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.26.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.26.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.27 Check continuous system stability in the Lyapunov sense

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

     20    1   0 3
     666           777
A =  660    0   1 77
     66           77
     4  1     2    35

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 AT P    PA  =   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}})//MatrixForm

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

f6:18272; 1:1149;0:202375g





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

2                         3
66 6:18272045921436      0:0i77
66                         77
66 1:11490451203192      0:0i77
66                         77
640:202375028753723     0:0i75



1.28 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.325556 I,0.543991 +0.325556 I}}}



  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
  set(gcf,'Position',[10,10,320,320]);

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.29 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}]

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.30 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.31 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.32 Draw the root locus from the open loop transfer function

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

                s2   2s    4
L}s~ = ----------------2-----------
       s}s   4~}s    6~}s    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.33 Find eAt   where A  is a matrix



Mathematica

  ClearAll["Global*"]
  mat={{0,1},
       {-2,-3}}
  MatrixExp[mat t];
  MatrixForm[%]

pict

Now verify the above result by solving for  At
e   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 jA      Ij

  m = mat-lambda IdentityMatrix[Length[mat]]
  Out[9]= {{-lambda,  1        },
           {-2,       -3-lambda}}
  Det[m]
  Out[14]= 2+3 lambda+lambda^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}]
  
  Out[20]= {b0->E^(-2 t) (-1+2 E^t),
            b1->-E^(-2 t)+E^-t}
  (*Now find e^At*)
  b0=b0/.sol[[1]]
  Out[24]= E^(-2 t) (-1+2 E^t)
  
  b1=b1/.sol[[2]]
  Out[25]= -E^(-2 t)+E^-t
  
  b0 IdentityMatrix[Length[mat]]+b1 mat;
  Simplify[%]
  MatrixForm[%]

pict

The answer above 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.34 Draw root locus for a discrete system

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

sys =  -52s----1---
      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.35 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}};
  MatrixForm[a]

pict



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

pict



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

pict



  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.41561 t) (0.000115693 E^(3.2078 t) UnitStep[t]-
   0.000231385 E^(6.41561 t) UnitStep[t]+
   0.000115693 E^(9.62341 t) UnitStep[t]+
   0.0238095 E^(6.41561 t) t^2 UnitStep[t]),
  
   E^(-6.41561 t) (-0.00242954 E^(3.2078 t) UnitStep[t]+
   0.00485909 E^(6.41561 t) UnitStep[t]-
   0.00242954 E^(9.62341 t) UnitStep[t])}



  (*The following obtained the response
  as interpolatingFunction
  since the full range of time
  domain is now specified*)
  
  sys=sys/.{g->9.8,m->1,M->20,L->1};
  y=OutputResponse[sys,UnitStep[t],{t,0,3}]

   {InterpolatingFunction[{{0.,3.}},<>][t],
    InterpolatingFunction[{{0.,3.}},<>][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.36 How to build and connect a closed loop control systems and show the response?

1.36.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



1.37 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~ = -2-------------2
        s     2  !ns    ! n

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

                 2
-Y }s~ = ------!-n-------
U  }s~    s2   2  !ns    !2n

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

y00}t~    2  !n  y0}t~    !2 y }t~ = !2 u }t~
                      n          n

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

2  03  2             32  3   2   3
66x 177= 66  0      1   7766x177    66 0 77u }t~
66x0 77  66   !2n    2  !n 7766x277   66!2n77
4  25  4             54  5   4   5

For a step input, u}t~ = 1   , we obtain

266x01377 = 266 0      1   377266x1 377   266 0 377
66x077   66  !2     2  !n 7766x2 77  66!2 77
4 25   4   n        54   5  4  n5

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

pict



1.38 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.39 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]

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

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

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

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 t1
  implicit none
  
  integer, parameter :: N=3
  real(8), 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(8), parameter, DIMENSION(N, 1) :: &
    x=DBLE(reshape([1.0, 2.0, 3.0], shape(x)))
  
  real(8), DIMENSION(N, 1) :: result
  integer, DIMENSION(2)    :: siz
  
  result = MATMUL(A,x)
  
  Write(*,*) result
  siz = SHAPE(result)
  print *, "number of rows = ",siz(1)
  print *, "number of columns = ",siz(2)
  
  end program
  
  #compile and run
  
  $ gfortran -std=f2003 -Wextra -Wall -pedantic
    -funroll-loops -ftree-vectorize -march=native
    -Wsurprising -Wconversion
    t1.f90  /usr/lib/liblapack.a /usr/lib/libblas.a
  $ ./a.exe
   14.000000000000000 32.000000000000000 50.000000000000000
   number of rows =            3
   number of columns =            1



Julia

  A=[1 2 3;4 5 6;7 8 9]
  x=[1;2;3]
  c=A*x
  
  3-element Array{Int64,1}:
   14
   32
   50

Maple

  A:=Matrix([[1,2,3],[4,5,6],[7,8,9]]):
  x:=Vector([1,2,3]): #default is column vector
  c:=A.x;
2  3
661477
663277
66  77
645075



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..]]);
h       i
 129934



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
  >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,..]>>;
26       37
661  2 3 77
664  5 6 77
66       77
66909192 77
667  8 9 77
4       5

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,..]
            ] ]);
2       3
661  2 3 77
664  5 6 77
66       77
666909192 777
667  8 9 77
4       5



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..]>>;
2619023 37
66      77
6649156 77
6679289 77
4      5

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..] ]]);
2      3
6619023 77
6649156 77
66      77
6479289 75



2.5 Generate a random 2D matrix from uniform (0 to 1) and from normal distributions



Mathematica

  (*from uniform over 0,1*)
  
  mat = RandomReal[
         UniformDistribution[{0,1}],{3,4}]
  
   {{0.100843,0.356115,0.700317,0.657852},
    {0.238019,0.59598,0.523526,0.576362},
    {0.339828,0.32922,0.0632487,0.815892}}
  
  (*from normal, zero mean, 1 std*)
  
  mat=RandomReal[NormalDistribution[0,1],{3,4}]
  
    {{-0.226424,1.19553,0.601433,1.28032},
    {1.66596,0.176225,-0.619644,0.371884},
    {0.895414,-0.394081,0.153507,-1.74643}}

Matlab

  mat=rand(3,4)
  
  mat =
   0.6787    0.3922    0.7060    0.0462
   0.7577    0.6555    0.0318    0.0971
   0.7431    0.1712    0.2769    0.8235
  
  mat=randn(3,4)
  
  mat =
   0.3252   -1.7115    0.3192   -0.0301
  -0.7549   -0.1022    0.3129   -0.1649
   1.3703   -0.2414   -0.8649    0.6277



Maple

  A:=ArrayTools:-RandomArray(3,4,distribution = uniform);
260:743132468124916179   0:933993247757550549    0:6557406991565868370:915735525189067090     37
66                                                                                          77
660:757740130578333448   0:849129305868777107    0:9594924263929029970:421761282626274991     77
660:6787351548577734710:03571167857418955370:7922073295595544180:141886338627215336         77
4                                                                                          5

  interface(displayprecision=5):
  A:=ArrayTools:-RandomArray(3,4,distribution = uniform);

pict



Fortran

  program t5
    implicit none
    INTEGER(4) ::i=0,j
    REAL ::A(3,3),mean,sd,pi,temp
  
    CALL RANDOM_SEED(i)
    CALL RANDOM_NUMBER(A)
  
    print *, 'from uniform distribution'
    do i=LBOUND(A, 1),UBOUND(A, 1)
       print *,A(i,:)
    end do
  
    !this block below uses the logic as explained in
    !http://rosettacode.org/wiki/Random_numbers#Fortran
    !
    mean=0.0
    sd=1.0
    pi = 4.0*ATAN(1.0)
  
    ! Now convert to normal distribution
      DO i = 1, SIZE(A,1)-1, 2
        DO j = 1, SIZE(A,2)-1, 2
           temp = sd * SQRT(-2.0*LOG(A(i,j))) * COS(2*pi*A(i+1,j+1)) + mean
           A(i+1,j+1) = sd * SQRT(-2.0*LOG(A(i,j))) * SIN(2*pi*A(i+1,j+1)) + mean
           A(i,j) = temp
        END DO
     END DO
  
   print *, 'from normal distribution'
   do i=LBOUND(A, 1),UBOUND(A, 1)
        print *,A(i,:)
   end do
  
  end program t5
  
  >gfortran -Wall -std=f2008 t5.f90
  >./a.out
   from uniform distribution
    0.99755955      0.74792767      7.37542510E-02
    0.56682467      0.36739087      5.35517931E-03
    0.96591532      0.48063689      0.34708124
  
   from normal distribution
   -4.70122509E-02  0.74792767      7.37542510E-02
    0.56682467      5.17370142E-02  5.35517931E-03
    0.96591532      0.48063689      0.34708124

Did not find a build-in support for random numbers from normal distribution.

2.6 Generate an n  by m  zero matrix



Mathematica

  mat=Table[0,{3},{4}]
  
      {{0,0,0,0},
       {0,0,0,0},
       {0,0,0,0}}
  

Matlab

  A=zeros(3,4)
  
  A =
       0     0     0     0
       0     0     0     0
       0     0     0     0



Maple

  LinearAlgebra:-ZeroMatrix(3,4);
2    3
66000077
66000077
66    77
64000075



2.7 Rotate a matrix by 90 degrees



Mathematica

  mat = {{1, 2, 3},
         {4, 5, 6},
         {7, 8, 9}};
  mat = Reverse[Transpose[mat]]
      {{3,6,9},
       {2,5,8},
       {1,4,7}}

Matlab

  A=[1 2 3;
     4 5 6;
     7 8 9]
  A=rot90(A)
       3     6     9
       2     5     8
       1     4     7



Maple

  A:=Matrix([ [1,2,3],[4,5,6],[7,8,9]]);
  A:=ArrayTools:-FlipDimension(A,2);
26321 37
66    77
66654 77
66987 77
4    5

  A:=A^%T;
2    3
66369 77
66258 77
66    77
64147 75



2.8 Generate a diagonal matrix with given values on the diagonal

Problem: generate diagonal matrix with 2;4;6;8   on the diagonal.



Mathematica

  DiagonalMatrix[2 Range[1, 4]]
       {{2,0,0,0},
        {0,4,0,0},
        {0,0,6,0},
        {0,0,0,8}}

Matlab

  diag(2*(1:4))
       2     0     0     0
       0     4     0     0
       0     0     6     0
       0     0     0     8



Maple

  A:=LinearAlgebra:-DiagonalMatrix([2,4,6,8]);
  #or
  A:=LinearAlgebra:-DiagonalMatrix([seq(i,i=2..8,2)]);
2    3
66200077
66    77
66040077
66006077
66    77
64000875



2.9 Sum elements in a matrix along the diagonal



Mathematica

  mat = {{1, 2, 3},
         {4, 5, 6},
         {7, 8, 9}}
  Tr[mat]
  Out[45]= 15

Matlab

  A=[1 2 3;
     4 5 6;
     7 8 9]
  sum(diag(A))
  
  ans =
      15



Maple

  A:=Matrix([ [1,2,3],[4,5,6],[7,8,9]]);
  MTM:-diag(A);
  add(x,x in %);
  
       15

Another ways to obtain the diagonal of a matrix is

  LinearAlgebra:-Diagonal(A);
  Student:-LinearAlgebra:-Diagonal(A);



2.10 Find the product of elements in a matrix along the diagonal



Mathematica

  mat= {{1, 2, 3},
        {4, 5, 6},
        {7, 8, 9}}
  
  Apply[Times,Diagonal[mat]]
  
  Out[49]= 45

Matlab

  A=[1 2 3;
     4 5 6;
     7 8 9]
  prod(diag(A))
  
  ans =
      45



Maple

  A:=Matrix([ [1,2,3],[4,5,6],[7,8,9]]);
  mul(x,x in LinearAlgebra:-Diagonal(A));

45



2.11 Check if a Matrix is diagonal

A diagonal matrix is one which has only zero elements off the diagonal. The Mathematica code was contributed by Jon McLoone.



Mathematica

  diagonalQ[m_List]/;ArrayDepth[m]===2&&Equal@@Dimensions[m]:=
    And@@Flatten[MapIndexed[#1===0||Equal@@#2&,m,{2}]];
  
  diagonalQ[m_]:=Return[False];
  matA = {{1, 2},
          {2, 4}};
  
  matB = {{1, 0},
          {0, 2}};
  
  matC = {{1, 0, 2},
          {0, 2, 4}};
  
  diagonalQ[matA]
  diagonalQ[matB]
  diagonalQ[matC]

  Out[59]= False
  Out[60]= True
  Out[61]= False



Maple

  A:=Matrix([[1,0],[0,2]]);
  Student:-NumericalAnalysis:-IsMatrixShape(A,'diagonal');
    true



2.12 Find all positions of elements in a Matrix that are larger than some value

The problem is to find locations or positions of all elements in a matrix that are larger or equal than some numerical value such as 2   in this example.



Mathematica

  mat = {{1, 2,3},
         {2,1,-5},
         {6,0,0}};
  
  p = Position[mat, _?(#1 >= 2&)]
  
     {{1,2},{1,3},{2,1},{3,1}}

Matlab

  A=[1 2 3;
     2 1 -5;
     6 0 0]
  [I,J]=find(A>=2)
  I =
       2
       3
       1
       1
  J =
       1
       1
       2
       3



Maple

  A:=Matrix( [ [1,2,3],[2,1,-5],[6,0,0]]);
  seq(seq(if(A[i,j]>=2,[i,j],NULL),i=1..3),j=1..3);
  
  Gives
  
   [2, 1], [3, 1], [1, 2], [1, 3]



2.13 Replicate a matrix

Given Matrix

       1     2
       3     4

Generate a new matrix of size 2   by 3   where each element of the new matrix is the above matrix. Hence the new matrix will look like

       1     2     1     2     1     2
       3     4     3     4     3     4
       1     2     1     2     1     2
       3     4     3     4     3     4

In Matlab, repmat() is used. In Mathematica, a Table command is used, followed by ArrayFlatten[]



Mathematica

  mat = {{1,2},
         {3,4}}
  
  r = Table[mat,{2},{3}];
  ArrayFlatten[r]
  
    {{1,2,1,2,1,2},
     {3,4,3,4,3,4},
     {1,2,1,2,1,2},
     {3,4,3,4,3,4}}
  

Matlab

  A=[1 2;3 4]
  repmat(A,2,3)
  
  ans =
       1     2     1     2     1     2
       3     4     3     4     3     4
       1     2     1     2     1     2
       3     4     3     4     3     4



2.14 Find the location of the maximum value in a matrix



Mathematica

  mat ={{2,4,-3},
        {7,9,5},
        {0,5.4,9}}
  
  m = Max[mat];  (*this gives values *)
  Position[mat,m](*this find locations*)
  
  Out[142]= {{2,2},{3,3}}

Matlab

  h = [2, 4, -3;
       7, 9, 5;
       0, 5.4, 9]
  
  m = max(h(:));
  
  [r,c]=find(h==m)
  r =
       2
       3
  c =
       2
       3



2.15 Swap 2 columns in a matrix

Give a matrix

       1     2
       3     4
       5     6

How to change it so that the second column becomes the first, and the first becomes the second? so that the result become

       2     1
       4     3
       6     5


Mathematica

  a={{1,2},
     {3,4},
     {5,6}}
  Reverse[a,2]
  Out[29]= {{2, 1},
            {4, 3},
            {6, 5}}

Matlab

  a =[1 2;
      3 4;
      5 6]
  [a(:,2) a(:,1)]
  
       2     1
       4     3
       6     5



2.16 Join 2 matrices side-by-side and on top of each others



Mathematica In Mathematica, to join 2 matrices side-by-side, use Join with '2' as the third argument. To join them one on top of the other, use '1' as the third argument

  a = RandomInteger[10, {3, 3}]
  Out[146]= {{8,5, 1},
             {8,10,8},
             {6,0, 8}}
  
  b = RandomInteger[10, {3, 3}]
  Out[147]= {{9, 0,1},
             {10,2,7},
             {8, 0,8}}
  
  (*join side-by-side as in Matlab [A B]*)
  Join[a, b, 2]
  Out[149]= {{8,5, 1, 9,0,1},
             {8,10,8,10,2,7},
             {6,0, 8, 8,0,8}}
  
  (*join vertically as in Matlab [A;B]*)
  Join[a, b, 1]
  Out[150]= {{8, 5 , 1},
             {8, 10, 8},
             {6, 0 , 8},
             {9, 0 , 1},
             {10,2 , 7},
             {8, 0 , 8}}

Matlab

  A=rand(3,3)
  B=rand(3,3)
  [A B]
  
   0.5472  0.2575  0.8143 0.3500 0.6160 0.8308
   0.1386  0.8407  0.2435 0.1966 0.4733 0.5853
   0.1493  0.2543  0.9293 0.2511 0.3517 0.5497
  
  [A;B]
  
      0.5472    0.2575    0.8143
      0.1386    0.8407    0.2435
      0.1493    0.2543    0.9293
      0.3500    0.6160    0.8308
      0.1966    0.4733    0.5853
      0.2511    0.3517    0.5497
  



2.17 Copy the lower triangle to the upper triangle of a matrix to make symmetric matrix

Question posted on the net

  please help me with simple to apply function that will construct
  symmetric matrix from given just a half matrix with diagonal.
  Eg:
  
  From:
  1  0  0  0
  2  3  0  0
  4  9  5  0
  2  2  3  4
  
  To give:
  
  1  2  4  2
  2  3  9  2
  4  9  5  3
  2  2  3  4

Many answers were given, below is my answer, and I also show how to do it in Matlab



Mathematica

  a = {{1, 0, 0, 0},
       {2, 3, 0, 0},
       {4, 9, 5, 0},
       {2, 2, 3, 4}};
  
  Transpose[LowerTriangularize[a, -1]] + a
  
   {{1, 2, 4, 2},
    {2, 3, 9, 2},
    {4, 9, 5, 3},
    {2, 2, 3, 4}
   }

Matlab

  a =
       1     0     0     0
       2     3     0     0
       4     9     5     0
       2     2     3     4
  
  a+tril(a,-1)'
  
  ans =
       1     2     4     2
       2     3     9     2
       4     9     5     3
       2     2     3     4



2.18 extract values from matrix given their index

Given a matrix A, and list of locations within the matrix, where each location is given by i;j  entry, find the value in the matrix at these locations

Example, given

  A={{1,2,3},
     {4,5,6},
     {7,8,9}}

obtain the entries at 1;1   and 3;3   which will be 1   and 9   in this example.



Mathematica

  mat={{1,2,3},
       {4,5,6},
       {7,8,9}}
  
  pos = { {1,1},{3,3}};
  Map[ mat[[Sequence@@ # ]] & , pos ]
  
     {1,9}

Another method (same really as above, but using Part explicit)

  Map[Part[A, Sequence @@ #] &, pos]
  {1,9}

Matlab

  A=[1 2 3;
     4 5 6;
     7 8 9];
  I=[1 3];  % setup the I,J indices
  J=[1 3];
  A(sub2ind(size(A),I,J))
  
  ans =
  
       1     9
  



2.19 Convert N  by M  matrix to a row of length N M

Given

  a=[1 2 3
     4 5 6
     7 8 9]

covert the matrix to one vector

  [1 2 3 4 5 6 7 8 9]


Mathematica

  a={{1,2,3},
     {4,5,6},
     {7,8,9}}
  
  Flatten[a]
  
  Out[17]= {1,2,3,4,5,6,7,8,9}

Matlab

  a=[1 2 3;
     4 5 6;
     7 8 9]
  
  a=a(:);
  a'
  >>
   1 4 7 2 5 8 3 6 9



2.20 find rows in a matrix based on values in different columns

Example, given

       1     9     0    10
       5     6     7     8
       3     9     2    10

Select rows which has 9   in the second column and 10   in the last column. Hence the result will be the first and third rows only

        1     9     0    10
        3     9     2    10


Mathematica

  mat={{1,9,0,10},
        {5,6,7,8},
        {3,9,2,10}};
  
  p = Position[mat[[All,{2,4}]],x_/;x==={9,10}]
  Out[184]= {{1},{3}}
  
  
  Extract[mat,p]
  Out[185]= {{1,9,0,10},
             {3,9,2,10}}

Matlab

  A=[1 9 0 10;
     5 6 7 8;
     3 9 2 10];
  
  r=bsxfun(@eq,A(:,[2 4]),[9 10]);
  ind=all(r,2);
  A(ind,:)
  
  ans =
       1     9     0    10
       3     9     2    10



2.21 Select entries in one column based on a condition in another column

Given

  A=[4  3
     6  4  ---->
     7  6  ---->
     2  1
     1  3
     9  2
     2  5  ---->
     1  2
     ]

Select elements in the first column only which has corresponding element in the second column greater than 3, hence the result will be

       [6 7 2]


Mathematica

  mat = {{4, 3},
         {6, 4},
         {7, 6},
         {2, 1},
         {1, 3},
         {9, 2},
         {2, 5},
         {1, 2}};
  
  r = Select[mat, #[[2]] > 3 &];
  r[[All, 1]]
  
        {6, 7, 2}

another way is to find the index using Position and then use Extract

  loc = Position[mat, {x_, y_} /; y > 3]
  r = Extract[mat, loc];
  r[[All, 1]]
  
      {6, 7, 2}

another way is to use Cases[]. This is the shortest way

  Cases[mat, {x_, y_} /; y > 3 :> x]
  
      {6, 7, 2}

Matlab

  A=[4, 3;
     6, 4;
     7, 6;
     2, 1;
     1, 3;
     9, 2;
     2, 5;
     1, 2];
  
  A(A(:,2)>3,1)
  
       6
       7
       2



2.22 Locate rows in a matrix with column being a string

The problem is to select rows in a matrix based on string value in the first column. Then sum the total in the corresponding entries in the second column. Given. For example, given

  mat = {'foobar', 77;
         'faabar', 81;
         'foobur', 22;
         'faabaa', 8;
         'faabian', 88;
         'foobar', 27;
         'fiijii', 52};

and given list

  {'foo', 'faa'}

The problem is to select rows in which the string in the list is part of the string in the first column in the matrix, and then sum the total in the second column for each row found. Hence the result of the above should be

      'foo'    'faa'
      [126]    [ 177]


Mathematica

  
  mat = {{"foobar", 77},
     {"faabar", 81},
     {"foobur", 22},
     {"faabaa", 8},
     {"faabian", 88},
     {"foobar", 27},
     {"fiijii", 52}};
  lst = {"foo", "faa"};
  
  foo[mat_, lst_] :=
    Select[mat, StringMatchQ[lst,
        StringTake[#[[1]], 3]] &];
  
  r = Map[foo[mat, #] &, lst];
  
  MapThread[{#1,
   Total[#2[[All, 2]]]}&,{lst, r}]
  
  {{"foo", 126}, {"faa", 177}}

Matlab

  A = {'foobar', 77;
     'faabar', 81;
     'foobur', 22;
     'faabaa', 8;
     'faabian', 88;
     'foobar', 27;
     'fiijii', 52};
  lst= {'foo', 'faa'};
  f=@(x) [x sum(cell2mat(A(strncmp(A(:,1),x,3),2)))];
  r=arrayfun(@(i) f(lst(i)),1:2,'UniformOutput',false)
  reshape([r{:}],2,2)
  ans =
      'foo'    'faa'
      [126]    [177]

But notice that in Matlab, the answer is a cellarray. To access the numbers above

  r{:}
  ans =
      'foo'    [126]
  ans =
      'faa'    [177]
  
  r{1}{2}
  ans =
     126
  
  r{2}{2}
  ans =
      177



2.23 Remove set of rows and columns from a matrix at once

Given: square matrix, and list which represents the index of rows to be removed, and it also represents at the same time the index of the columns to be removed (it is square matrix, so only one list is needed).

output: the square matrix, with BOTH the rows and the columns in the list removed.

Assume valid list of indices.

This is an example: remove the second and fourth rows and the second and fourth columns from a square matrix.

pict

I asked this question at SO, and more methods are shown there at HTML



Mathematica Three methods are shown.

method 1: (credit for the idea to Mike Honeychurch at stackoverflow). It turns out it is easier to work with what we want to keep instead of what we want to delete so that Part[] can be used directly.

Hence, given a list of row numbers to remove, such as

  pos = {2, 4};

Start by generating list of the rows and columns to keep by using the command Complement[], followed by using Part[]

  a = {{0, 5, 2, 3, 1, 0},
       {4, 3, 2, 5, 1, 3},
       {4, 1, 3, 5, 3, 2},
       {4, 4, 1, 1, 1, 5},
       {3, 4, 4, 5, 3, 3},
       {5, 1, 4, 5, 2, 0}
     };
  pos = {2, 4};
  keep = Complement[Range[Length[a]], pos];
  a[[keep, keep]]
  
  {{0, 2, 1, 0},
   {4, 3, 3, 2},
   {3, 4, 3, 3},
   {5, 4, 2, 0}
  }

method 2: (due to Mr Wizard at stackoverflow)

  ReplacePart[a, {{2},{4},{_, 2},{_, 4}}
        :> Sequence[]]
  
  {{0, 2, 1, 0},
   {4, 3, 3, 2},
   {3, 4, 3, 3},
   {5, 4, 2, 0}
  }

method 3: (me)

use Pick. This works similar to Fortran pack(). Using a mask matrix, we set the entry in the mask to False for those elements we want removed. Hence this method is just a matter of making a mask matrix and then using it in the Pick[] command.

  {nRow,nCol} = Dimensions[a];
  mask = Table[True,{nRow},{nCol}];
  pos  = {2,4};
  mask[[All,pos]]=False;
  mask[[pos,All]]=False;
  r=Pick[a,mask]/.{}->Sequence[]
  
  Out[39]= {{0,2,1,0},
            {4,3,3,2},
            {3,4,3,3},
            {5,4,2,0}}

Matlab

  a = [0, 5, 2, 3, 1, 0;
       4, 3, 2, 5, 1, 3;
       4, 1, 3, 5, 3, 2;
       4, 4, 1, 1, 1, 5;
       3, 4, 4, 5, 3, 3;
       5, 1, 4, 5, 2, 0
     ]
  
  list=[2 4];
  a(list,:)=[];
  a(:,list)=[];
  
  a =
       0     2     1     0
       4     3     3     2
       3     4     3     3
       5     4     2     0



2.24 Convert list of separated numerical numbers to strings

Problem: given a list of numbers such as

  {{1, 2},
   {5, 3, 7},
   {4}
   }

convert the above to list of strings

  {"12",
   "537",
   "4"}


Mathematica

  a={{1,2},{5,3,7},{4}};
  Map[ToString[#]&,a,{2}];
  Map[StringJoin[#]&,%]
  
  Out[402]//FullForm= List["12","537","4"]

Matlab

  a={[1,2],[5,3,7],[4]};
  r=arrayfun(@(i) ...
     cellstr(num2str(a{i})),1:length(a));
  r'
  
  ans =
      '1  2'
      '5  3  7'
      '4'

This answer below is due to Bruno Luong posted at the Matlab newsgroup

  a={[1,2],[5,3,7],[4]};
  map='0':'9';
  cellfun(@(x) map(x+1), a, 'uni', 0)
  
  ans =
  
      '12'    '537'    '4'
  



2.25 Obtain elements that are common to two vectors

Given vector or list d = [  9;1;3;  3;50; 7;19 ];  t = [0;7;2;50 ]   , find the common elements.



Mathematica

  d = {-9,1,3,-3,50,7,19};
  t = {0,7,2,50};
  Intersection[d,t]
  
  Out[412]= {7,50}

Matlab

  d=[-9 1 3 -3 50 7 19];
  t =[0 7 2 50];
  intersect(t,d)
  
  ans =
       7    50



2.26 Sort each column (on its own) in a matrix

Given

       4     2     5
       2     7     9
      10     1     2

Sort each column on its own, so that the result is

       2     1     2
       4     2     5
      10     7     9

In Matlab, the sort command is used. But in the Mathematica, the Sort command is the same the Matlab's sortrows() command, hence it can't be used as is. Map is used with Sort to accomplish this.



Mathematica

  mat={{4,2,5},
       {2,7,9},
       {10,1,2}};
  
  Map[Sort[#]&,Transpose[mat]];
  mat=Transpose[%]
  
  Out[176]= {{2, 1, 2},
             {4, 2, 5},
             {10,7, 9}}

Matlab

  A=[4  2 5;
     2  7 9;
     10 1 2];
  sort(A,1)
  
  ans =
       2     1     2
       4     2     5
      10     7     9



2.27 Sort each row (on its own) in a matrix

Given

       4     2     5
       2     7     9
      10     1     2

Sort each row on its own, so that the result is

       2     4     5
       2     7     9
       1     2     10


Mathematica

  mat={{4, 2,5},
       {2, 7,9},
       {10,1,2}};
  Map[Sort[#]&,mat]
  
  Out[180]= {{2, 4, 5},
             {2, 7, 9},
             {1, 2, 10}}
  

Matlab

  A=[4  2 5;
     2  7 9;
     10 1 2];
  sort(A,2)
  
  ans =
       2     4     5
       2     7     9
       1     2    10



2.28 Sort a matrix row-wise using first column as key

Given

       4     2     5
       2     7     9
      10     1     2

Sort the matrix row-wise using first column as key so that the result is