1.
UP
2.
This document in PDF

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

Nasser M. Abbasi

July 18, 2014

Contents

1 Obtain the step response of an LTI from its transfer function
2 plot the impulse and step responses of a system from its transfer function
3 Obtain the response of a transfer function for an arbitrary input
4 Obtain the poles and zeros of a transfer function
5 Obtain the continuous time transfer function given the poles, zeros and gain
6 Convert transfer function to state space representation
7 Create a state space representation from A,B,C,D and find the step response
8 Generate Bode plot of a transfer function
9 How to Convert continuous time to discrete time transfer function, compute the gain and phase margins, and show the discrete system step response
10 Convert transfer function to state space
11 Obtain partial-fraction expansion
12 Obtain Laplace transform for a piecewise functions
13 Obtain Inverse Laplace transform of a transfer function
14 Display the response to a unit step of an under, critically, and over damped system
15 Display the impulse response to under, critically, and over damped system
16 Convert differential equation to transfer functions and state space
17 View the steady state error of 2nd order LTI system as the undamped natural frequency changes
 17.1 Mathematica
  17.1.1 ramp input
  17.1.2 step input
18 Convert a system from continuous time Laplace transfer function to discrete time Z transfer function
19 Show the use of the inverse Z transform
 19.1 example 1
 19.2 example 2
20 Find the Z transform of sequence x(n)
 20.1 example 1
 20.2 example 2
21 Sample a continuous time system
22 Find closed loop transfer function from the open loop transfer function for a unity feedback
23 Linear convolution of 2 sequences
24 Circular convolution of two sequences
 24.1 example 1. The sequences are of equal length
 24.2 example 2. The sequences are of unequal length
25 Compute the Jordan canonical/normal form of a matrix A
26 Plot the power spectral density of a signal
27 Solve the continuous-time algebraic Riccati equation
28 Solve the discrete-time algebraic Riccati equation
29 Display impulse response of discrete time system H  (z )  along with the impulse response of its continuous time approximation H (s)
30 Find the system type given an open loop transfer function
31 Find the eigenvalues and eigenvectors of a matrix
32 Find the characteristic polynomial of a matrix
33 Verify the Cayley-Hamilton theorem that every matrix is zero of its characteristic polynomial
34 Solve Ax  = b  and display the solution
35 Determine if a set of linear equations Ax  = b  has a solution and what type of solution
36 Given a set of linear equations automatically generate the matrix A  and vector b  and solve Ax  = b
37 Convert a matrix to row echelon form and to reduced row echelon form
38 How to check for stability of a continuous system represented as a transfer function and state space
 38.1 Checking stability using transfer function poles
 38.2 Checking stability using state space A matrix
39 Check continuous system stability in the Lyapunov sense
40 Given a closed loop block diagram, generate the closed loop Z transform and check its stability
41 Check if a matrix is Hermite
42 Determine the state response of a system to only initial conditions in state space
43 Determine the response of a system to only initial conditions in state space
44 Determine the response of a system to step input with nonzero initial conditions
45 Draw the root locus from the open loop transfer function
46 Find the cross correlation between two sequences
47 Find the different norms of a vector
48 Find orthonormal vectors that span the range of matrix A
49 Plot the surface described by f (y, x)
50 Find the point of intersection of 3 surfaces
51 Obtain the LU decomposition of a matrix
52 Use FFT to display the power spectrum of the content of an audio wav file
53 Solve homogeneous 1st order linear differential equation with constant coefficients and initial conditions
54 Solve homogeneous 2nd order linear differential equation with constant coefficients and initial conditions
55 Solve non-homogeneous 2nd order linear differential equation with constant coefficients and initial conditions
56 Solve homogeneous 2nd order linear differential equation with constant coefficients with boundary values (BVP)
57 Solve the 1-D heat partial differential equation (PDE)
58 Show the effect of different boundary/initial conditions on the solution of the 1-D heat partial differential equation (PDE)
59 Generate uniform distributed random numbers
 59.1 How to generate 5 uniform distributed random numbers from 0 to 1?
 59.2 How to generate 5 random numbers from a to b?
 59.3 How to generate MATRIX of random numbers from a to b?
60 Determine and plot the Fourier Transform for a continuous time function
61 Solve the 2-D Laplace PDE   2
∇  T (x,y) = 0  for a rectangular plate with Dirichlet boundary conditions
62 Generate and plot one pulse signal of different width and amplitude
63 Generate and plot train of pulses of different width and amplitude
64 Find the discrete Fourier transform of a real sequence of numbers
65 Find eAt  where A  is a matrix
66 How to do some common operations on vectors and matrices?
 66.1 Multiply matrix with a vector
 66.2 Insert a number at specific position in a vector or list
 66.3 Insert a row into a matrix
 66.4 Insert a column into a matrix
 66.5 Generate a random 2D matrix from uniform (0 to 1) and from normal distributions
 66.6 Generate an n  by m  zero matrix
 66.7 Rotate a matrix by 90 degrees
 66.8 Generate a diagonal matrix with given values on the diagonal
 66.9 Sum elements in a matrix along the diagonal
 66.10 Find the product of elements in a matrix along the diagonal
 66.11 Check if a Matrix is diagonal
 66.12 Find all positions of elements in a Matrix that are larger than some value
 66.13 Replicate a matrix
 66.14 Find the location of the maximum value in a matrix
 66.15 Swap 2 columns in a matrix
 66.16 Join 2 matrices side-by-side and on top of each others
 66.17 Copy the lower triangle to the upper triangle of a matrix to make symmetric matrix
 66.18 extract values from matrix given their index
 66.19 Convert N  by M  matrix to a row of length N M
 66.20 find rows in a matrix based on values in different columns
 66.21 Select entries in one column based on a condition in another column
 66.22 Locate rows in a matrix with column being a string
 66.23 Remove set of rows and columns from a matrix at once
 66.24 Convert list of separated numerical numbers to strings
 66.25 Obtain elements that are common to two vectors
 66.26 Sort each column (on its own) in a matrix
 66.27 Sort each row (on its own) in a matrix
 66.28 Sort a matrix row-wise using first column as key
 66.29 Sort a matrix row-wise using non-first column as key
 66.30 Replace the first nonzero element in each row in a matrix by some value
 66.31 Perform outer product and outer sum between two vector
 66.32 Find the rank and the bases of the Null space for a matrix A
 66.33 Find the singular value decomposition (SVD) of a matrix
 66.34 Solve Ax  = b
 66.35 Find all nonzero elements in a matrix
 66.36 evaluate f (x)  on a vector of values
 66.37 generates equally spaced N points between x1   and x2
 66.38 evaluate and plot a f(x, y)  on 2D grid of coordinates
 66.39 Find determinant of matrix
 66.40 Generate sparse matrix with n  by n  matrix repeated on its diagonal
 66.41 Generate sparse matrix for the tridiagonal representation of the classic second difference operator on 1D grid
 66.42 Generate sparse matrix for the Laplacian differential operator ∇2u  for 2D grid
 66.43 Generate sparse matrix for the Laplacian differential operator ∇2u  for 3D grid
 66.44 Generate the adjugate matrix for square matrix
 66.45 Multiply each column by values taken from a row
 66.46 extract submatrix from a larger matrix by removing row/column
 66.47 delete one row from a matrix
 66.48 delete one column from a matrix
 66.49 generate random matrix so that each row adds to 1
 66.50 generate random matrix so that each column adds to 1
 66.51 sum all rows in a matrix
 66.52 sum all columns in a matrix
 66.53 find in which columns values that are not zero
 66.54 How to remove values from one vector that exist in another vector
 66.55 How to find mean of equal sized segments of a vector
 66.56 find first value in column larger than some value and cut matrix from there
 66.57 make copies of each value into matrix into a into a larger matrix
 66.58 How to apply a function to each value in a matrix?
 66.59 How to sum all numbers in a list (vector)?
67 Display spectrum of 2D image
68 Draw root locus for a discrete system
69 Plot the response of the inverted pendulum problem using state space
70 Convert a Laplace transfer function to an ordinary differential equation
71 Solve 2nd order ODE (Van Der Pol) and generate phase plot
72 Generate a plot using x  and y  data
73 Compare the effect on the step response of a standard second order system as ζ  changes
74 Find the linear convolution of 2 sequences where the origin is located at an arbitrary position in the sequence
75 Generate direction field plot of a first order differential equation
76 Obtain Fourier Series coefficients for a periodic function?
77 Plot the constant energy levels for a nonlinear pendulum
78 Obtain the statistical maximum likelihood estimates (MLE) of probability distributions
79 Make a histogram of data sampled from some probability distribution
80 Solve numerically the ODE  ′′′′
u   + u = f  using the point collocation method
81 Plot the dynamic response factor R
  d  of a system as a function of r = -ω-
    ωn  for different damping ratios
82 Visualize a 2D matrix
83 Find the particular and homogenous solution to undetermined system of equations
84 Numerically integrate f(x)  on the real line
85 Numerically integrate f(x,y)  in 2D
86 apply a filter on 1D numerical data (a vector)
87 apply an averaging Laplacian filter on 2D numerical data (a matrix)
88 How to compute ∑10   1
    k(k+1)
 k=1
89 How to find closed loop step response to a plant with a PID controller?
90 How to make table of x, y  values
91 How to find the moving average of a 1D sequence?
92 How to select N random values from a set of numbers?
93 How to build and connect a closed loop control systems and show the response?
 93.1 example 1, single input, single output closed loop
94 How to numerically solve a set of non-linear equations?
95 How to numerically solve Poisson PDE on 2D using Jacobi iteration method?
96 How to solve BVP second order ODE using finite elements with linear shape functions and using weak form formulation?
97 How to solve Poisson PDE in 2D using finite elements methods using rectangular element?
 97.1 Integrating the force vector
  97.1.1 f (x,y) = xy
 97.2 f(x,y) = 3 (cos(4πx ) + sin(3πy ))
  97.2.1 f (x,y) = 1
 97.3 f(x,y) = 1
 97.4 f(x,y) = xy
 97.5 f(x,y) = 3 (cos(4πx ) + sin(3πy ))
 97.6 Solving for reactangle grid
 97.7 f(x,y) = 1
 97.8 f(x,y) = xy
 97.9 f(x,y) = 3 (cos(4πx ) + sin(3πy ))
98 How to solve Poisson PDE in 2D using finite elements methods using triangle element?
 98.1 f (x,y ) = 1
 98.2 Program output
 98.3 f(x,y) = 1
99 How to draw an ellipse?
 99.1 draw an ellipse centered at origin and given its semi-major/minor (a,b)
 99.2 draw an ellipse centered at (x,y  ) and given its semi-major/minor (a,b)
 99.3 draw an ellipse centered at origin and given its semi-major/minor (a,b)  and rotated at angle 𝜃
100 How to solve wave equation using leapfrog method?

Introduction

This is a collection of useful 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.

If a bug is found in these examples, please let me know so I can correct them.

1 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

H (s) = -----25-----
        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]

PIC





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

PIC





Maple

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

PIC



2 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])}

PIC





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

PIC



This Maple solution will need to be updated so as to use the new DynamicSystems package.



Maple

  restart;
  with(inttrans):
  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"]);

PIC



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

PIC





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

PIC





Maple

  #This Maple solution will need to be
  #updated so as to use the new
  #DynamicSystems package.
  
  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"]);

PIC



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



5 Obtain the 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



6 Convert transfer function to state space representation

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

PIC

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

PIC

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



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



8 Generate Bode plot of a transfer function

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

        -----5s-----
H (s) = 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}}
   ]

PIC





Matlab

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

PIC



9 How to 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"]

PIC



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

PIC





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

PIC



10 Convert transfer function to state space

Problem: Given the continuous time S transfer function defined by

        ---1 +-s--
H (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]]

PIC



TransferFunctionModel[ss]

PIC





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



11 Obtain partial-fraction expansion

Problem: Given the continuous time S transfer function defined by

        s4-+-8s3-+-16s3-+-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



12 Obtain Laplace transform for a piecewise functions

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

PIC

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



13 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[%]]

PIC



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


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

PIC





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

PIC



15 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) =  s2 + 2ξω s + ω2
                 n     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
  ]

PIC





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

PIC



16 Convert differential equation to transfer functions and state space

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

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



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]

PIC



tf = Simplify[TransferFunctionModel[ss]]

PIC





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



17 View the steady state error of 2nd order LTI system as the undamped natural frequency changes

Problem: Given the transfer function

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

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

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

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

Use maximum time of 10  seconds and ξ = 0.707  and change omegan  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.

17.1 Mathematica

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

PIC



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

PIC



18 Convert a system from continuous time Laplace transfer function to discrete time Z transfer function

Problem: Convert

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

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



Mathematica

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

PIC





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.



19 Show the use of the inverse Z transform

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

19.1 example 1

Problem: Given

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

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



Mathematica

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

PIC





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

PIC



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

PIC





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

PIC



20 Find the Z transform of sequence x(n)

20.1 example 1

Find the Z transform for the unit step discrete function

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



Mathematica

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

Matlab

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



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



21 Sample a continuous time system

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

PIC

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

PIC



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

PIC





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

PIC





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

PIC



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

PIC

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

  First[Flatten[closedLoop[[1,1]]/closedLoop[[1,2]]]]

PIC

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
  



23 Linear convolution of 2 sequences

Problem: Given 2 sequences x1[n]  and x2 [m ]  , determine the linear convolution of the 2 sequences. Assume x1 = [1,2,3,4, 5]  and x2 = [3,4,5]  .



Mathematica

  Clear ["Global*"]
  x1 = {1, 2, 3, 4, 5};
  x2 = {4, 5, 6};
  ListConvolve[x2, x1, {1, -1}, 0]
  
  {4, 13, 28, 43, 58, 49, 30}

Matlab

  clear all; close all;
  x1=[1 2 3 4 5];
  x2=[4 5 6];
  conv(x2,x1)
  
  ans =
   4    13    28    43    58    49    30



24 Circular convolution of two sequences

Problem: Given 2 sequences x1 [n]  and x2[m ]  , determine the circular convolution of the 2 sequences.

24.1 example 1. The sequences are of equal length



Mathematica

  x1 = {2, 1, 2, 1};
  x2 = {1, 2, 3, 4};
  X1 = Chop[Fourier[x1, FourierParameters -> {1, -1} ]];
  X2 = Chop[Fourier[x2, FourierParameters -> {1, -1}]];
  Re[InverseFourier[X1*X2, FourierParameters -> {1, -1}]]
  
  {14., 16., 14., 16.}

Matlab

  clear all; close all;
  x1=[2 1 2 1];
  x2=[1 2 3 4];
  X1=fft(x1);
  X2=fft(x2);
  y=real(ifft(X1.*X2))
  
  y =
  
      14    16    14    16



24.2 example 2. The sequences are of unequal length



Mathematica

  Clear["Global*"]
  x1 = {1, 2, 3, 4, 5};
  x2 = {4, 5, 6};
  x2 = Append[x2, {0, 0}] // Flatten;
  X1 = Chop[Fourier[x1, FourierParameters -> {1, -1} ]];
  X2 = Chop[Fourier[x2, FourierParameters -> {1, -1}]];
  Re[InverseFourier[X1*X2, FourierParameters -> {1, -1}]]
  
  {53., 43., 28., 43., 58.}

Matlab

  x1=[1 2 3 4 5];
  x2=[4 5 6];
  N=max(length(x1),length(x2));
  X1=fft(x1,N);
  X2=fft(x2,N);
  y=real(ifft(X1.*X2))
  
  y =
  
      53    43    28    43    58



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

PIC

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

PIC

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



26 Plot the power spectral density of a signal

Problem: Given a signal

3sin(2πf1t ) + 4 cos(2πf2t) + 5cos (2πf3t)

for the following frequency values (in Hz)

f1 = 20, f2 = 30,f3 = 50

and plot the signal power spectral density so that 3 spikes will show at the above 3 frequencies. Use large enough sampling frequency to obtain a sharper spectrum



Mathematica

  Clear ["Global*"];
  f1 = 20;
  f2 = 30;
  f3 = 100;
  maxTime = 0.2;
  
  y[t_]:=2 Sin[2 Pi f1 t]+4 Cos[2 Pi f2 t]+
         5 Cos[2 Pi f3 t];
  
  Plot[y[t],{t,0,maxTime},
    Frame->True,
    FrameLabel->{{"y(t)",None},
        {"t (sec)","signal in time domain"}},
    RotateLabel->False,
    GridLines->Automatic,
    GridLinesStyle->Dashed,
    PlotRange->All]

PIC



  fs  = 7.0*Max[{f1,f2,f3}];
  yn  = Table[y[n],{n,0,maxTime,(1/fs)}];
  len = Length[yn];
  py  = Fourier[yn];
  nUniquePts = Ceiling[(len+1)/2];
  py  = py[[1;;nUniquePts]];
  py  = Abs[py];
  py  = 2*(py^2);
  py[[1]] = py[[1]]/2;
  f   = (Range[0,nUniquePts-1] fs)/len;
  
  ListPlot[Transpose[{f,py}],Joined->False,
    FrameLabel->{{"|H(f)|",None},
                 {"hz","Magnitude spectrum"}},
    ImageSize->400,
    Filling->Axis,
    FillingStyle->{Thick,Red},
    Frame->True,
    RotateLabel->False,
    GridLines->Automatic,
    GridLinesStyle->Dashed,
    PlotRange->All]
  

PIC



Here is the same plot as above, but using InterpolationOrder -> 0

PIC



And using InterpolationOrder -> 2

PIC





Matlab

  clear all; close all;
  maxTime = 0.2;
  t  = 0:0.001:maxTime;
  f1 =20; %Hz
  f2 =30; %Hz
  f3 =100; %Hz
  y = @(t) 2*sin(2*pi*f1.*t)+4*cos(2*pi*f2.*t)+...
           5*cos(2*pi*f3.*t);
  plot(t,y(t));
  set(gcf,'Position',[10,10,320,320]);
  grid on
  title('signal in time domain');
  xlabel('time(msec)');
  ylabel('Amplitude');

PIC



  fs  = 7.0*max([f1 f2 f3]);
  delT=1/fs;
  yn  = y(linspace(0,maxTime,maxTime/delT));
  len = length(yn);
  py  = fft(yn);
  nUniquePts = ceil((len+1)/2);
  py  = py(1:nUniquePts);
  py  = abs(py);
  py  = 2*(py.^2);
  py(1) = py(1)/2;
  f   = 0:nUniquePts-1;
  f   = f*fs/len;
  plot(f,py);
  title('power spectral density');
  xlabel('freq (Hz)');
  ylabel('Energy content');
  grid on
  set(gcf,'Position',[10,10,320,320]);

PIC



27 Solve the continuous-time algebraic Riccati equation

Problem: Solve for X  in the Riccati equation

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

given

pict


Mathematica

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

PIC





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



28 Solve the discrete-time algebraic Riccati equation

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

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

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

The Riccati equation is given by

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

Let R = [3]  .



Mathematica

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

PIC



ss = StateSpaceModel[dsys]

PIC



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

PIC





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



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

PIC



Plot the impulse response of the continuous system

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

PIC



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

PIC



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

PIC





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

PIC



30 Find the system type given an open loop transfer function

Problem: Find the system type for the following transfer functions

1.
s+21-
s −s
2.
-s+1-
s3−s2
3.
s+s51

To find the system type, the transfer function is put in the form -k∑i-(s−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



31 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

PIC



  N[Eigenvalues[a]]//MatrixForm

PIC



  N[Eigenvectors[a]]//MatrixForm

PIC





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



32 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



33 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

PIC

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



34 Solve Ax  = b  and display the solution

Problem: Solve for x given that Ax =  b  where

     (       )
        1  2
A =     1  3

and

    (    )
b =    5
       8

These 2 equations represent 2 straight lines in a 2D space. An exact solution exist if the 2 lines intersect at a point. The 2 lines are x + 2y = 5  and x + 3y = 8  .



Mathematica

  Remove["Global*"];
  ContourPlot[{x+2y==5,x+3y==8},
    {x,-3,2},
    {y,1.5,4},
    ContourStyle->{{Red,Thickness[0.01]},
                   {Blue,Thickness[0.01]}},
    GridLines->Automatic,
    GridLinesStyle->Dashed,
    ImageSize->300]

PIC



  mat = {{1,2},{1,3}};
  b   = {5,8};
  x   = LinearSolve[mat,b]

Out[58]= {-1,3}





Matlab

  A = [1 2;1 3];
  b = [5;8];
  x = A\b

  x =
      -1
       3



35 Determine if a set of linear equations Ax = b  has a solution and what type of solution

Problem: Given a general non homogeneous set of linear equations Ax =  b  how to test if it has no solution (inconsistent), or one unique solution, or an infinity number of solutions?

The following algorithm summarizes all the cases

  
  Let [A|b] be the augmented matrix, where b is appended to A.
  
  Assume A is an M by N matrix. i.e. M equations and N unknowns.
  
  IF rank A < rank [A|b]  THEN  -- system is inconsistent
     -- NO exact solution exist, but can use least square approximation
  
     x= A/b  -- Matlab.
     x= PseudoInverse[A].b  -- Mathematica uses SVD to computer pseduoInverse
  
  ELSE -- we must have rank A == rank[A|b]  -- system is consistent
       IF rank(A) == N  -- we have one solution.
          IF M==N -- one unique solution, can use crammer rule.
             x=A/b  -- Matlab. Here we get exact solution from \ operator
             x=LinearSolve[A,b] -- Mathematica
          ELSE  -- infinite solutions, pick one.
             x=A/b  -- Matlab. Here we get exact solution from \ operator
             x=LinearSolve[A,b] -- Mathematica
          END
       ELSE
          -- rank A  < N, i.e. rank deficient, infinite solutions. will pick one
          x=A/b
          x=LinearSolve[A,b]
       END
  END

Let the system of equations be

pict

So

     (         )
A  =     1  1
        − 2 1

and

    (      )
        1
b =    − 1



Mathematica

  Remove["Global*"];
  ContourPlot[{y==x-1,y==2 x+1} ,
              {x,-3,2},{y,-6,3},
       ContourStyle->{{Red,Thickness[0.01]},
                      {Blue,Thickness[0.01]}},
       GridLinesStyle->Dashed,
       ImageSize->300]

PIC



  a = {{-2,1},{-1,1}};
  b = {1,-1};
  {nRow, nCol} = Dimensions[a];
  aRank=MatrixRank[a]

  Out[]=         2



  abRank=MatrixRank[Insert[a,b,-1]]

  Out[]=         2



The above algorithm can now be run as follows

  If[aRank<abRank,
    Print["System is no consistent, no exact solution"];
    x=PseudoInverse[a].b,
    If[aRank==nCol,
     Print["System is consistent. Exact solution."];
     x=LinearSolve[a,b]
     ,
     Print["consistent, rank deficient,infinite solutions"];
     x=LinearSolve[a,b]
    ]
  ];
  
  Print["Solution is x=",x];

The output of the above is

  System is consistent. Exact solution.
  Solution is x={-2,-3}





Matlab

  A=[-2 1;-1 1];
  b=[1; -1];
  
  [nRow,nCol]=size(A);
  aRank=rank(A);
  abRank=rank([A b]);
  
  fprintf('A rank=%d, [A b] rank=%d\n',aRank,abRank);
  fprintf('Number of unknowns=%d\n',nCol);
  
  x=A\b;
  
  if aRank<abRank
   fprintf('System not consistent. no exact solution\n');
  else
   if aRank==nCol
    fprintf('System consistent. Exact solution.\n');
   else
    fprintf('consistent,rank deficient,infinite solutions\n');
   end
  end
  
  fprintf('solution is\n');
  x

  
  A rank=2, [A b] rank=2
  Number of unknowns=2
  
  System is consistent. Exact solution.
  solution is
  x =
      -2
      -3



36 Given a set of linear equations automatically generate the matrix A  and vector b  and solve Ax  = b

Problem: Given

4x + 4y + 2z   =   12
5x + 6y + 3z   =   7
9x + 7y + 10z  =   9

Automatically convert it to the form Ax  = b  and solve



Mathematica

  Remove["Global*"]
  eqs = {4x+3y+2z==12,
         5x+6y+3z  ==7,
         9x+7y+10z ==9};
  
  {b,a} = CoefficientArrays[eqs,{x,y,z}];
  
  Normal[a]//MatrixForm
  Normal[b]//MatrixForm

PIC



  LinearSolve[a,-b]//N

  Out[27]= {6.71429,-2.85714,-3.14286}



37 Convert a matrix to row echelon form and to reduced row echelon form

Problem: Given a matrix A, convert it to REF and RREF. Below shows how to convert the matrix A to RREF. To convert to REF (TODO).



Mathematica

  Remove["Global*"]
  (mat={{1, 1, -2,  1},
        {3, 2,  4, -4},
        {4, 3,  3, -4}})//MatrixForm

PIC



  MatrixForm[RowReduce[mat]]

PIC





Matlab

  clear all;
  A=[1 1 -2 1
     3 2 4  -4
     4 3 3  -4];
  rref(A)

  ans =
       1     0     0     2
       0     1     0    -3
       0     0     1    -1



38 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 = s2 + 4s + 25

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



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



39 Check continuous system stability in the Lyapunov sense

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

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

      − 1  − 2  − 3

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



Mathematica

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

PIC



  p = -LyapunovSolve[mat,IdentityMatrix[Length[mat]]];
  MatrixForm[N[p]]

PIC



  N[Eigenvalues[p]]

  {5.52534,1.78449,0.190164}



  r = Select[Re[N[Eigenvalues[p]]],# <= 0 &]

  Out[42]= {}



  If[Length[r]==0,
    Print["System is asymptotically stable"],
    Print["System is not asymptotically stable"]
  ]

   System is asymptotically stable





Matlab

  clear all;
  
  A=[0   1   0
     0   0   1
     -1 -2  -3];
  
  p=lyap(A,eye(length(A)))

  p =
  
      5.0000   -0.5000   -1.5000
     -0.5000    1.5000   -0.5000
     -1.5000   -0.5000    1.0000



  e=eig(p)

  e =
  
      0.1902
      1.7845
      5.525



  indx=find(real(e)<=0)

  indx =
  
     Empty matrix: 0-by-1



  if isempty(indx)
    fprintf('System is asymptotically stable\n');
  else
    fprintf('System is not asymptotically stable\n');
  end

   System is asymptotically stable



40 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

PIC

System block diagram.



Mathematica

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

PIC



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

PIC



  sys=SystemsModelSeriesConnect[d,plantd]

PIC



  loopBack=SystemsModelFeedbackConnect[sys]

PIC



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]

PIC



  ss=StateSpaceModel[loopBack]

PIC



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

PIC



  [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
  



41 Check if a matrix is Hermite

Problem: Given a matrix A, check to see if it is Hermite.

A Matrix is Hermite if it is the same as the conjugate of its transpose. One way to check is to take the difference and check that all values in the resulting difference matrix are zero.

To account for numerical values, the check is done using machine epsilon.



Mathematica

  mat = {{-1,1-3*I,0},
         {1+3*I,0,-6*I},
         {0,6*I,1}};
  MatrixForm[mat]

PIC



  mat2=Conjugate[Transpose[mat]];
  (diff=mat-mat2)//MatrixForm

PIC



Now check that each element in the difference matrix is less than MachineEpsilon

  Map[Abs[#]<$MachineEpsilon&,diff,{2}]

  Out[32]= {{True,True,True},
            {True,True,True},
            {True,True,True}}



  Count[%,False]

  Out[33]= 0





Matlab

  clear all;
  i=sqrt(-1);
  mat=[-1     1-3*i    0;
      1+3*i    0      -6*i;
        0      6*i     1];
  
  mat2=conj(transpose(mat));
  diff = mat-mat2

  diff =
       0     0     0
       0     0     0
       0     0     0



  r=abs(diff)<eps('double')

  r =
       1     1     1
       1     1     1
       1     1     1



  all(r(:))

       1



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

PIC



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

PIC





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

PIC



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

PIC



  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]

PIC





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

PIC



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

PIC



  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]

PIC





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

PIC



45 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) =  s(s +-4)(s-+-6)(s2 +-1.4s-+-1)

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



Mathematica

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

PIC





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

PIC



46 Find the cross correlation between two sequences

Problem: Given

pict

Notice that the output sequence generated by Mathematica and Matlab are reversed with respect to each others.

Also, MATLAB uses the length 2N  − 1  as the length of cross correlation sequence, which in this example is 23 because N  is taken as the length of the larger of the 2 sequences if they are not of equal length which is the case in this example.

In Mathematica, the length of the cross correlation sequence was 22, which is 2N  .



Mathematica

  Clear["Global*"];
  
  a={0,0,2,-1,3,7,1,2,-3,0,0};
  b={0,0,1,-1,2,-2,4,1,-2,5,0,0};
  c=Reverse[ListCorrelate[a,b,{-1,1},0]]
  
  Out[31]= {0,
            0,
            0,
            0,
            10,
            -9,
            19,
            36,
            -14,
            33,
            0,
            7,
            13,
            -18,
            16,
            -7,
            5,
            -3,
            0,
            0,
            0,
            0}
  

Matlab In MATLAB use the xcross in the signal processing toolbox

  clear all; close all;
  A=[0,0,2,-1,3,7,1,2,-3,0,0];
  B=[0,0,1,-1,2,-2,4,1,-2,5,0,0];
  C=xcorr(A,B);
  format long
  C'
  
  ans =
  
     0.000000000000003
     0.000000000000002
    -0.000000000000002
                     0
     9.999999999999998
    -9.000000000000002
    19.000000000000000
    36.000000000000000
   -14.000000000000000
    33.000000000000000
    -0.000000000000002
     6.999999999999998
    13.000000000000000
   -18.000000000000000
    16.000000000000004
    -7.000000000000000
     4.999999999999999
    -2.999999999999998
    -0.000000000000000
     0.000000000000001
     0.000000000000002
    -0.000000000000004
                     0



47 Find the different norms of a vector

Problem: Given the vector say

v =  1,2,4

Find its norm for p = 1, 2,∞



Mathematica

  Remove["Global*"];
  
  v = {{1}, {2}, {4}};
  p = {1, 2, Infinity};
  Map[Norm[v,#]&,p]
  
  Out[8]= {7,Sqrt[21],4}

Matlab

  clear all; close all;
  v=[1 2 4];
  p=[1,2,inf];
  arrayfun(@(i) norm(v,p(i)),1:length(p))
  
  >>>
  ans =
  
      7.0000    4.5826    4.0000



48 Find orthonormal vectors that span the range of matrix A

Problem: Given the matrix A  whose columns represents some vectors, find the set of orthonormal vectors that span the same space as A  and verify the result. Let

     ⌊           ⌋
      0   1  1  2
A  = ⌈1   2  3  4⌉

      2   0  2  0

Notice that A  has rank 2, so we should get no more than 2 vectors in the orthonormal set.

With MATLAB use the orth(A) function, With Mathematica, use {u,s,v}=SingularValueDecomposition[A] , and since the rank is 2, then the first 2 columns of matrix u will give the answer needed (any 2 columns of u will also give a set of orthonormal vectors).



Mathematica

  Remove["Global*"];
  mat = {{0, 1, 1, 2},
         {1, 2, 3, 4},
         {2, 0, 2, 0}};
  r = MatrixRank[mat]
  
  Out[15]= 2
  
  {u,s,v}=SingularValueDecomposition[mat];
  
  orth=N[u[[All,{1,r}]]]
  
  Out[75]= {{0.378151, -0.308379},
            {0.887675, -0.146825},
            {0.262747, 0.939864}}
  
  Chop[Transpose[orth].orth]
  
  Out[74]= {{1.,0},
            {0,1.}}

Matlab

  clear all;
  A=[0 1 1 2
     1 2 3 4
     2 0 2 0];
  R=orth(A)
  >>
  R =
     -0.3782    0.3084
     -0.8877    0.1468
     -0.2627   -0.9399
  
  R'*R
  >>
      1.0000    0.0000
      0.0000    1.0000



49 Plot the surface described by f (y,x )

Problem: Plot the function

f(x,y) = x2 − 2xy +  3y + 2



Mathematica

  Remove["Global*"];
  f[x_,y_]:=x^2-2x y+3y+2;
  Plot3D[f[x,y],{x,-4,4},{y,-3,3},
      PlotLabel->"Plot of x^2-2x y+3y+2",
      ImageSize->300]

PIC



Matlab

  clear all; close all;
  x = -4:.3:4;
  y = -3:.3:3;
  [X,Y] = meshgrid(x,y);
  Z = X.^2 - 2*(X.*Y) + 3*Y + 2;
  surf(X,Y,Z)
  title('plot of z=x^2-2(x y)+3y+2');
  set(gcf,'Position',[10,10,420,420]);

PIC



50 Find the point of intersection of 3 surfaces

Problem: Given 3 surfaces find the point where they intersect by solving the linear system and display the solution.

Let the 3 surfaces be

pict

In Mathematica, plot the surfaces using the Plot3D command, and then use LinearSolve to solve the system of equations in order to find point of intersection, then use ScatterPlot to plot the point solution. In Matlab, use surf command to do the plotting, and solve the system of equations using Matlab



Mathematica

  Remove["Global*"];
  eq1 = z-2x-y ==0;
  eq2 = z-2y ==0;
  eq3 = z- 2 ==0;
  p1 = ContourPlot3D[Evaluate@eq1,
           {x,-5,5},{y,-5,5},{z,-10,10},
           AxesLabel->{"x","y","z"},
           PlotLabel->"z=2x+y",
           ImageSize->250];
  p2 = ContourPlot3D[Evaluate@eq2,
          {x,-5,5},{y,-5,5},{z,-10,10},
          AxesLabel->{"x","y","z"},
          PlotLabel->"z=2x",
          ImageSize->250];
  p3 = ContourPlot3D[Evaluate@eq3,
          {x,-5,5},{y,-5,5},{z,-10,10},
          AxesLabel->{"x","y","z"},
          PlotLabel->"z=2",
          ImageSize->250];
  Grid[{{p1,p2,p3}},Frame->All]

PIC



  {b,a} = CoefficientArrays[{eq1,eq2,eq3},
           {x,y,z}];
  sol   = LinearSolve[a,-b]//N
  
  Out[301]= {0.5,1.,2.}
  
  Show[p1,p2,p3,
    Graphics3D[{PointSize[0.2],
         Red,Point[sol]}],
    PlotLabel->"point where surfaces meet)"]

PIC





Matlab

  clear all; close all;
  x = -4:.5:4;
  y = -3:.5:3;
  
  [X,Y] = meshgrid(x,y);
  
  subplot(1,3,1);
  
  Z = 2*X + Y;  surf(X,Y,Z);
  title('z=2x+y');
  
  subplot(1,3,2);
  Z = 2*Y;  surf(X,Y,Z);
  title('z=2y');
  
  subplot(1,3,3);
  Z = zeros(length(y),length(x));
  Z(1:end,1:end)=2;
  surf(X,Y,Z);
  title('z=2');
  
  A=[-2 -1 1;0 -2 1;0 0 1];
  b=[0;0;2];
  sol=A\b
  
  figure;
  surf(X,Y,2*X + Y); hold on;
  surf(X,Y,2*Y);
  surf(X,Y,Z);
  
  %now add the point of intersection
  %from the solution
  
  plot3(sol(1),sol(2),sol(3),...
        'k.','MarkerSize',30)
  
  title('solution to system of linear equations.');

  sol =
  
      0.5000
      1.0000
      2.0000
  

PIC

PIC



51 Obtain the LU decomposition of a matrix

Problem: Given a matrix A  , find matrix L  and U  such that LU  =  A  . The matrix L  will have 1  in all the diagonal elements and zeros in the upper triangle. The matrix U  will have 0  in the lower triangle as shown in this diagram.

PIC



Mathematica

  Remove["Global*"]
  mat = {{1, -3, 5},
         {2, -4, 7},
         {-1, -2, 1}};
  MatrixForm[mat]

PIC

  {lu,pivotVector,conditionNumber}=
       LUDecomposition[mat]
  
  Out[32] = {{{1,-3,5},
              {2,2,-3},
              {-1,-(5/2),-(3/2)}},
              {1,2,3},1}
  
  lower=lu*SparseArray[{i_,j_}/;j<i->1,{3,3}]+
           IdentityMatrix[3]
  Out[]= {{1,0,0},{2,1,0},{-1,-(5/2),1}}
  
  MatrixForm[lower]

PIC

  upper=lu SparseArray[{i_,j_}/;j>=i->1,{3,3}]
  Out[] = {{1,-3,5},{0,2,-3},{0,0,-(3/2)}}
  MatrixForm[upper]

PIC

  lower.upper
  Out[8]= {{1,-3,5},{2,-4,7},{-1,-2,1}}
  MatrixForm[%]

PIC

Matlab

  clear all;
  
  A=[1 -3 5;
     2 -4 7;
     -1 -2 1];
  
  [L,U,p]=lu(A)
  >>>
  L =
      1.0000         0         0
     -0.5000    1.0000         0
      0.5000    0.2500    1.0000
  
  U =
      2.0000   -4.0000    7.0000
           0   -4.0000    4.5000
           0         0    0.3750
  p =
       0     1     0
       0     0     1
       1     0     0
  >>
  L*U
  
  ans =
       2    -4     7
      -1    -2     1
       1    -3     5



52 Use FFT to display the power spectrum of the content of an audio wav file

Problem: download a wav file and display the frequency spectrum of the audio signal using FFT. The Matlab example was based on Matheworks tech note 1702.



Mathematica

  Remove["Global*"];
  fname = "stereo.wav";
  SetDirectory[ToFileName[Extract
   ["FileName"/.NotebookInformation[
   EvaluationNotebook[]],{1},FrontEndFileName]]];
  
  ele = Import[fname,"Elements"]

  {AudioChannels,AudioEncoding,
   Data,SampledSoundList,SampleRate,Sound}



  (*listen to it *)
  sound = Import[fname,"Sound"];
  EmitSound[sound]
  
  fs = Import[fname,"SampleRate"]

  Out[113]= 22050



  (*now load the samples and do power spectrum*)
  data = Import[fname,"Data"];
  {nChannel,nSamples}=Dimensions[data]

  Out[115]= {2,29016}



  py = Fourier[data[[1,All]],
               FourierParameters->{1,-1}];
  nUniquePts = Ceiling[(nSamples+1)/2]

  Out[117]= 14509



  py = py[[1;;nUniquePts]];
  py = Abs[py];
  py = py/nSamples;
  py = py^2;
  
  If[ OddQ[nSamples],
    py[[2;;-1]]=2*py[[2;;-1]],
    py[[2;;-2]]=2*py[[2;;-2]]
  ];
  
  f = (Range[0,nUniquePts-1] fs)/nSamples;
  ListPlot[Transpose[{f,py}],
    Joined->True,
    FrameLabel->{{"|H(f)|",None},
                 {"hz","Magnitude spectrum"}},
    ImageSize->400,
    Frame->True,
    RotateLabel->False,
    GridLines->Automatic,
    GridLinesStyle->Dashed,
    PlotRange->{{0,12000},All}]

PIC





Matlab

  [data,Fs,nbits] = wavread('stereol.wav');
  [nSamples,nChannels] = size(data);
  waveFileLength = nSamples/Fs;
  N = 2^(nextpow2(length(data(:,1))));
  Y = fft(data(:,1),N);
  NumUniquePts = ceil((N+1)/2);
  Y = Y(1:NumUniquePts);
  P = abs(Y);
  P = P/length(data(:,1));
  P = P.^2;
  
  if rem(N, 2)
    P(2:end) = P(2:end)*2;
  else
    P(2:end -1) = P(2:end -1)*2;
  end
  
  f = (0:NumUniquePts-1)*Fs/N;
  plot(f,P);
  title(sprintf('Power Spectrum of a wave file'));
  xlabel('Frequency Hz');
  ylabel('|H(f)|');

PIC



53 Solve homogeneous 1st order linear differential equation with constant coefficients and initial conditions

Problem: Solve

 ′
y (t) = 3y (t)

with initial condition y (0) = 1  and plot the solution for t = 0⋅⋅⋅2  seconds.



Mathematica

  Remove["Global*"];
  sol = First@DSolve[{y'[t]==3y[t],
           y[0]==1},y[t],t];
  y   = y[t]/.sol
  
  Out[26]= E^(3 t)
  
  Plot[y,{t,0,2},
    FrameLabel->{{"y(t)",None},
     {"t","Solution of y'=3y, y(0)=1"}},
    Frame->True,
    GridLines->Automatic,
    GridLinesStyle->Automatic,
    RotateLabel->False,
    ImageSize->300,
    AspectRatio->1]

PIC

Matlab

  function e53
  t         = 0:0.001:2;   % time
  initial_y = 1;
  
  [t,y] = ode45( @rhs, t, initial_y);
  
  plot(t,y)
  title('Solution of y''=3y , y(0)=1, using ode45');
  xlabel('time');
  ylabel('y(t)');
  grid on
  set(gcf,'Position',[10,10,420,420]);
  
      function dydt=rhs(t,y)
          dydt = 3*y;
      end
  end

PIC



54 Solve homogeneous 2nd order linear differential equation with constant coefficients and initial conditions

Problem: Solve

 ′′         ′
y (t) − 1.5y (t) + 5y(t) = 0

with initial conditions

           ′
y(0) = 1,y  (0 ) = 0

To use Matlab ode45, the second order ODE is first converted to state space formulation as follows

Given y′′(t) − 1.5y′(t) + 5y (t) = 0  let

pict

hence

x′=  x2
 1

and

pict

Hence we can now write

[  ]   [        ] [  ]
 x′1       0   1    x1
 x′  =   − 5  1.5   x
  2                  2

Now Matlab ODE45 can be used.



Mathematica

  Remove["Global*"];
  eq  = y''[t]-1.5y'[t]+5y[t]==0;
  ic  = {y'[0]==0,y[0]== 1};
  sol = First@DSolve[{eq,ic},y[t],t];
  y   = y[t]/.sol
  
  Out[51]= 1. E^(0.75 t) Cos[2.10654 t]-
           0.356034 E^(0.75 t) Sin[2.10654 t]
  
  Plot[y,{t,0,10},
       FrameLabel->{{"y(t)",None},
           {"t","Solution"}},
       Frame->True,
       GridLines->Automatic,
       GridLinesStyle->Automatic,
       RotateLabel->False,
       ImageSize->300,
       AspectRatio->1,
       PlotRange->All,
       PlotStyle->{Thick,Red}]

PIC

Matlab

  function e54
  
  t0 = 0; %initial time
  tf = 10; %final time
  ic =[1 0]'; %initial conditions [y(0)  y'(0)]
  
  [t,y] = ode45(@rhs, [t0 tf], ic);
  
  plot(t,y(:,1),'r')
  title('Solution using ode45');
  xlabel('time');
  ylabel('y(t)');
  grid on
  set(gcf,'Position',[10,10,320,320]);
  
      function dydt=rhs(t,y)
          dydt=[y(2) ;
               -5*y(1)+1.5*y(2)];
      end
  end

PIC



55 Solve non-homogeneous 2nd order linear differential equation with constant coefficients and initial conditions

Problem: Solve

  ′′         ′
y  (t) − 1.5y (t) + 5y (t) = 4cos (t)

with initial conditions

           ′
y(0) = 1,y  (0 ) = 0

To use Matlab ode45, the second order ODE is converted to state space as follows

Given y′′(t) − 1.5y′(t) + 5y (t) = 4 cos(t)  , let

pict

hence

x′=  x2
 1

and

pict

Hence we can now write

[   ]   [        ] [  ]   [  ]
  x′      0    1    x1     0
   1′  =                 +     cos (t)
  x2     − 5  1.5   x2     4



Mathematica

  Remove["Global*"];
  eq  = y''[t]-1.5y'[t]+5y[t]==4  Cos[t];
  ic  = {y'[0]==0,y[0]== 1};
  sol = First@DSolve[{eq,ic},y[t],t];
  
  y=Chop[y[t]/.sol]
  
  Out[65]= 0.123288 E^(0.75 t) Cos[2.10654 t]+
           0.876712 Cos[t] Cos[2.10654 t]^2-
           0.328767 Cos[2.10654 t]^2 Sin[t]+
           0.112175 E^(0.75 t) Sin[2.10654 t]+
           0.876712 Cos[t] Sin[2.10654 t]^2-
           0.328767 Sin[t] Sin[2.10654 t]^2
  
  
  Plot[y,{t,0,10},FrameLabel->{{"y(t)",None},
        {"t","Solution"}},
       Frame->True,
       GridLines->Automatic,
       GridLinesStyle->Automatic,
       RotateLabel->False,
       ImageSize->300,
       AspectRatio->1,
       PlotRange->All,
       PlotStyle->{Thick,Red}]

PIC

Matlab

  function e55
  t0 = 0; %initial time
  tf = 10; %final time
  ic =[1 0]'; %initial conditions [y(0)  y'(0)]
  
  [t,y] = ode45(@rhs, [t0 tf], ic);
  
  plot(t,y(:,1),'r')
  title('Solution using ode45');
  xlabel('time');
  ylabel('y(t)');
  grid on
  set(gcf,'Position',[10,10,320,320]);
  
      function dydt=rhs(t,y)
        dydt=[y(2) ; -5*y(1)+1.5*y(2)+4*cos(t)];
      end
  end
  

PIC



56 Solve homogeneous 2nd order linear differential equation with constant coefficients with boundary values (BVP)

Problem: Solve

y ′′(t) + t y(t) = 0

with the boundary conditions

y(0) = 3,y (20) = − 1

For solving with Matlab, the ODE is first converted to state space as follows

Given  ′′
y  (t) + t y (t) = 0  , let

pict

Therefore

 ′
x1 = x2

And

pict

This results in

[  ]   [      ] [  ]
 x′1      0   1   x1
 x′  =   − t 0   x
  2                2

Now bvp4c() can be used.



Mathematica

  Clear[y,t];
  eq  = y''[t]+t y[t]==0;
  ic  = {y[0]==3,y[20]==-1};
  sol = First@DSolve[{eq,ic},y[t],t];
  
  y   = y[t]/.sol
  
  Out[69]=(-Sqrt[3] AiryAi[(-1)^(1/3) t]+
   AiryBi[(-1)^(1/3) t]-3 3^(2/3)
   AiryAi[(-1)^(1/3) t]
   AiryBi[20 (-1)^(1/3)] Gamma[2/3]+3 3^(2/3)
   AiryAi[20 (-1)^(1/3)] AiryBi[(-1)^(1/3) t]
   Gamma[2/3])/(Sqrt[3] AiryAi[20 (-1)^(1/3)]-
      AiryBi[20 (-1)^(1/3)])
  
  Plot[Re[y],{t,0,20},
       FrameLabel->{{"y(t)",None},
        {"t","Solution"}},
       Frame->True,
       GridLines->Automatic,
       GridLinesStyle->Automatic,
       RotateLabel->False,
       ImageSize->300,
       AspectRatio->1,
       PlotRange->All,
       PlotStyle->{Thick,Red},
       Exclusions->None]
  

PIC

Matlab

  function e56
  
  t0 = 0; %initial time
  tf = 20; %final time
  
  solinit = bvpinit(linspace(t0,tf,5),[3 -1]);
  
  sol = bvp4c(@odefun,@bcfun,solinit);
  
  plot(sol.x(1,:),sol.y(1,:),'r')
  title('solution');
  xlabel('time');
  ylabel('y(t)');
  grid;
  set(gcf,'Position',[10,10,320,320]);
  
  function dydt=odefun(t,y)
  dydt=[y(2); -t.*y(1)];
  
  function res=bcfun(yleft,yright)
  res=[yleft(1)-3
      yright(1)+1];
  

PIC



57 Solve the 1-D heat partial differential equation (PDE)

The PDE is

∂T  (x, t)    ∂2T (x, t)
---------=  k------2---
   ∂t           ∂x

Problem: given a bar of length L  and initial conditions T (x, 0) = sin (πx)  and boundary conditions T (0,t) = 0,T (L,t) = 0  , solve the above PDE and plot the solution on 3D.

Use bar length of 4  and k = 0.5  and show the solution for 1  second.



Mathematica

  Clear[y,x,t];
  barLength = 4;
  timeDuration = 1;
  eq = D[y[x, t], t] == k*D[y[x, t], x, x];
  eq = eq /. k -> 0.5;
  boundaryConditions = {y[0,t]==0,
                       y[barLength,t]==0};
  initialConditions = y[x,0]==Sin[Pi x];
  sol=First@NDSolve[{eq,
                     boundaryConditions,
                     initialConditions},
                     y[x,t],
                     {x,0,barLength},
                     {t,0,timeDuration}];
  y = y[x,t]/.sol
  Out[138]= InterpolatingFunction[
               {{0.,4.},{0.,1.}},<>][x,t]
  
  Plot3D[y,{x,0,barLength},{t,0,timeDuration},
         PlotPoints->30,
         PlotRange->All,
         AxesLabel->{"x","time","T[x,t]"},
         ImageSize->300]

PIC

  ContourPlot[y,{x,0,barLength},
                {t,0,timeDuration},
              PlotPoints->15,
              Contours->15,
              ContourLines->False,
              ColorFunction->Hue,
              PlotRange->All,
              ImageSize->300]

PIC

Matlab

  function e57
  
  m = 0;
  x = linspace(0,4,30);
  t = linspace(0,1,30);
  
  sol = pdepe(m,...
              @pde,...
              @pdeInitialConditions,...
              @pdeBoundaryConditions,x,t);
  
  u = sol(:,:,1);
  
  surf(x,t,u)
  title('Numerical PDE solution')
  xlabel('x')
  ylabel('Time t')
  set(gcf,'Position',[10,10,320,320]);
  
  % ----------------
  function [c,f,s] = pde(x,t,u,DuDx)
  k=0.5;
  c = 1/k;
  f = DuDx;
  s = 0;
  % ---------------
  function T0 = pdeInitialConditions(x)
  T0 = sin(pi*x);
  % ---------------
  function [pl,ql,pr,qr] = ...
     pdeBoundaryConditions(xl,ul,xr,ur,t)
  pl = ul;
  ql = 0;
  pr = 0;
  qr = 1;

PIC



58 Show the effect of different boundary/initial conditions on the solution of the 1-D heat partial differential equation (PDE)

The PDE is

∂T  (x, t)    ∂2T (x, t)
---------=  k------2---
   ∂t           ∂x

Problem: given a bar of length L  , solve the above 1-D heat PDE for 4 different boundary/initial condition to show that the solution depends on these.

Mathematica

  Clear[y,x,t,k];
  barLength = 4*Pi;
  timeDuration = 10;
  eq = D[y[x, t], t] == k*D[y[x, t], x, x];
  eq = eq /. k -> 0.5
  solveHeat[eq_,bc_,ic_,y_,x_,t_,len_,timeDuration_]:=Module[{sol},
  
    sol=First@NDSolve[{eq,bc,ic},
                      y[x,t],
                      {x,0,len},
                      {t,0,timeDuration}];
  
    Plot3D[y[x,t]/.sol,{x,0,barLength},{t,0,timeDuration},
           PlotPoints->30,
           PlotRange->All,
           AxesLabel->{"x","time","y[x,t]"},
           ImageSize->200,
           PlotLabel->bc]
  ];
  
  bc = {{y[0,t] == 1, y[barLength,t] == 1},
        {y[0,t] == 0, y[barLength,t] == 0},
        {y[0,t] == 1, y[barLength,t] == Exp[-t barLength]},
        {y[0,t] == 0, y[barLength,t] == 0}
       };
  
  ic = {y[x,0] == Cos[x],
        y[x,0] == Sin[x],
        y[x,0] == 1,
        y[x,0] == Cos[x]Sin[x]
        };
  
  sol = MapThread[
          solveHeat[eq,#1,#2,y,x,t,barLength,timeDuration]&,
            {bc,ic}
        ];
  
  Grid[Partition[sol,2],Frame->All]

Each plot shows the boundary conditions used.

PIC

59 Generate uniform distributed random numbers

59.1 How to generate 5 uniform distributed random numbers from 0 to 1?



Mathematica

  SeedRandom[1];
  Table[RandomVariate[
     UniformDistribution[{0,1}]],{5}]
  
  Out[66]= {0.817389,
            0.11142,
            0.789526,
            0.187803,
            0.241361}

Matlab

  rand(5,1)
  
  ans =
  
     0.814723686393179
     0.905791937075619
     0.126986816293506
     0.913375856139019
     0.632359246225410



Fortran

  program t3
  implicit none
  real :: x(5)
  
   CALL RANDOM_SEED()
   CALL random_number(x)
   print *,x
  
  end program
  
  #compile and run
  
  $ gfortran -std=f95 -Wextra -Wall -pedantic -funroll-loops
    -ftree-vectorize -march=native  -Wsurprising -Wconversion
    t3.f90  /usr/lib/liblapack.a /usr/lib/libblas.a
  $ ./a.exe
    0.99755955      0.56682467      0.96591532      0.74792767      0.36739087
  

59.2 How to generate 5 random numbers from a to b?

Generate uniform numbers from a to b, say a=-2 and b=5



Mathematica

  SeedRandom[1];
  Table[RandomVariate[
      UniformDistribution[{-2,5}]],{5}]
  
  Out[70]= {3.72173,
           -1.22006,
            3.52668,
           -0.685378,
           -0.310473}

Matlab

  -2 + (5+2)*rand(5,1)
  
  ans =
    -1.317217165004133
    -0.050512467930661
     1.828170634434887
     4.702547848040084
     4.754219746394936



Fortran

  program t3_2
  implicit none
  integer ::i
  real, parameter :: a=-1.0,b=1.0
  real :: x(5)
  
   CALL RANDOM_SEED()
   DO i=1,2
      CALL random_number(x)
      x = a+(b-a)*x
      print *,x
   END DO
  
  end program
  
  #compile and run
  
  $ gfortran -std=f95 -Wextra -Wall -pedantic -funroll-loops
    -ftree-vectorize -march=native  -Wsurprising -Wconversion
    t3_2.f90  /usr/lib/liblapack.a /usr/lib/libblas.a
  $ ./a.exe
    0.99511909      0.13364935      0.93183064      0.49585533     -0.26521826
   -3.87262106E-02 -0.85249150     -0.98928964     -0.30583751     -0.31551242
  $

59.3 How to generate MATRIX of random numbers from a to b?

Let a = − 2  and b = 5  , matrix of size 5  by 5



Mathematica

  SeedRandom[1];
  Table[RandomVariate[
     UniformDistribution[{-2,5}]],{5},{5}]
  
  {{3.72173,-1.22006,3.52668,-0.685378,-0.310473},
   {-1.53983,1.79573,-0.381918,0.772043,2.90332},
   {-0.517218,3.2406,0.959955,-0.267537,4.8402},
   {3.77614,4.47693,2.04639,0.0500882,-0.543643},
   {2.06332,-1.09825,0.144992,2.98408,0.734073}}

Matlab

  -2 + (5+2)*rand(5,5)
  
  ans =
   3.7642    1.0712    1.4284   -0.0678    1.4885
   2.8638    0.6709    1.1191    2.7579    4.7182
   0.2197    3.3586    2.5242    2.5857    0.3827
   4.6516    3.5664    2.9656   -0.8617    2.0969
   -1.7589   -0.6919    3.2828   -1.1670   -0.4333



Fortran

  program t3_3
  implicit none
  integer ::i
  real, parameter :: a=-1.0,b=1.0
  real :: x(5,5)
  
   CALL RANDOM_SEED()
   DO i=1,2
      CALL random_number(x)
      x = a+(b-a)*x
      CALL write(x)
   END DO
  
   !--- internal functions below ------
       contains
         SUBROUTINE write(A)
         implicit none
         REAL, DIMENSION(:,:) :: A
         integer :: i,j
  
         WRITE(*,*)
         DO i = lbound(A,1), ubound(A,1)
            WRITE(*,*) (A(i,j), j = lbound(A,2), ubound(A,2))
         END DO
  
    END SUBROUTINE write
  
  end program
  
  #compile and run
  
  $ gfortran -std=f95 -Wextra -Wall -pedantic -funroll-loops
      -ftree-vectorize -march=native  -Wsurprising -Wconversion
      t3_3.f90  /usr/lib/liblapack.a /usr/lib/libblas.a
  
  $ ./a.exe
  
    0.99511909     -3.87262106E-02 -0.56409657      0.32386434      0.71138477
    0.13364935     -0.85249150     -0.73367929     -0.96778345     -0.19742620
    0.93183064     -0.98928964      0.80104899      0.30170965     -0.58625138
    0.49585533     -0.30583751     -0.22646809      0.29281759      0.93707883
   -0.26521826     -0.31551242     -0.10903549     -0.35402548      0.19679904
  
    0.34596145      0.21138644      0.22462976      0.31809497      0.45771694
   -8.62354040E-02  0.43809581      0.95732045      0.10801017     -0.19508958
   -0.33996975      0.79466915      0.99828446      0.95552015      0.85725522
   -0.79923415      0.31645823     -0.48640406      0.80384660     -0.70432973
    0.51090658     -0.69856644      0.10173070      0.31584930      0.34905851

60 Determine and plot the Fourier Transform for a continuous time function

Plot the Fourier Transform of 3 sin(t)



Mathematica

  Clear["Global*"];
  f = 3*Sin[t];
  y = FourierTransform[f, t, w, FourierParameters -> {1, -1}]
  
  Out[138]= -3 I Pi DiracDelta[-1+w]+3 I Pi DiracDelta[1+w]
  
  tab = Table[{k,y/.w->k},{k,-3,3}]
  tab = tab/.DiracDelta[0]->1
  tab[[All,2]] = Map[Sign[Im[#]]Abs[#]&,tab[[All,2]]]
  
  ListPlot[tab,PlotStyle->AbsolutePointSize[5],
           AxesLabel->{"frequency \[CapitalOmega] rad/sec",
           "|F(\[CapitalOmega]|"},
           PlotRange->{All,{-12,12}},
           Filling->Axis,
           FillingStyle->{{Thick,Red}}]

PIC

Matlab

  syms t;
  F=fourier(3*sin(t))
  
  >>
  F =
  
  -pi*(dirac(w-1)-dirac(w+1))*3*i
  

Need to do the plot.



61 Solve the 2-D Laplace PDE ∇2T  (x,y ) = 0  for a rectangular plate with Dirichlet boundary conditions

Problem: Given the following plate, with height h =  30  , and width w  = 10  , and with its edges held at fixed temperature as shown, find the steady state temperature distribution

PIC

System boundary conditions.

Mathematica NDSolve[] does not currently support Laplace PDE as it is not an initial value problem.

Jacobi iterative method is used below to solve it. 100 iterations are made and then the resulting solution plotted in 3D.



Mathematica

  n = 32;  (*grid size*)
  h = 1/(n-1);(*grid spacing*)
  u = Table[0.,{n},{n}]; (*init to zero*)
  u[[1,All]] = 100; (*B.C.*)
  
  Do[ (*iterate 100 times*)
    tmp=u;
    For[i=2,i<=n-1,i++,
      For[j=2,j<=n-1,j++,
        tmp[[i,j]]=
          (1/4)(u[[i-1,j]]+u[[i+1,j]]+
             u[[i,j-1]]+u[[i,j+1]])
      ]
    ];
  
    u=tmp,
    {100}
  ];
  
  ListPlot3D[u,PlotRange->All,
   ImagePadding->20,
   Mesh->{n,n},
   ImageSize->300,
   PlotLabel->"solution of Laplace PDE on 2D"
  ]

PIC

Matlab

  clear all; close all;
  n = 32;
  h = 1/(n-1); %grid spacing
  u = zeros(n,n);
  
  u(1,:) = 100; %B.C.
  % coordinates
  [X,Y]  = meshgrid(0:h:1,0:h:1);
  
  for k=1:100
      tmp = u;
  
      for i=2:n-1
          for j=2:n-1
            tmp(i,j)=...
             (1/4)*(u(i-1,j)+u(i+1,j)+...
               u(i,j-1)+u(i,j+1));
          end
      end
  
      u=tmp;
  end
  
  mesh(X,Y,u);
  title('solution of Laplace PDE on 2D');
  set(gcf,'Position',[10,10,320,320]);

PIC



62 Generate and plot one pulse signal of different width and amplitude

Problem: Generate one signal of some width and height.



Mathematica

  Remove["Global*"]
  f[x_,w_,h_]:=If[0<=x<w,
        h*If[x-IntegerPart[x]<=0.5,
        SquareWave[x],
        -SquareWave[x]
        ],
        0]
  w = {0.5,1,2};   (*pulse width *)
  h = {1,1.5,0.4}; (*pulse height *)
  plots = MapThread[Plot[
            f[x,#1,#2],{x,-2,2.5},
              PlotRange->{{-1,2.5},{-0.3,2.5}},
              Frame->True,
              Axes->False,
              GridLines->Automatic,
              GridLinesStyle->Dashed,
              PlotStyle->{Thick,Red},
              FrameLabel->{{"y(t)",None},
      {"t","pulse of "<>ToString[#1]<>
      " sec width and amplitude"<>
      ToString[#2]}},
              ImageSize -> 160,
              AspectRatio->1]&,{w,h}];
  
  Grid[{plots},Frame->All]

PIC

Matlab

  clear all; close all;
  f = [1 0.5 0.25];
  w = [0.5 1 2];
  h = [1 1.5 0.4];
  
  figure;
  
  for i=1:3
    t=0:0.01:w(i);
    y=h(i)*square(2*pi*f(i)*t);
    subplot(3,1,i);
    plot(t,y,'r','LineWidth',3);
    grid on;
    xlim([-2 2.5]);
    ylim([-.6 2.5]);
    title(sprintf...
   ('pulse of %2.1f sec width and amplitude=%2.1f',...
    w(i),h(i)));
  
  end
  

PIC



63 Generate and plot train of pulses of different width and amplitude

Problem: Generate a pulse train of pulses of certain width and height.



Mathematica

  Remove["Global*"]
  
  w = {0.25,0.125};
  h = {1,1.5};
  plots = MapThread[Plot[#2 SquareWave[#1 x],
           {x,-10,10},
    PlotRange->{All,{-2,2}},
    Frame->True,
    Axes->False,
    PlotStyle->{Thick,Red},
    FrameLabel->{{"y(t)",None},
     {"t","pulse of "<>ToString[1/(2*#1)]<>
     " sec width and amplitude"<>ToString[#2]}},
    ImageSize->200,
    AspectRatio->1,
    ExclusionsStyle->Dotted]&,{w,h}];
  
  Grid[{plots},Frame->All]

PIC

Matlab

  clear all; close all;
  f = [1 0.5];
  h = [1 1.5];
  t = -4:0.01:4;
  
  figure;
  
  for i=1:2
    y = h(i)*square(2*pi*f(i)*t);
    subplot(2,1,i);
    plot(t,y,'r','LineWidth',3);
    grid on;
    xlim([-2 2.5]);
    ylim([-2 2.5]);
    title(sprintf...
   ('pulse of %2.1f sec width and amplitude=%2.1f',..
     1/(2*f(i)),h(i)));
  end
  

PIC



64 Find the discrete Fourier transform of a real sequence of numbers

Given the sequence of numbers x (n ) = [1,2,3 ]  , Find          N∑−1         2π
X  (k) =    x (m )e−jN mk
         m=0  where N  is the length of the data sequence x (m )  and k =  0⋅⋅⋅N  − 1



Mathematica

  Chop[Fourier[{1, 2, 3},
        FourierParameters ->{1, -1}]]
  
  Out[5]= {6.,
           -1.5 + 0.8660254037844386*I,
           -1.5 - 0.8660254037844386*I}

Matlab

  fft([1,2,3])'
  
  ans =
      6.0000
     -1.5000 - 0.8660i
     -1.5000 + 0.8660i



Ada I posted the code below on comp.lang.ada on June 8,2010. Thanks to Ludovic Brenta for making improvement to the Ada code.

  --
  -- dtf.adb, compiled with GNAT 4.3.4 20090804 (release) 1
  -- under CYGWIN 1.7.5
  -- $ gnatmake dft.adb
  -- gcc -c dft.adb
  -- gnatbind -x dft.ali
  -- gnatlink dft.ali
  -- $ ./dft.exe
  -- ( 6.00000E+00, 0.00000E+00)
  -- (-1.50000E+00, 8.66026E-01)
  -- (-1.50000E+00,-8.66025E-01)
  
  with Ada.Text_IO; use Ada.Text_IO;
  with Ada.Numerics.Complex_Types; use  Ada.Numerics.Complex_Types;
  
  with Ada.Numerics; use  Ada.Numerics;
  
  with Ada.Numerics.Complex_Elementary_Functions;
  use  Ada.Numerics.Complex_Elementary_Functions;
  with Ada.Complex_Text_IO; use Ada.Complex_Text_IO;
  
  procedure dft is
       N : constant := 3; -- named number, no conversion to Float needed
       X : array(0 .. N-1) of Complex  := (others=>(0.0,0.0));
       data : constant array(0 .. N-1) of float :=(1.0,2.0,3.0);
       Two_Pi_Over_N : constant := 2 * Pi / N;
        -- named number, outside the loop, like in ARM 3.3.2(9)
  begin
       FOR k in X'range LOOP
           FOR m in data'range LOOP
               X(k) := X(k) + data(m)*exp(-J*Two_Pi_Over_N * float(m*k));
           END LOOP;
           put(X(k)); new_line;
       END LOOP;
  end dft;

Fortran Thanks to Vincent Lafage for making improvment to the Fortran code.

  ! dtf.f90, compiled with GCC 4.3.4
  ! under CYGWIN 1.7.5
  !  $ gfortran -Wall -Wsurprising -Wconversion dft.f90
  !  $ ./a.exe
  ! (  6.0000000    ,  0.0000000    )
  ! ( -1.4999999    , 0.86602557    )
  ! ( -1.5000005    ,-0.86602497    )
  !  $
  
  PROGRAM dft
  
    IMPLICIT NONE
  
    INTEGER, parameter :: N = 3
    COMPLEX, parameter :: J =(0.0,1.0)
    REAL,    parameter :: Pi = ACOS(-1.0)
    INTEGER                   :: k,m
    COMPLEX, dimension(0:N-1) :: X
    REAL,    dimension(0:N-1) :: data=(/1.0,2.0,3.0/)
    REAL,    parameter        :: Two_Pi_Over_N = 2.0*Pi/real(N)
  
  DO k = lbound(X, 1), ubound(X, 1)
     X(k)=(0.0,0.0)
     DO m = lbound(data, 1), ubound(data, 1)
        X(k) = X(k) + complex(data(m),0.0)                   &
               * EXP(-J*complex(Two_Pi_Over_N*real(m*k),0.0))
     END DO
     print *,X(k)
  END DO
  
  END PROGRAM dft

65 Find eAt  where A  is a matrix



Mathematica

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

PIC 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 |A  − λI|

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

PIC 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)|
  +-                                             -+



66 How to do some common operations on vectors and matrices?

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

66.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
  c =
      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



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



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