home
PDF letter size
PDF legal size

General report

How to solve Basic Engineering and Mathematics Problems

using Mathematica, Matlab and Maple

Nasser M. Abbasi

June 11, 2020   Compiled on June 11, 2020 at 6:48pm

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

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

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

1 Control systems, Linear systems, transfer functions, state space related problems
 1.1 Creating tf and state space and different Conversion of forms
 1.2 Obtain the step response of an LTI from its transfer function
 1.3 plot the impulse and step responses of a system from its transfer function
 1.4 Obtain the response of a transfer function for an arbitrary input
 1.5 Obtain the poles and zeros of a transfer function
 1.6 Generate Bode plot of a transfer function
 1.7 How to check that state space system \(x'=Ax+Bu\) is controllable?
 1.8 Obtain partial-fraction expansion
 1.9 Obtain Laplace transform for a piecewise functions
 1.10 Obtain Inverse Laplace transform of a transfer function
 1.11 Display the response to a unit step of an under, critically, and over damped system
 1.12 View steady state error of 2nd order LTI system with changing undamped natural frequency
 1.13 Show the use of the inverse Z transform
 1.14 Find the Z transform of sequence x(n)
 1.15 Sample a continuous time system
 1.16 Find closed loop transfer function from the open loop transfer function for a unity feedback
 1.17 Compute the Jordan canonical/normal form of a matrix A
 1.18 Solve the continuous-time algebraic Riccati equation
 1.19 Solve the discrete-time algebraic Riccati equation
 1.20 Display impulse response of H(z) and the impulse response of its continuous time approximation H(s)
 1.21 Find the system type given an open loop transfer function
 1.22 Find the eigenvalues and eigenvectors of a matrix
 1.23 Find the characteristic polynomial of a matrix
 1.24 Verify the Cayley-Hamilton theorem that every matrix is zero of its characteristic polynomial
 1.25 How to check for stability of system represented as a transfer function and state space
 1.26 Check continuous system stability in the Lyapunov sense
 1.27 Given a closed loop block diagram, generate the closed loop Z transform and check its stability
 1.28 Determine the state response of a system to only initial conditions in state space
 1.29 Determine the response of a system to only initial conditions in state space
 1.30 Determine the response of a system to step input with nonzero initial conditions
 1.31 Draw the root locus from the open loop transfer function
 1.32 Find \(e^{A t}\) where \(A\) is a matrix
 1.33 Draw root locus for a discrete system
 1.34 Plot the response of the inverted pendulum problem using state space
 1.35 How to build and connect a closed loop control systems and show the response?
 1.36 Compare the effect on the step response of a standard second order system as \(\zeta \) changes
 1.37 Plot the dynamic response factor \(R_{d}\) of a system as a function of \(r=\frac{\omega }{\omega _{n}}\) for different damping ratios
 1.38 How to find closed loop step response to a plant with a PID controller?
 1.39 How to make Nyquist plot?
2 Linear algebra, Linear solvers, common operations on vectors and matrices
 2.1 Multiply matrix with a vector
 2.2 Insert a number at specific position in a vector or list
 2.3 Insert a row into a matrix
 2.4 Insert a column into a matrix
 2.5 Build matrix from other matrices and vectors
 2.6 Generate a random 2D matrix from uniform (0 to 1) and from normal distributions
 2.7 Generate an n by m zero matrix
 2.8 Rotate a matrix by 90 degrees
 2.9 Generate a diagonal matrix with given values on the diagonal
 2.10 Sum elements in a matrix along the diagonal
 2.11 Find the product of elements in a matrix along the diagonal
 2.12 Check if a Matrix is diagonal
 2.13 Find all positions of elements in a Matrix that are larger than some value
 2.14 Replicate a matrix
 2.15 Find the location of the maximum value in a matrix
 2.16 Swap 2 columns in a matrix
 2.17 Join 2 matrices side-by-side and on top of each others
 2.18 Copy the lower triangle to the upper triangle of a matrix to make symmetric matrix
 2.19 extract values from matrix given their index
 2.20 Convert \(N\) by \(M\) matrix to a row of length \(N M\)
 2.21 find rows in a matrix based on values in different columns
 2.22 Select entries in one column based on a condition in another column
 2.23 Locate rows in a matrix with column being a string
 2.24 Remove set of rows and columns from a matrix at once
 2.25 Convert list of separated numerical numbers to strings
 2.26 Obtain elements that are common to two vectors
 2.27 Sort each column (on its own) in a matrix
 2.28 Sort each row (on its own) in a matrix
 2.29 Sort a matrix row-wise using first column as key
 2.30 Sort a matrix row-wise using non-first column as key
 2.31 Replace the first nonzero element in each row in a matrix by some value
 2.32 Perform outer product and outer sum between two vector
 2.33 Find the rank and the bases of the Null space for a matrix A
 2.34 Find the singular value decomposition (SVD) of a matrix
 2.35 Solve \(Ax=b\)
 2.36 Find all nonzero elements in a matrix
 2.37 evaluate f(x) on a vector of values
 2.38 generates equally spaced N points between \(x_1\) and \(x_2\)
 2.39 evaluate and plot a \(f(x,y)\) on 2D grid of coordinates
 2.40 Find determinant of matrix
 2.41 Generate sparse matrix with \(n\) by \(n\) matrix repeated on its diagonal
 2.42 Generate sparse matrix for the tridiagonal representation of second difference operator in 1D
 2.43 Generate sparse matrix for the Laplacian differential operator \(\nabla ^{2}u\) for 2D grid
 2.44 Generate sparse matrix for the Laplacian differential operator \(\nabla ^{2}u\) for 3D grid
 2.45 Generate the adjugate matrix for square matrix
 2.46 Multiply each column by values taken from a row
 2.47 extract submatrix from a larger matrix by removing row/column
 2.48 delete one row from a matrix
 2.49 delete one column from a matrix
 2.50 generate random matrix so that each row adds to 1
 2.51 generate random matrix so that each column adds to 1
 2.52 sum all rows in a matrix
 2.53 sum all columns in a matrix
 2.54 find in which columns values that are not zero
 2.55 How to remove values from one vector that exist in another vector
 2.56 How to find mean of equal sized segments of a vector
 2.57 find first value in column larger than some value and cut matrix from there
 2.58 make copies of each value into matrix into a larger matrix
 2.59 repeat each column of matrix number of times
 2.60 How to apply a function to each value in a matrix?
 2.61 How to sum all numbers in a list (vector)?
 2.62 How to find maximum of each row of a matrix?
 2.63 How to find maximum of each column of a matrix?
 2.64 How to add the mean of each column of a matrix from each column?
 2.65 How to add the mean of each row of a matrix from each row?
 2.66 Find the different norms of a vector
 2.67 Check if a matrix is Hermite
 2.68 Obtain the LU decomposition of a matrix
 2.69 Linear convolution of 2 sequences
 2.70 Circular convolution of two sequences
 2.71 Linear convolution of 2 sequences with origin at arbitrary position
 2.72 Visualize a 2D matrix
 2.73 Find the cross correlation between two sequences
 2.74 Find orthonormal vectors that span the range of matrix A
 2.75 Solve \(A x= b\) and display the solution
 2.76 Determine if a set of linear equations \(A x= b\) has a solution and what type of solution
 2.77 Given a set of linear equations automatically generate the matrix \(A\) and vector \(b\) and solve \(A x=b\)
 2.78 Convert a matrix to row echelon form and to reduced row echelon form
 2.79 Convert 2D matrix to show the location and values
 2.80 Find rows in matrix with zeros in them, and then remove the zeros
 2.81 How to apply a function to two lists at the same time?
 2.82 How to apply a function to two lists are the same time, but with change to entries?
 2.83 How to select all primes numbers from a list?
 2.84 How to collect result inside a loop when number of interation is not known in advance?
 2.85 How flip an array around?
 2.86 How to divide each element by its position in a list?
3 signal processing, Image processing, graphics, random numbers
 3.1 Generate and plot one pulse signal of different width and amplitude
 3.2 Generate and plot train of pulses of different width and amplitude
 3.3 Find the discrete Fourier transform of a real sequence of numbers
 3.4 Generate uniform distributed random numbers
 3.5 Obtain Fourier Series coefficients for a periodic function
 3.6 Determine and plot the CTFT (continuous time Fourier Transform) for continuous time function
 3.7 Determine the DTFT (Discrete time Fourier Transform) for discrete time function
 3.8 Determine the Inverse DTFT (Discrete time Fourier Transform)
 3.9 Use FFT to display the power spectrum of the content of an audio wav file
 3.10 Plot the power spectral density of a signal
 3.11 Display spectrum of 2D image
 3.12 Obtain the statistical maximum likelihood estimates (MLE) of probability distributions
 3.13 Make a histogram of data sampled from some probability distribution
 3.14 apply a filter on 1D numerical data (a vector)
 3.15 apply an averaging Laplacian filter on 2D numerical data (a matrix)
 3.16 How to find cummulative sum
 3.17 How to find the moving average of a 1D sequence?
 3.18 How to select N random values from a set of numbers?
 3.19 How to sample a sin signal and plot it?
 3.20 How find the impulse response of a difference equation?
4 Differential, PDE solving, integration, numerical and analytical solving of equations
 4.1 Generate direction field plot of a first order differential equation
 4.2 Solve the 2-D Laplace PDE for a rectangular plate with Dirichlet boundary conditions
 4.3 Solve homogeneous 1st order ODE, constant coefficients and initial conditions
 4.4 Solve homogeneous 2nd order ODE with constant coefficients
 4.5 Solve non-homogeneous 2nd order ODE, constant coefficients
 4.6 Solve homogeneous 2nd order ODE, constant coefficients (BVP)
 4.7 Solve the 1-D heat partial differential equation (PDE)
 4.8 Show the effect of boundary/initial conditions on 1-D heat PDE
 4.9 Find particular and homogenous solution to undetermined system of equations
 4.10 Plot the constant energy levels for a nonlinear pendulum
 4.11 Solve numerically the ODE \(u^{''''}+u=f\) using point collocation method
 4.12 How to numerically solve a set of non-linear equations?
 4.13 Solve 2nd order ODE (Van Der Pol) and generate phase plot
 4.14 How to numerically solve Poisson PDE on 2D using Jacobi iteration method?
 4.15 How to solve BVP second order ODE using finite elements with linear shape functions and using weak form formulation?
 4.16 How to solve Poisson PDE in 2D using finite elements methods using rectangular element?
 4.17 How to solve Poisson PDE in 2D using finite elements methods using triangle element?
 4.18 How to solve wave equation using leapfrog method?
 4.19 Numerically integrate \(f(x)\) on the real line
 4.20 Numerically integrate \(f(x,y)\) in 2D
 4.21 How to solve set of differential equations in vector form
 4.22 How to implement Runge-Kutta to solve differential equations?
 4.23 How to differentiate treating a combined expression as single variable?
 4.24 How to solve Poisson PDE in 2D with Neumann boundary conditions using Finite Elements
5 plotting how to, ellipse, 3D
 5.1 Generate a plot using x and y data
 5.2 Plot the surface described by 2D function
 5.3 Find the point of intersection of 3 surfaces
 5.4 How to draw an ellipse?
 5.5 How to plot a function on 2D using coordinates of grid locations?
 5.6 How to make table of \(x,y\) values and plot them
 5.7 How to make contour plot of \(f(x,y)\) ?
6 Some symbolic operations
 6.1 How to find all exponents of list of expressions?

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

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

1.1 Creating tf and state space and different Conversion of forms

  1.1.1 Create continuous time transfer function given the poles, zeros and gain
  1.1.2 Convert transfer function to state space representation
  1.1.3 Create a state space representation from A,B,C,D and find the step response
  1.1.4 Convert continuous time to discrete time transfer function, gain and phase margins
  1.1.5 Convert differential equation to transfer functions and to state space
  1.1.6 Convert from continuous transfer function to discrete time transfer function
  1.1.7 Convert a Laplace transfer function to an ordinary differential equation

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

   1.1.1.1 Mathematica
   1.1.1.2 Matlab
   1.1.1.3 Maple

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

1.1.1.1 Mathematica



Out[30]= TransferFunctionModel[{
     {{5*(1 + s)*(2 + s)}},
     s*(4 + s)*(6 + s)}, s]



 

1.1.1.2 Matlab



num/den =
    5 s^2 + 15 s + 10
   -------------------
   s^3 + 10 s^2 + 24 s



 

1.1.1.3 Maple



\[ \text{tf} = \left [{\begin{array}{c}{\frac{5\,{s}^{2}+15\,s+10}{{s}^{3}+10\,{s}^{2}+24\,s}}\end{array}} \right ] \]




\[ 5*s^2+15*s+10 \]




\[ s*(s^2+10*s+24) \]



 

1.1.2 Convert transfer function to state space representation

   1.1.2.1 problem 1
   1.1.2.2 Mathematica
   1.1.2.3 Matlab
   1.1.2.4 Maple
   1.1.2.5 problem 2

1.1.2.1 problem 1

Problem: Find the state space representation for the continuous time system defined by the transfer function \[ H(s)=\frac{5s}{s^{2}+4s+25}\]

1.1.2.2 Mathematica



pict




pict



 

1.1.2.3 Matlab



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



 

1.1.2.4 Maple



\[ \left [{\begin{array}{cc} 0 & 1 \\ -25 & -4 \end{array}} \right ] \]




\[ \left [{\begin{array}{c} 0\\ 1\end{array}} \right ] \]




\[ \left [{\begin{array}{cc} 0&5\end{array}} \right ] \]




\[ \left [{\begin{array}{c} 0\end{array}} \right ] \]



 

1.1.2.5 problem 2

Problem: Given the continuous time S transfer function defined by \[ H(s)=\frac{1+s}{s^{2}+s+1}\] convert to state space and back to transfer function.



Mathematica


\[ \left ({\begin{array}{cc|c|c} 0 & 1 & 0 \\ -1 & -1 & 1 \\ \hline 1 & 1 & 0 \\ \end{array}} \right )_{} \]


\[ \frac{s+1}{s^2+s+1} \]



 


Matlab









 

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

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




Mathematica


pict



 


Matlab


pict



 


Maple


pict



 

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

Problem: Compute the gain and phase margins of the open-loop discrete linear time system sampled from the continuous time S transfer function defined by \[ H(s)=\frac{1+s}{s^{2}+s+1}\] Use sampling period of 0.1 seconds.



Mathematica


pict




pict








pict



 


Matlab






pict



 

1.1.5 Convert differential equation to transfer functions and to state space

Problem: Obtain the transfer and state space representation for the differential equation \[ 3\frac{d^{2}y}{dt^{2}}+2\frac{dy}{dt}+y\left ( t\right ) = u(t) \]



Mathematica


pict




pict



 


Matlab









 


Maple


\[ \frac{1}{3{s}^{2}+2\,s+1} \]




\[ \left \{ \left [{\begin{array}{cc} 0&1\\ \noalign{\medskip }-1/3&-2/3 \end{array}} \right ] , \left [{\begin{array}{c} 0\\ \noalign{\medskip } 1\end{array}} \right ] , \left [{\begin{array}{cc} 1/3&0\end{array}} \right ] , \left [{\begin{array}{c} 0\end{array}} \right ] \right \} \]



 

1.1.6 Convert from continuous transfer function to discrete time transfer function

Problem: Convert \[ H\left (s\right ) = \frac{1}{s^2+10 s+10} \]

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



Mathematica


pict



 


Matlab





 

1.1.7 Convert a Laplace transfer function to an ordinary differential equation

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



Mathematica





 

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

Problem: Find the unit step response for the continuous time system defined by the transfer function \[ H(s) = \frac{25}{s^{2}+4s+25} \]



Mathematica


pict



 


Matlab


pict



 


Maple


pict



 

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

Problem: Find the impulse and step responses for the continuous time system defined by the transfer function \[ H\left ( s\right ) =\frac{1}{s^{2}+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










pict



 


Matlab


pict



 


Maple

Using Maple DynamicSystems


pict



Using Laplace transform method:


pict



 

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

Problem: Find the response for the continuous time system defined by the transfer function \[ H(s)=\frac{1}{s^{2}+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


pict



 


Matlab


pict



 


Maple


pict



Using DynamicSystem package


pict



 

1.5 Obtain the poles and zeros of a transfer function

Problem: Find the zeros, poles, and gain for the continuous time system defined by the transfer function \[ H(s)=\frac{25}{s^{2}+4s+25}\]



Mathematica









 


Matlab





 


Maple





 

1.6 Generate Bode plot of a transfer function

Problem: Generate a Bode plot for the continuous time system defined by the transfer function \[ H(s)=\frac{5s}{s^{2}+4s+25}\]



Mathematica


pict



 


Matlab


pict



 


Maple


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


pict



 

1.7 How to check that state space system \(x'=Ax+Bu\) is controllable?

A system described by \begin{align*} x' &= Ax+Bu \\ y &= Cx+Du \end{align*}

Is controllable if for any initial state \(x_0\) and any final state \(x_f\) there exist an input \(u\) which moves the system from \(x_0\) to \(x_f\) in finite time. Only the matrix \(A\) and \(B\) are needed to decide on controllability. If the rank of \[ [B \> AB\> A^2B\> \ldots \> A^{n-1}B] \] is \(n\) which is the number of states, then the system is controllable. Given the matrix \[ A=\left ({\begin{array}{cccc} 0&1&0&0\\ 0&0&-1&0\\ 0&0&0&1\\ 0&0&5&0 \end{array}} \right ) \] And \[ B=\left ({\begin{array}{c} 0\\ 1\\ 0\\ -2 \end{array}} \right ) \]



Mathematica


\[ \left ({\begin{array}{cccc} 0 & 1 & 0 & 2 \\ 1 & 0 & 2 & 0 \\ 0 & -2 & 0 & -10 \\ -2 & 0 & -10 & 0 \\ \end{array}} \right ) \]




True




4



 


Matlab






4



 


Maple


\[ \left [{\begin{array}{cccc} 0&1&0&2\\ \noalign{\medskip }1&0&2&0 \\ \noalign{\medskip }0&-2&0&-10\\ \noalign{\medskip }-2&0&-10&0 \end{array}} \right ] \]




true




true




4



 

1.8 Obtain partial-fraction expansion

Problem: Given the continuous time S transfer function defined by \[ H(s)=\frac{s^{4}+8s^{3}+16s^{2}+9s+9}{s^{3}+6s^{2}+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





 


Matlab





 


Maple


\[ s+2- \frac{7}{s+2}-{\frac{9}{2\,s+6}}+{\frac{9}{2\,s+2 }} \]


\[ [s,2,\frac{7}{s+2},-{\frac{9}{2\,s+6}},{\frac{9}{2\,s+2}}] \]



 

1.9 Obtain Laplace transform for a piecewise functions

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

pict

Function f(t) to obtain its Laplace transform

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



Mathematica





 


Matlab





 


Maple

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


\[{\frac{1-{{e}^{-sT}}}{{s}^{2}}} \]



 

1.10 Obtain Inverse Laplace transform of a transfer function

Problem: Obtain the inverse Laplace transform for the function \[ H(s)=\frac{s^{4}+5s^{3}+6s^{2}+9s+30}{s^{4}+6s^{3}+21s^{2}+46s+30}\]

Mathematica


\[ \delta (t)+\left (\frac{1}{234}+\frac{i}{234}\right ) e^{(-1-3 i) t} \left ((73+326 i) e^{6 i t}+(-326-73 i)\right )-\frac{3 e^{-3 t}}{26}+\frac{23 e^{-t}}{18} \]

Matlab





Maple


\[ Dirac \left ( t \right ) -{\frac{3\,{{\rm e}^{-3\,t}}}{26}}+{ \frac{ \left ( -506\,\cos \left ( 3\,t \right ) -798\,\sin \left ( 3\,t \right ) +299 \right ){{\rm e}^{-t}}}{234}} \]

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

Problem: Obtain unit step response of the second order system given by the transfer function \[ H\left ( s\right ) =\frac{\omega _{n}^{2}}{s^{2}+2\xi \omega _{n}s+\omega _{n}^{2}}\] in order to illustrate the response when the system is over, under, and critically damped. use \(\omega _{n}=1\) and change \(\xi \) over a range of values that extends from under damped to over damped.

Mathematica


pict

Matlab


pict

Maple


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


pict

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

  1.12.1 Mathematica

Problem: Given the transfer function \[ H\left ( s\right ) =\frac{\omega _{n}^{2}}{s^{2}+2\xi \omega _{n}s+\omega _{n}^{2}}\]

Display the output and input on the same plot showing how the steady state error changes as the un damped natural frequency \(\omega _{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 \(\xi =0.707\) and change \(\omega _{n}\) from \(0.2\) to \(1.2\).

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

1.12.1 Mathematica

   1.12.1.1 ramp input
   1.12.1.2 step input

1.12.1.1 ramp input

pict

1.12.1.2 step input

pict

1.13 Show the use of the inverse Z transform

  1.13.1 example 1
  1.13.2 example 2

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

1.13.1 example 1

Problem: Given \[ F\left (z\right ) =\frac{z}{z-1}\] find \(x[n]=F^{-1}\left (z\right )\) which is the inverse Ztransform.



Mathematica


pict



 


Matlab


pict



 

1.13.2 example 2

Problem: Given \[ F\left ( z\right ) =\frac{5z}{\left ( z-1\right ) ^{2}}\] find \(x[n]=F^{-1}\left ( z\right ) \)

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



Mathematica


Inverse Z is 5 n


pict



 


Matlab


pict



 

1.14 Find the Z transform of sequence x(n)

  1.14.1 example 1
  1.14.2 example 2

1.14.1 example 1

Find the Z transform for the unit step discrete function

Given the unit step function \(x[n]=u[n]\) defined as \(x=\{1,1,1,\cdots \}\,\ \) for \(n\geq 0\,\), find its Z transform.



Mathematica





 


Matlab





 

1.14.2 example 2

Find the Z transform for \(x[n]=\left ( \frac{1}{3}\right ) ^{n}u\left ( n\right ) +\left ( 0.9\right ) ^{n-3}u\left ( n\right ) \)



Mathematica





 


Matlab





 

1.15 Sample a continuous time system

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

pict

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

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

Mathematica


pict


pict


pict

Matlab




pict

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

Problem: Given \[ L(s)=\frac{s}{(s+4)(s+5)}\] as the open loop transfer function, how to find \(G(s)\), the closed loop transfer function for a unity feedback?

pict



Mathematica


pict




pict



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


\[ \frac{s}{s^2+10 s+20} \]



 


Matlab





 

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



Mathematica


\[ \left ({\begin{array}{cccccc} 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 \\ \end{array}} \right ) \]




\[ \left ({\begin{array}{cccccc} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & 1 & 0 & 0 & 0 \\ 0 & 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 2 \\ \end{array}} \right ) \]



 


Matlab





 

1.18 Solve the continuous-time algebraic Riccati equation

Problem: Solve for \(X\) in the Riccati equation \[ A^{\prime }X+XA-XBR^{-1}B^{\prime }X+C^{\prime }C=0 \] given \begin{align*} A & =\begin{pmatrix} -3 & 2\\ 1 & 1 \end{pmatrix} \\ B & =\begin{pmatrix} 0\\ 1 \end{pmatrix} \\ C & =\begin{pmatrix} 1 & -1 \end{pmatrix} \\ R & =3 \end{align*}



Mathematica


\[ \left ({\begin{array}{cc} 0.589517 & 1.82157 \\ 1.82157 & 8.81884 \\ \end{array}} \right ) \]



 


Matlab





 

1.19 Solve the discrete-time algebraic Riccati equation

Problem: Given a continuous-time system represented by a transfer function \[ \frac{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^{\prime }X+XA-XBR^{-1}B^{\prime }X+C^{\prime }C=0 \]

Let \(R=[3]\).



Mathematica


pict




pict




\[ \left ({\begin{array}{cc} 0.671414 & -0.977632 \\ -0.977632 & 2.88699 \\ \end{array}} \right ) \]



 


Matlab













 

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

Plot the impulse response of \(H(z)=z/(z^2-1.4 z+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


pict

Find its impulse response



approximate to continuous time, use ZeroOrderHold


pict

Find the impulse response of the continuous time system


\[ -1.25559 e^{-0.693147 t} (-13.3012 \theta (t) \sin (0.283794 t)-1. \theta (t) \cos (0.283794 t)) \]


pict

Plot the impulse response of the discrete system


pict

Plot the impulse response of the continuous system


pict

Plot both responses on the same plot


pict

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


pict

Matlab


pict

1.21 Find the system type given an open loop transfer function

Problem: Find the system type for the following transfer functions

  1. \(\frac{s+1}{s^{2}-s}\)
  2. \(\frac{s+1}{s^{3}-s^{2}}\)
  3. \(\frac{s+1}{s^{5}}\)

To find the system type, the transfer function is put in the form \(\frac{k\sum _{i}\left ( s-s_{i}\right ) }{s^{M}\sum _{j}\left ( s-s_{j}\right ) }\), 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 \(s^{1}\left ( s-1\right )\) and the second system has type 2 since the denominator can be written as \(s^{2}\left ( s-1\right ) \) and the third system has type 5. The following computation determines the type



Mathematica


Out[171]= 1




Out[173]= 2




Out[175]= 5



 


Matlab













 

1.22 Find the eigenvalues and eigenvectors of a matrix

Problem, given the matrix \[ \left ({\begin{array} [c]{ccc}1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 9 \end{array}} \right ) \] Find its eigenvalues and eigenvectors.



Mathematica


\[ \left ({\begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ \end{array}} \right ) \]




\[ \left \{\frac{3}{2} \left (5+\sqrt{33}\right ),\frac{3}{2} \left (5-\sqrt{33}\right ),0\right \} \] \[ \{16.1168,-1.11684,0.\} \]




\[ \left ({\begin{array}{ccc} -\frac{-15-\sqrt{33}}{33+7 \sqrt{33}} & \frac{4 \left (6+\sqrt{33}\right )}{33+7 \sqrt{33}} & 1 \\ -\frac{15-\sqrt{33}}{-33+7 \sqrt{33}} & \frac{4 \left (-6+\sqrt{33}\right )}{-33+7 \sqrt{33}} & 1 \\ 1 & -2 & 1 \\ \end{array}} \right ) \] \[ \left ({\begin{array}{ccc} 0.283349 & 0.641675 & 1. \\ -1.28335 & -0.141675 & 1. \\ 1. & -2. & 1. \\ \end{array}} \right ) \]



 


Matlab

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





 

1.23 Find the characteristic polynomial of a matrix

\[ \left ({\begin{array} [c]{ccc}1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 0 \end{array}} \right ) \]



Mathematica


\[ -x^3+6 x^2+72 x+27 \]



 


Matlab

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





 

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

Problem, given the matrix \[ \left ({\begin{array} [c]{cc}1 & 2\\ 3 & 2 \end{array}} \right ) \] 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


\(x^2-3 x-4\)




\[ \left ({\begin{array}{cc} 0 & 0 \\ 0 & 0 \\ \end{array}} \right ) \]



Another way is as follows


\[ \left ({\begin{array}{cc} 0 & 0 \\ 0 & 0 \\ \end{array}} \right ) \]



 

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



Matlab





 

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

  1.25.1 Checking stability using transfer function poles
  1.25.2 Checking stability using state space A matrix

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 \[ sys=\frac{5s}{s^{2}+4s+25}\]

1.25.1 Checking stability using transfer function poles



Mathematica













 


Matlab









 

1.25.2 Checking stability using state space A matrix



Mathematica

















 


Matlab













 

1.26 Check continuous system stability in the Lyapunov sense

Problem: Check the stability (in Lyapunov sense) for the state coefficient matrix \[ A=\begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ -1 & -2 & -3 \end{bmatrix} \]

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

We must transpose the matrix \(A\) when calling these functions, since the Lyapunov equation is defined as \(A^TP+PA=-Q\) and this is not how the software above defines them. By simply transposing the \(A\) matrix when calling them, then the result will be correct.



Mathematica


\[ \left ({\begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ -1 & -2 & -3 \\ \end{array}} \right ) \]




\[ \left ({\begin{array}{ccc} 2.3 & 2.1 & 0.5 \\ 2.1 & 4.6 & 1.3 \\ 0.5 & 1.3 & 0.6 \\ \end{array}} \right ) \]




\(\{6.18272,1.1149,0.202375\}\)



 


Matlab









 


Maple


\[ \left [{\begin{array}{c} 6.18272045921436+ 0.0\,i \\ \noalign{\medskip } 1.11490451203192+ 0.0\,i\\ \noalign{\medskip } 0.202375028753723+ 0.0\,i\end{array}} \right ] \]



 

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

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

pict

System block diagram.

Mathematica


pict


pict


pict


pict

Now generate unit step response


pict


pict


\[{\begin{array}{c} \{0.543991\, -0.325556 i,0.543991\, +0.325556 i\} \\ \end{array}} \]


\[ \left ({\begin{array}{c} \{0.633966,0.633966\} \\ \end{array}} \right ) \]

Poles are inside the unit circle, hence stable.

Matlab




pict



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

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



Mathematica


pict




pict



 


Matlab


pict



 

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

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



Mathematica


pict




pict



 


Matlab


pict



 

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

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



Mathematica


pict




pict



 


Matlab


pict



 

1.31 Draw the root locus from the open loop transfer function

Problem: Given \(L(s)\), the open loop transfer function, draw the root locus. Let \[ L(s)=\frac{s^2+2 s+4}{s(s+4)(s+6)(s^2+1.4 s+1)}\]

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



Mathematica


pict



 


Matlab


pict



 

1.32 Find \(e^{A t}\) where \(A\) is a matrix

Mathematica


\[ \left ({\begin{array}{cc} -e^{-2 t}+2 e^{-t} & -e^{-2 t}+e^{-t} \\ 2 e^{-2 t}-2 e^{-t} & 2 e^{-2 t}-e^{-t} \\ \end{array}} \right ) \]

Now verify the result by solving for \(e^{At}\) 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-\lambda I|\)


\[ \left ({\begin{array}{cc} -\lambda & 1 \\ -2 & -\lambda -3 \\ \end{array}} \right ) \]


\[ \lambda ^2+3 \lambda +2 \]






\[ \left \{\text{b0}\to e^{-2 t} \left (2 e^t-1\right ),\text{b1}\to e^{-2 t} \left (e^t-1\right )\right \} \]


\[ e^{-2 t} \left (2 e^t-1\right ) \]


\[ e^{-2 t} \left (e^t-1\right ) \]


\[ \left ({\begin{array}{cc} e^{-2 t} \left (-1+2 e^t\right ) & e^{-2 t} \left (-1+e^t\right ) \\ -2 e^{-2 t} \left (-1+e^t\right ) & -e^{-2 t} \left (-2+e^t\right ) \\ \end{array}} \right ) \] The answer is the same given by Mathematica’s command MatrixExp[]

Matlab





1.33 Draw root locus for a discrete system

Problem: Given the open loop for a continuous time system as \[ sys=\frac{5s+1}{s^{2}+2s+3}\]

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



Mathematica


pict




pict



 


Matlab


pict



 

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

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

pict

Mathematica


\[ \left ({\begin{array}{cccc} 0 & 1 & 0 & 0 \\ 0 & 0 & -\frac{g m}{M} & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & \frac{g (m+M)}{L M} & 0 \\ \end{array}} \right ) \]


\[ \left ({\begin{array}{c} 0 \\ \frac{1}{M} \\ 0 \\ -\frac{1}{L M} \\ \end{array}} \right ) \]


\[ \left ({\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ \end{array}} \right ) \]


pict

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

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

Here is example of obtaining the analytical expression of the response


\({\begin{array}{l} e^{-6.41561 t} (0.0238095 e^{6.41561 t} t^2 \theta (t)+0.000115693 e^{3.2078 t} \theta (t)-0.000231385 e^{6.41561 t} \theta (t)+0.000115693 e^{9.62341 t} \theta (t))\\,e^{-6.41561 t} (-0.00242954 e^{3.2078 t} \theta (t)+0.00485909 e^{6.41561 t} \theta (t)-0.00242954 e^{9.62341 t} \theta (t)) \end{array}} \)


pict


pict

Matlab


pict


pict

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

  1.35.1 example 1, single input, single output closed loop

1.35.1 example 1, single input, single output closed loop

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

pict

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

Mathematica


pict

Matlab


Another way to do the above is


pict

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

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

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

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

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

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

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

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

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

Mathematica


pict

Matlab


pict

1.37 Plot the dynamic response factor \(R_{d}\) of a system as a function of \(r=\frac{\omega }{\omega _{n}}\) for different damping ratios

Problem: Plot the standard curves showing how the dynamic response \(R_{d}\) changes as \(r=\frac{\omega }{\omega _{n}}\) changes. Do this for different damping ratio \(\xi \). Also plot the phase angle.

These plots are result of analysis of the response of a second order damped system to a harmonic loading. \(\omega \) is the forcing frequency and \(\omega _{n}\) is the natural frequency of the system.



Mathematica


pict




pict



 

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

Find and plot the step response of the plant \(\frac{1}{s^2+2s+1}\) connected to a PID controller with \(P=10, I=3.7, D=0.7\). Use negative closed loopback.



Mathematica


pict



 


Matlab


pict



 

1.39 How to make Nyquist plot?

  1.39.1 Example 1
  1.39.2 Example 2

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

Here are two examples.

1.39.1 Example 1

Given \(G(s)=\frac{s}{1-0.2s}\) generate Nyquist plot.



Mathematica


pict



 


Matlab


pict



 

1.39.2 Example 2

Given \(G(s)=\frac{5(1-0.5s)}{s(1+0.1s)(1-0.25s)}\) generate Nyquist plot.



Mathematica


pict



 


Matlab






pict



 

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




pict



 

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

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

Mathematica and Matlab allow 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). So a Mathematica matrix is stored in a list of lists. This is similar in a way to Matlab cell data structure, since each raw can have different length. In a standard matrix each row must have the same length.

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 \[ A = \left ({\begin{array}{cccc} 1 & 2 & 3 & 4 \\ 2 & 3 & 1 & 5 \\ 5 & 6 & 7 & 2 \end{array}} \right )\] Then the result of using the find() command in Matlab is



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


which gives


Mathematica searched row-wise.

I found Mathematica for matrix operations takes more time to get familiar with compared to Matlab’s, since Mathematica data structure is more general and little more complex (ragged arrays, or list of lists is its main data structure) compared to Matlab’s. This is because Mathematica has to also support symbolics in its commands and not just numbers

In Maple the following short cuts can be used enter vectors and matrices: For row vector:  v:=<1|2|3|4> and for column vector  v:=<1,2,3,4> and for a Matrix of say 2 rows and 3 columns  A:=<1,2|3,4|5,6>. Here | acts as a column separator. For transpose of matrix,  A^%T is the short cut notation.

2.1 Multiply matrix with a vector

How to perform the following matrix/vector multiplication?


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



Mathematica





 


Matlab





 

Ada


compile and run


Fortran


compile and run




Maple


\[ \left [{\begin{array}{c} 14\\ 32\\ 50 \end{array}} \right ] \]



 


Python





 

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

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



Mathematica





 


Matlab





 


Fortran





 


Maple


Using <> notation





 


Python

Python uses zero index.





 

2.3 Insert a row into a matrix

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



Mathematica





 


Matlab





 


Maple

Using <<>> notation


Using Matrix/Vector





 


Python





 


Fortran


Compile and run




 

2.4 Insert a column into a matrix

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



Mathematica





 


Matlab





 


Fortran





 


Maple


Using Matrix/Vector





 


Python





 

2.5 Build matrix from other matrices and vectors

  2.5.1 First example
  2.5.2 second example

2.5.1 First example

Given column vectors \(v1= \left [{\begin{array}{c} 1\\ 2\\ 3 \end{array}} \right ]\) and \(v2= \left [{\begin{array}{c} 4\\ 5\\ 6 \end{array}} \right ]\) and \(v3= \left [{\begin{array}{c} 7\\ 8\\ 9 \end{array}} \right ]\) and \(v4=\left [{\begin{array}{c} 10\\ 11\\ 12 \end{array}} \right ]\) generate the matrix \(m=\left [{\begin{array}{cc} v_1&v_2\\ v_3&v_4 \end{array}} \right ] = \left [{\begin{array}{cc} 1&4\\ \noalign{\medskip }2&5 \\ \noalign{\medskip }3&6\\ \noalign{\medskip }7&10\\ \noalign{\medskip } 8&11\\ \noalign{\medskip }9&12\end{array}} \right ] \)

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

The command ArrayFlatten is essential for this in Mathematica.

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



Mathematica


\(\left ({\begin{array}{cc} 1 & 4 \\ 2 & 5 \\ 3 & 6 \\ 7 & 10 \\ 8 & 11 \\ 9 & 12 \\ \end{array}} \right ) \)



 


Matlab





 


Maple


\(\left [{\begin{array}{cc} 1&4\\ \noalign{\medskip }2&5 \\ \noalign{\medskip }3&6\\ \noalign{\medskip }7&10\\ \noalign{\medskip } 8&11\\ \noalign{\medskip }9&12\end{array}} \right ] \)



 


Python


Another way





 


Fortran


Using the RESHAPE command





 

2.5.2 second example

Given mix of matrices and vectors, such as \(v1= \left [{\begin{array}{c} 1\\ 2\\ 3 \end{array}} \right ]\) and \(v2= \left [{\begin{array}{cc} 4& 5\\ 6& 7\\ 8& 9 \end{array}} \right ]\) and \(v3= \left [{\begin{array}{c} 10\\ 11\\ 12 \end{array}} \right ]\) and \(v4=\left [{\begin{array}{c} 13\\ 14\\ 15 \end{array}} \right ]\) and

\(v5=\left [{\begin{array}{c} 16\\ 17\\ 18 \end{array}} \right ]\)

generate the matrix 6 by 3 matrix \(m= \left [{\begin{array}{ccc} v1&v2\\v3&v4&v5\end{array}}\right ]= \left [{\begin{array}{ccc} 1&4&5\\ \noalign{\medskip }2&6&7 \\ \noalign{\medskip }3&8&9\\ \noalign{\medskip }10&13&16 \\ \noalign{\medskip }11&14&17\\ \noalign{\medskip }12&15&18\end{array}} \right ] \)

Mathematica, thanks for function by Kuba at Mathematica stackexachage, this becomes easy to do

Mathematica


Maple


Matlab


Fortran



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



Mathematica









 


Matlab









 


Maple





 

Or


\[ \left [{\begin{array}{cccc} 0.970592781760615697& 0.957506835434297598& 0.0975404049994095246& 0.126986816293506055\\ \noalign{\medskip } 0.157613081677548283& 0.546881519204983846& 0.632359246225409510& 0.905791937075619225\\ \noalign{\medskip } 0.964888535199276531& 0.278498218867048397& 0.913375856139019393& 0.814723686393178936\end{array}} \right ] \]

Fortran



Did not find a build-in support for random numbers from normal distribution, need to look more.

2.7 Generate an n by m zero matrix



Mathematica





 


Matlab





 


Maple


\[ \left [{\begin{array}{cccc} 0&0&0&0\\ 0&0&0&0 \\ 0&0&0&0\end{array}} \right ] \]



 


Fortran





 


Ada




 

2.8 Rotate a matrix by 90 degrees



Mathematica





 


Matlab





 


Maple









 


Fortran

Using additional copy for the matrix





 

2.9 Generate a diagonal matrix with given values on the diagonal

Problem: generate diagonal matrix with \(2,4,6,8\) on the diagonal.



Mathematica





 


Matlab





 


Maple





 

2.10 Sum elements in a matrix along the diagonal



Mathematica





 


Matlab





 


Maple


Another ways


15



 

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



Mathematica


Out[49]= 45



 


Matlab


ans = 45



 


Maple


45



 

2.12 Check if a Matrix is diagonal

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



Mathematica





 


Maple


true



 

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

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



Mathematica





 


Matlab





 


Maple





 

2.14 Replicate a matrix

Given Matrix


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


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



Mathematica





 


Matlab





 

Another way is to use kron() in matlab, and KroneckerProduct in Mathematica and LinearAlgebra[KroneckerProduct] in Maple, which I think is a better way. As follows



Mathematica





 


Matlab





 


Maple





 

2.15 Find the location of the maximum value in a matrix



Mathematica





 


Matlab





 


Maple

This below finds position of first max.





Maple support for such operations seems to be not as strong as Matlab. One way to find locations of all elements is by using explicit loop





 

2.16 Swap 2 columns in a matrix

Give a matrix


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




Mathematica





 


Matlab





 


Maple





 

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



Mathematica

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

















 


Matlab









 


Maple




 

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

Question posted on the net


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



Mathematica





 


Matlab





 


Maple





 

2.19 extract values from matrix given their index

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

Example, given


obtain the entries at \(1,1\) and \(3,3\) which will be \(1\) and \(9\) in this example.



Mathematica





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





 


Matlab





 


Maple





 

2.20 Convert \(N\) by \(M\) matrix to a row of length \(N M\)

Given


covert the matrix to one vector




Mathematica





 


Matlab





 


Maple

Maple reshapes along columns, like Matlab. To get same result as Mathematica, we can transpose the matrix first. To get same result as Matlab, do not transpose.


Notice the result is a row matrix and not a vector. To get a vector


They look the same on the screen, but using whattype we can find the type.




 

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

Example, given Matrix


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




Mathematica









 


Matlab





 


Maple





 

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

Given


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




Mathematica





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





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





 


Matlab





 


Maple





 

2.23 Locate rows in a matrix with column being a string

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


and given list


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




Mathematica





 


Matlab





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





 

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

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

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

Assume valid list of indices.

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

pict

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



Mathematica Three methods are shown.

method 1:

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

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


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





method 2: (due to Mr Wizard at stackoverflow)





method 3: (me)

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





 


Matlab





 

2.25 Convert list of separated numerical numbers to strings

Problem: given a list of numbers such as


convert the above to list of strings




Mathematica





 


Matlab





answer below is due to Bruno Luong at Matlab newsgroup





 

2.26 Obtain elements that are common to two vectors

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



Mathematica





 


Matlab





 

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

Given


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


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



Mathematica





 


Matlab





 

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

Given


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




Mathematica





 


Matlab





 

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

Given


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


In Matlab, the sortrows() command is used. In Mathematica the Sort[] command is now used as is.



Mathematica





 


Matlab





 

2.30 Sort a matrix row-wise using non-first column as key

Given


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


In Matlab, the sortrows() command is used, but now we tell it to use the second column as key.

In Mathematica the SortBy[] command is now used but we tell it to use the second slot as key.



Mathematica





 


Matlab





 

2.31 Replace the first nonzero element in each row in a matrix by some value

Problem: Given a matrix, replace the first nonzero element in each row by some a specific value. This is an example. Given matrix \(A\) below, replace the first non-zero element in each row by \(-1\), then

\(A=\begin{pmatrix} 50 & 75 & 0\\ 50 & 0 & 100\\ 0 & 75 & 100\\ 75 & 100 & 0\\ 0 & 75 & 100\\ 0 & 75 & 100 \end{pmatrix} \) will become \(B=\begin{pmatrix} -1 & 75 & 0\\ -1 & 0 & 100\\ 0 & -1 & 100\\ -1 & 100 & 0\\ 0 & -1 & 100\\ 0 & -1 & 100 \end{pmatrix} \)



Mathematica


Solution due to Bob Hanlon (from Mathematica newsgroup)


Solution by Fred Simons (from Mathematica newsgroup)


Solution due to Adriano Pascoletti (from Mathematica newsgroup)


Solution due to Oliver Ruebenkoenig (from Mathematica newsgroup)


Solution due to Szabolcs Horvát (from Mathematica newsgroup)





 


Matlab


This solution due to Bruno Luong (from matlab newsgroup)


This solution due to Jos from matlab newsgroup





 

2.32 Perform outer product and outer sum between two vector

Problem: Given 2 vectors, perform outer product and outer sum between them. The outer operation takes the first element in one vector and performs this operation on each element in the second vector. This results in first row. This is repeated for each of the elements in the first vector. The operation to perform can be any valid operation on these elements.



Mathematica

using symbolic vectors. Outer product





Outer sum





using numerical vectors. Outer product





Outer sum





 


Matlab

Outer product





Outer sum





 


Maple

Due to Carl Love from the Maple newsgroup









 

2.33 Find the rank and the bases of the Null space for a matrix A

Problem: Find the rank and nullities of the following matrices, and find the bases of the range space and the Null space.

\( A=\begin{pmatrix} 2 & 3 & 3 & 4 \\ 0 & -1 & -2 & 2 \\ 0 & 0 & 0 & 1 \end{pmatrix} \)



Mathematica


3,4















 


Matlab





 

2.34 Find the singular value decomposition (SVD) of a matrix

Problem: Find the SVD for the matrix \(\begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6 \end{bmatrix} \) Notice that in Maple, the singular values matrix, normally called S, is returned as a column vector. So need to call DiagonalMatrix() to format it as expected.



Mathematica









 


Matlab









 


Maple





















 

2.35 Solve \(Ax=b\)

Solve for x in the following system of equations




Mathematica





 


Matlab





 

Fortran


compile and run


2.36 Find all nonzero elements in a matrix

Given a matrix, find the locations and the values of all nonzero elements. Hence given the matrix

\[ \left ({\begin{array}{ccc} 0 & 0 & 1 \\ 10 & 0 & 2 \\ 3 & 0 & 0 \end{array}} \right )\]

the positions returned will be \((1,3),(2,1),(2,3),(3,1)\) and the corresponding values are \(1,10,2,3\).



Mathematica

In Mathematica, standard Mathematica matrix operations can be used, or the matrix can be converted to SparseArray and special named operation can be used on it.









Or standard list operations can be used









 


Matlab









 

2.37 evaluate f(x) on a vector of values

Given a function \(f(x)\) evaluate it for each value contained in a vector. For example, given \(f(x)=x^2\) evaluate it on \((1,2,3)\) such that the result is \((1,4,9)\).



Mathematica





 


Matlab





 

2.38 generates equally spaced N points between \(x_1\) and \(x_2\)



Mathematica


Matlab




2.39 evaluate and plot a \(f(x,y)\) on 2D grid of coordinates

Evaluate \(x \exp ^{-x^2-y^2}\) on 2D cartesian grid between \(x=-2 \cdots 2\) and \(y=-4 \cdots 4\) using \(h=0.4\) for grid spacing.

Mathematica


pict

The above can also be done using Plot3D


pict

I need to sort out the orientation difference between the two plots above.

Matlab


pict

2.40 Find determinant of matrix

Given a square matrix, find its determinant. In Mathematica, the Det[] command is used. In Matlab the det() command is used.



Mathematica





 


Matlab





 

2.41 Generate sparse matrix with \(n\) by \(n\) matrix repeated on its diagonal

Given matrix \(\begin{pmatrix} 1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 9 \end{pmatrix} \), generate the following sparse matrix with this matrix on the diagonal

\[\begin{pmatrix} 1 & 2 & 3 & 0 & 0 & 0 & 0 & 0 & 0\\ 4 & 5 & 6 & 0 & 0 & 0 & 0 & 0 & 0\\ 7 & 8 & 9 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 2 & 3 & 0 & 0 & 0\\ 0 & 0 & 0 & 4 & 5 & 6 & 0 & 0 & 0\\ 0 & 0 & 0 & 7 & 8 & 9 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 2 & 3\\ 0 & 0 & 0 & 0 & 0 & 0 & 4 & 5 & 6\\ 0 & 0 & 0 & 0 & 0 & 0 & 7 & 8 & 9 \end{pmatrix} \]



Mathematica


\[ \left ({\begin{array}{ccccccccc} 1. & 2. & 3. & 0 & 0 & 0 & 0 & 0 & 0 \\ 4. & 5. & 6. & 0 & 0 & 0 & 0 & 0 & 0 \\ 7. & 8. & 9. & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1. & 2. & 3. & 0 & 0 & 0 \\ 0 & 0 & 0 & 4. & 5. & 6. & 0 & 0 & 0 \\ 0 & 0 & 0 & 7. & 8. & 9. & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1. & 2. & 3. \\ 0 & 0 & 0 & 0 & 0 & 0 & 4. & 5. & 6. \\ 0 & 0 & 0 & 0 & 0 & 0 & 7. & 8. & 9. \\ \end{array}} \right ) \]



 


Matlab









 

2.42 Generate sparse matrix for the tridiagonal representation of second difference operator in 1D

The second derivative \(\frac{d^{2}u}{dx^{2}}\) is approximated by \(\frac{u_{i-1}-2u_{i}+u_{i+1}}{h^{2}}\) where \(h\) is the grid spacing. Generate the \(A\) matrix that represent this operator for \(n=4\) where \(n\) is the number of internal grid points on the line



Mathematica


\[ \left ({\begin{array}{cccccccccccc} -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 \\ \end{array}} \right ) \]



 


Matlab





 

2.43 Generate sparse matrix for the Laplacian differential operator \(\nabla ^{2}u\) for 2D grid

\(\nabla ^{2}u=f\) in 2D is defined as \(\frac{\partial ^{2}u}{\partial x^{2}}+\frac{\partial ^{2}u}{\partial y^{2}}=f\) and the Laplacian operator using second order standard differences results in \(\left ( \frac{u_{i-1,j}+u_{i+1,j}+u_{i,j-1}+u_{i,j+1}-4u_{i,j}}{h^{2}}\right ) =f_{i,j}\) where \(h\) is the grid size. The above is solved for \(x\) in the form \(Ax=f\) by generating the \(A\) matrix and taking into consideration the boundary conditions. The follows show how to generate the sparse representation of \(A\). Assuminh the number of unknowns \(n=3\) in one direction. Hence there are 9 unknowns to solve for and the \(A\) matrix will be \(9\) by \(9\).



Matlab





 

2.44 Generate sparse matrix for the Laplacian differential operator \(\nabla ^{2}u\) for 3D grid

The goal is to solve\[ \frac{\partial ^{2}u}{\partial x^{2}}+\frac{\partial ^{2}u}{\partial y^{2}}+\frac{\partial ^{2}u}{\partial z^{2}}=-f(x,y,z) \] On the unit cube. The following diagram was made to help setting up the 3D scheme to approximate the above PDE

pict

The discrete approximation to the Laplacian in 3D is\[ \frac{\partial ^{2}u}{\partial x^{2}}+\frac{\partial ^{2}u}{\partial y^{2}}+\frac{\partial ^{2}u}{\partial z^{2}}=\frac{1}{h^{2}}\left (U_{i-1,j,k}+U_{i+1,j,k}+U_{i,j-1,k}+U_{i,j+1,k}+U_{i,k,k-1}+U_{i,j,k+1}-6U_{i,j,k}\right ) \] For the direct solver, the \(A\) matrix needs to be formulated. From\[ \frac{1}{h^{2}}\left ( U_{i-1,j,k}+U_{i+1,j,k}+U_{i,j-1,k}+U_{i,j+1,k}+U_{i,k,k-1}+U_{i,j,k+1}-6 U_{i,j,k}\right ) =f_{i,j,k}\] Solving for \(U_{i,j,k}\) results in\[ U_{i,j,k}=\frac{1}{6}\left ( U_{i-1,j,k}+U_{i+1,j,k}+U_{i,j-1,k}+U_{i,j+1,k}+U_{i,k,k-1}+U_{i,j,k+1}-h^{2}f_{i,j,k}\right ) \] To help make the \(A\) matrix, a small example with \(n=2,\) is made. The following diagram uses the standard numbering on each node

pict

By traversing the grid, left to right, then inwards into the paper, then upwards, the following \(A\) matrix results

pict

The recursive pattern involved in these \(A\) matrices can now be seen. Each \(A\) matrix contains inside it a block on its diagonal which repeats \(n\) times. Each block in turn contain inside it, on its diagonal, smaller block, which also repeats \(n\) times.

It is easier to see the pattern of building \(A\) by using numbers for the grid points, and label them in the same order as they would be visited, this allowes seeing the connection between each grid point to the other easier. For example, for \(n=2\),

pict

The connections are now more easily seen. Grid point 1 has connection to only points \(2,3,5\). This means when looking at the \(A\) matrix, there will be a \(1\) in the first row, at columns \(2,3,5\). Similarly, point \(2\) has connections only to \(1,4,6\), which means in the second row there will be a \(1\) at columns \(1,4,6\). Extending the number of points to \(n=3\) to better see the pattern of \(A\) results in

pict

From the above it is seen that for example point \(1\) is connected only to \(2,4,10\) and point \(2\) is connected to \(1,3,5,11\) and so on.

The above shows that each point will have a connection to a point which is numbered \(n^{2}\) higher than the grid point itself. \(n^{2}\) is the size of the grid at each surface. Hence, the general \(A\) matrix, for the above example, can now be written as

pict

The recursive structure can be seen. There are \(n=3\) main repeating blocks on the diagonal, and each one of them in turn has \(n=3\) repeating blocks on its own diagonal. Here \(n=3\), the number of grid points along one dimension.

Now that the \(A\) structure is understood, the A matrix can be generated.

Example 1: Using \(n_x=2\), \(n_y=2\), \(n_z=2\). These are the number of grid points in the \(x,y,z\) directions respectively.



Matlab





 

Example 2: Using \(n_x=2\), \(n_y=2\), \(n_z=3\).



Matlab





 

Example 3: Using \(n_x=3\), \(n_y=3\), \(n_z=3\).

Matlab


Matlab result is below.

pict

2.45 Generate the adjugate matrix for square matrix

Given square matrix such as

\[\begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ \end{pmatrix} \]

find the adjugate matrix which is

\[\begin{pmatrix} -3 & 6 & -3 \\ 6 & -12 & 6 \\ -3 & 6 & -3 \\ \end{pmatrix} \]



Mathematica





 


Matlab

Will try to find function in Matlab to do this. But for non-singular matrices only the direct method of finding the inverse and then multiplying by the determinant to recover the adjunct can be used.





The following is due to Matt J from the matlab newsgroup





 

Fortran

Thanks goes to James Van Buskirk and Louisa from comp.lang.fortran for the review and suggestions which lead to improvements of the code.



Ada


compile and run


2.46 Multiply each column by values taken from a row

  2.46.1 Mathematica
  2.46.2 Matlab
  2.46.3 Fortran
  2.46.4 Octave
  2.46.5 Maple

Given a row of the same number of values as the number of columns in a matrix, how to scale each column by the corresponding value in that row? This picture explains more the problem

pict

In Matlab, bsxfun is used.

2.46.1 Mathematica

credit for this solution goes to Bob Hanlon, Adriano Pascoletti, Kurt Tekolste, David Park, and Peter. J. C. Moses from the Math group


   {{2,  6},
    {6,  12},
    {10, 18}}

Another way is to use Inner[] command. Credit for this solution goes to Sswziwa Mukasa and Peter. J. C. Moses from the Math group


Out[66]= {{2,  6},
          {6,  12},
          {10, 18}}

2.46.2 Matlab


ans =
     2     6
     6    12
    10    18

2.46.3 Fortran


#compile and run
>gfortran -std=f2008 t45.f90
>./a.out
    2    6
    6   12
   10   18

2.46.4 Octave

Use automatic broadcasting


2.46.5 Maple


Matrix(3, 2, [[2, 6], [6, 12], [10, 18]])

\[ \left [ \begin{array}{cc} 2&6\\ \noalign{\medskip }6&12 \\ \noalign{\medskip }10&18\end{array} \right ] \]

2.47 extract submatrix from a larger matrix by removing row/column

  2.47.1 Mathematica
  2.47.2 Matlab
  2.47.3 Fortran
  2.47.4 Ada
  2.47.5 Maple

Given 2D matrix \(A\) extract all submatrices by removing row/column.

For example, given the matrix \(\begin{pmatrix} 1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 9 \end{pmatrix} \) then the submatrix for \(\left ( 1,1\right ) \) is\(\begin{pmatrix} 5 & 6\\ 8 & 9 \end{pmatrix} \) obtained by removing the first row and the first column.

In Mathematica, ReplacePart can be used. In Matlab the [] operator can be used. In Fortran, the pack() function can be used.

2.47.1 Mathematica




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



 

2.47.2 Matlab




B =
     5     6
     8     9
B =
     4     6
     7     9
B =
     4     5
     7     8
B =
     2     3
     8     9
B =
     1     3
     7     9
B =
     1     2
     7     8
B =
     2     3
     5     6
B =
     1     3
     4     6
B =
     1     2
     4     5



 

2.47.3 Fortran


compile and run

>gfortran -std=f2008 -fcheck=all -Wall -Wextra -Wconversion-extra t46.f90
>./a.out
   5.0   6.0
   8.0   9.0

   4.0   6.0
   7.0   9.0

   4.0   5.0
   7.0   8.0

   2.0   3.0
   8.0   9.0

   1.0   3.0
   7.0   9.0

   1.0   2.0
   7.0   8.0

   2.0   3.0
   5.0   6.0

   1.0   3.0
   4.0   6.0

   1.0   2.0
   4.0   5.0

2.47.4 Ada


compile and run

>gnatmake foo.adb
gcc-4.6 -c foo.adb
gnatbind -x foo.ali
gnatlink foo.ali
>./foo
 >./foo
  5.00000E+00 6.00000E+00
  8.00000E+00 9.00000E+00

  4.00000E+00 6.00000E+00
  7.00000E+00 9.00000E+00

  4.00000E+00 5.00000E+00
  7.00000E+00 8.00000E+00

  2.00000E+00 3.00000E+00
  8.00000E+00 9.00000E+00

  1.00000E+00 3.00000E+00
  7.00000E+00 9.00000E+00

  1.00000E+00 2.00000E+00
  7.00000E+00 8.00000E+00

  2.00000E+00 3.00000E+00
  5.00000E+00 6.00000E+00

  1.00000E+00 3.00000E+00
  4.00000E+00 6.00000E+00

  1.00000E+00 2.00000E+00
  4.00000E+00 5.00000E+00

2.47.5 Maple


\[ [ \left [ \begin{array}{cc} 5&6\\ \noalign{\medskip }8&9\end{array} \right ] , \left [ \begin{array}{cc} 4&6\\ \noalign{\medskip }7&9 \end{array} \right ] , \left [ \begin{array}{cc} 4&5 \\ \noalign{\medskip }7&8\end{array} \right ] , \left [ \begin{array}{cc} 2&3\\ \noalign{\medskip }8&9\end{array} \right ] , \left [ \begin{array}{cc} 1&3\\ \noalign{\medskip }7&9\end{array} \right ] , \left [ \begin{array}{cc} 1&2\\ \noalign{\medskip }7&8\end{array} \right ] , \left [ \begin{array}{cc} 2&3\\ \noalign{\medskip }5&6 \end{array} \right ] , \left [ \begin{array}{cc} 1&3 \\ \noalign{\medskip }4&6\end{array} \right ] , \left [ \begin{array}{cc} 1&2\\ \noalign{\medskip }4&5\end{array} \right ] ] \]

2.48 delete one row from a matrix

  2.48.1 Mathematica
  2.48.2 Matlab
  2.48.3 Maple

Example, Given the folowing 2D matrix \(A\) delete the second row

\(\begin{pmatrix} 1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 9 \end{pmatrix} \)

2.48.1 Mathematica




Out[81]= {{1,2,3},
          {7,8,9}}

or a little longer solution using Pick


Out[70]= {{1,2,3},
          {7,8,9}}



 

2.48.2 Matlab




A =
     1     2     3
     7     8     9



 

2.48.3 Maple




     [1  2  3]
A := [       ]
     [7  8  9]



 

2.49 delete one column from a matrix

  2.49.1 Mathematica
  2.49.2 Matlab
  2.49.3 Maple

Example, Given the folowing 2D matrix \(A\) delete say the second column

\(\begin{pmatrix} 1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 9 \end{pmatrix}\)

2.49.1 Mathematica




Out[93]= {{1,3},
          {4,6},
          {7,9}}

or a little longer solution using Pick


Out[98]= {{1,3},
          {4,6},
          {7,9}}



 

2.49.2 Matlab




A =
     1     3
     4     6
     7     9



 

2.49.3 Maple




\[ \left [ \begin{array}{cc} 1&3\\ \noalign{\medskip }4&6 \\ \noalign{\medskip }7&9\end{array} \right ] \]



 

2.50 generate random matrix so that each row adds to 1

  2.50.1 Mathematica
  2.50.2 Matlab

Generate the random matrix and divide each row by its total

2.50.1 Mathematica


\[ \left ( \begin{array}{cccc} 0.76862 & 0.917299 & 0.606119 & 0.63735 \\ 0.610077 & 0.208307 & 0.0337861 & 0.473017 \\ 0.388772 & 0.432688 & 0.475881 & 0.68523 \\ \end{array} \right ) \]


{1.15201, 1.99068, 3.05063}


\[ \left ( \begin{array}{cccc} 0.15919 & 0.0553157 & 0.448352 & 0.337142 \\ 0.358594 & 0.264968 & 0.123659 & 0.25278 \\ 0.216269 & 0.248665 & 0.278871 & 0.256195 \\ \end{array} \right ) \]


Out[24]= {1.,1.,1.}

 

2.50.2 Matlab


A =
 0.6787    0.3922    0.7060    0.0462
 0.7577    0.6555    0.0318    0.0971
 0.7431    0.1712    0.2769    0.8235

B =
 0.3723    0.2151    0.3873    0.0253
 0.4913    0.4250    0.0206    0.0630
 0.3689    0.0850    0.1375    0.4087

ans =
    1.0000
    1.0000
    1.0000
 

2.51 generate random matrix so that each column adds to 1

  2.51.1 Mathematica
  2.51.2 Matlab
  2.51.3 Maple

Generate the random matrix and divide each column by its total

2.51.1 Mathematica

This method due to Bob Hanlon from the Math group


Out[31]= {{0.440393,0.945076,0.0527301,0.537288},
          {0.515868,0.565691,0.800959,0.0302484},
          {0.509004,0.143124,0.519455,0.264779}}

Out[32]= {1.46526,1.65389,1.37314,0.832315}


Out[36]= {{0.300555,0.571426,0.038401,0.645535},
          {0.352065,0.342036,0.583303,0.0363425},
          {0.34738,0.0865376,0.378296,0.318123}}

Out[37]= {1.,1.,1.,1.}

Or can use Transpose


{1.,1.,1.,1.}

Another way of doing the above, without the need to transpose 2 times is the following


{1.,1.,1.,1.}
 

2.51.2 Matlab


A =
    0.6948    0.0344    0.7655    0.4898
    0.3171    0.4387    0.7952    0.4456
    0.9502    0.3816    0.1869    0.6463

B =
    0.3541    0.0403    0.4380    0.3097
    0.1616    0.5133    0.4550    0.2817
    0.4843    0.4464    0.1069    0.4086

ans =
    1.0000    1.0000    1.0000    1.0000

 

2.51.3 Maple


                                [ 9  -95   51  24]
                                [                ]
                           A := [99  -20   76  65]
                                [                ]
                                [60  -25  -44  86]

                          b := [168, -140, 83, 175]

                                [3   19   51  24 ]
                                [--  --   --  ---]
                                [56  28   83  175]
                                [                ]
                                [33   1   76   13]
                           B := [--   -   --   --]
                                [56   7   83   35]
                                [                ]
                                [5   5   -44  86 ]
                                [--  --  ---  ---]
                                [14  28  83   175]

                                [1, 1, 1, 1]

2.52 sum all rows in a matrix

  2.52.1 Mathematica
  2.52.2 Matlab

Given the folowing 2D matrix \(A\) find the sum of each row

\[ \begin{pmatrix} 1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 9 \end{pmatrix} \]

2.52.1 Mathematica


Out[2]= {6,
        15,
        24}
 

2.52.2 Matlab


ans =
     6
    15
    24
 

2.53 sum all columns in a matrix

Given the folowing 2D matrix \(A\) find the sum of each column

\(\begin{pmatrix} 1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 9 \end{pmatrix}\)



Mathematica


Matlab




2.54 find in which columns values that are not zero

given matrix \(\begin{pmatrix} 0 & 39 & 0\\ 55 & 100 & 0\\ 34 & 0 & 0\\ 0 & 9 & 0\\ 0 & 0 & 50 \end{pmatrix} \) find the column index which contain values \(>0\) but do it from top row to the bottom row. Hence the result should be \(\left ( \overset{1st\ \operatorname{row}}{\overbrace{2}},\overset{2nd\ \operatorname{row}}{\overbrace{1,2}},\overbrace{1},\overbrace{2},\overbrace{3}\right ) \) this is becuase on the first row, nonzero is in second column, and on the second row, nonzero is at first and third column, and on the third row, nonzero is at the first column, and so on.



Mathematica


Out[18]= {2,1,2,1,2,3}



 


Matlab


2     1     2     1     2     3



 


Maple


[2,1,2,1,2,3]



 

2.55 How to remove values from one vector that exist in another vector

  2.55.1 Mathematica
  2.55.2 Matlab
  2.55.3 Fortran
  2.55.4 Maple

Given vector \(v1=\{1,3,4,5,8,9,20.2,30,-44\}\) remove from it any element that exist in \(v2=\{1,7,4\}\). Notice that Complement[] in Mathematica or setdiff() in Matlab can be used, but they will sort the result. Hence they are not used here in order to keep the order the same.

2.55.1 Mathematica




{3,5,8,9,20.2,30,-44}



 



{3,5,8,9,20.2,30,-44}



 

2.55.2 Matlab




v1 =
3.0000 5.0000 8.0000 9.0000 20.2000 30.0000 -44.0000



 

2.55.3 Fortran


>gfortran f2.f90
>./a.out
3.0000000 5.0000000 8.0000000  9.0000000  20.200001  30.000000 -44.000000

2.55.4 Maple




[3, 5, 8, 9, 20.2, 30, -44]



 

2.56 How to find mean of equal sized segments of a vector

Given vector \(V\) of length \(m\) find the mean of segments \(V(1:n),V(n+1:2n),V(2n+1:3n)...\)In otherwords, equall length segments



Mathematica





 


Matlab





 

2.57 find first value in column larger than some value and cut matrix from there

Given matrix \[ \left ({\begin{array}{ccc} 1 & 5 \\ 2 & 3 \\ 4 & 8 \\ 7 & 2 \end{array}} \right ) \]

search in column 2 of matrix for the first time a value exceeds 6 and return the matrix up to that row. The result should be

\[ \left ({\begin{array}{ccc} 1 & 5 \\ 2 & 3 \end{array}} \right ) \]



Mathematica





 


Matlab





 

2.58 make copies of each value into matrix into a larger matrix

Do the following transformation. We want to take each element of matrix, and replace it by a 2 by 2 matrix with its values in places.

kron() in Matlab and KroneckerProduct in Mathematica can be used for this.

pict



Mathematica


Another method that can be used in this, but I prefer the kron method above:





 


Matlab





 

2.59 repeat each column of matrix number of times

Gives a matrix, repeate each column a number of time, say 3 times, in place to produce a matrix 3 times as wide as the original.

kron() in Matlab and KroneckerProduct in Mathematica can be used for this.



Mathematica


Another method that can be used in this, but the above is better





 


Matlab





 

2.60 How to apply a function to each value in a matrix?

Apply function \(f(x,n)\) to ragged matrix where \(n\) is value taken from the matrix and \(x\) is some other value defined somewhere else. The result of the function should replace the same position in the matrix where \(n\) was.

For example, given matrix \[ mat = \left ({\begin{array}[c]{ccc}1 & 2 & 3\\ 4 & 5 & \\ 7 & 8 & 9 \end{array}}\right ) \]

generate the matrix

\[ \left ({\begin{array} [c]{ccc}f[x,1] & f[x,2] & f[x,3]\\ f[x,4] & f[x,5] & \\ f[x,7] & f[x,8] & f[x,9] \end{array}}\right ) \]



Mathematica

The trick is to use Map with 2 at end.


Matlab

to do



 

2.61 How to sum all numbers in a list (vector)?

Given a vector or list \(1,2,3,4,5,6,7,8,9,10\) how to find the sum of all its elements?



Mathematica






55



 


Matlab


55



 


Maple






55



 

2.62 How to find maximum of each row of a matrix?

Given \[ \left ({\begin{array}{cc} 1 & 2 \\ 3 & 4 \end{array}} \right ) \] Find the maximum of each row. The result should be \((2,4)\). (to find the min, just change Max with Min below.



Mathematica





 


Matlab





 

2.63 How to find maximum of each column of a matrix?

Given \[ \left ({\begin{array}{cc} 1 & 2 \\ 3 & 4 \end{array}} \right ) \] Find the maximum of each column. The result should be \((3,4)\). (to find the min, just change Max with Min below.



Mathematica





 


Matlab





 

2.64 How to add the mean of each column of a matrix from each column?

Given \[ \left ({\begin{array}{cc} 1 & 2 \\ 3 & 4 \end{array}} \right ) \] Find the mean of each column, and add this mean to each element of the corresponding column. The result should be

\[ \left ({\begin{array}{cc} 3 & 5 \\ 5 & 7 \end{array}} \right ) \]

To subtract the mean, just replace Plus with Subtract below for Mathematica and replace @plus with @minus for Matlab. This shows that Matlab bsxfun is analogue to Mathematica’s MapThread.



Mathematica





 


Matlab





 

2.65 How to add the mean of each row of a matrix from each row?

Given \[ \left ({\begin{array}{cc} 1 & 2 \\ 3 & 4 \end{array}} \right ) \] Find the mean of each row, and add this mean to each element of the corresponding row. The result should be

\[ \left ({\begin{array}{cc} 2.5 & 3.5 \\ 6.5 & 7.5 \end{array}} \right ) \]

To subtract the mean, just replace Plus with Subtract below for Mathematica and replace @plus with @minus for Matlab. This shows that Matlab bsxfun is analogue to Mathematica’s MapThread.. The main difference is that Matlab is a little bit more consistent in this example, since in Matlab all operations are done column wise. Hence mean(A) takes the mean of each column (same as Matematica in this one case), but bsxfun also acts column wise on the matrix A, while Mathematica Map and MapThread act row wise (list of lists). One just needs to be careful about this order difference.



Mathematica





 


Matlab





 

2.66 Find the different norms of a vector

Problem: Given the vector say \[ v={1,2,4} \]

Find its norm for \(p=1,2,\infty \)



Mathematica





 


Matlab





 

2.67 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


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









 


Matlab










1



 

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

pict



Mathematica

















 


Matlab









 

2.69 Linear convolution of 2 sequences

Problem: Given 2 sequences \(x_{1}[n]\) and \(x_{2}[m]\), determine the linear convolution of the 2 sequences. Assume \(x_{1}=[1,2,3,4,5]\) and \(x_{2}=[3,4,5]\).



Mathematica





 


Matlab





 

2.70 Circular convolution of two sequences

  2.70.1 example 1. The sequences are of equal length
  2.70.2 example 2. The sequences are of unequal length

Problem: Given 2 sequences \(x_{1}[n]\) and \(x_{2}[m]\), determine the circular convolution of the 2 sequences.

2.70.1 example 1. The sequences are of equal length



Mathematica





 


Matlab





 

2.70.2 example 2. The sequences are of unequal length



Mathematica





 


Matlab





 

2.71 Linear convolution of 2 sequences with origin at arbitrary position

For simple case where the 2 sequences are assumed to start at n=0, then this can be done in the time domain using the ListConvolve in Mathematica and using conv in Matlab.

The harder case is when the origin is not located at the first element in the sequence. I.e. one or both of the sequences is not causal.

Mathematica

case 1

Convolve 2 sequences where the origin is assumed to be at the first element in each sequence.



case 2

Convolve 2 sequences where the origin is located at different location from the first element of the sequence, use DiracDelta function, and the DTFT approach.

Example convolve x=1,2,0,2,1,4,01,2,2 with h=1/4,1/2,1/4 where the origin of h is under the second element 1/2.


\[ \frac{1}{4} e^{-i w} \left (1+e^{i w}\right )^2 \]

Write down the x sequence and take its DTFT


\[ e^{-7 i w} \left (e^{i w}-e^{2 i w}+4 e^{3 i w}+e^{4 i w}-2 e^{5 i w}+2 e^{6 i w}+e^{7 i w}+2\right ) \]

Now multiply the DTFT’s and take the inverse


pict

Now convolve z with h again, where z is the convolution of x and h found above. This can be done as follows in one command


pict

2.72 Visualize a 2D matrix

Problem: Given a 2 dimensional matrix, say \(m\times n\,\), how to visualize its content?

These are some examples showing how to visualize a matrix as a 3D data, where the height is taken as the values of the matrix entries, and the \(x,y\,\ \)indices as the coordinates.

Mathematica


pict


pict

Matlab


pict


pict

Maple


pict


pict

2.73 Find the cross correlation between two sequences

Problem: Given \begin{align*} 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] \end{align*}

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


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


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



 

Maple


\[ [ 7.0, 0.0, 33.0,- 14.0, 36.0, 19.0,- 9.0, 10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] \]

Not able to find out now why Maple result is different. May be definition used is different, no time now to find out.

2.74 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 \[ A=\begin{bmatrix} 0 & 1 & 1 & 2\\ 1 & 2 & 3 & 4\\ 2 & 0 & 2 & 0 \end{bmatrix} \]

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


2











 


Matlab









 

2.75 Solve \(A x= b\) and display the solution

Problem: Solve for x given that \(A x=b\) where \[ A=\left ({\begin{array} [c]{cc}1 & 2\\ 1 & 3 \end{array}} \right ) \] and \[ b=\left ({\begin{array} [c]{c}5\\ 8 \end{array}} \right ) \]

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


pict







 


Matlab





 

2.76 Determine if a set of linear equations \(A x= 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 the system of equations be \begin{align*} y & =x-1\\ y & =2x+1 \end{align*}

So \[ A=\left ({\begin{array} [c]{cc}1 & 1\\ -2 & 1 \end{array}} \right ) \] and \[ b=\left ({\begin{array} [c]{c}1\\ -1 \end{array}} \right ) \]

Mathematica


pict


2


2

The above algorithm can now be run as follows


The output of the above is


Matlab


Output is


2.77 Given a set of linear equations automatically generate the matrix \(A\) and vector \(b\) and solve \(A x=b\)

Problem: Given \[{\begin{array} [c]{ccc}4x+4y+2z & = & 12\\ 5x+6y+3z & = & 7\\ 9x+7y+10z & = & 9 \end{array}} \]

Automatically convert it to the form \(Ax=b\) and solve



Mathematica


\[ \left ({\begin{array}{ccc} 4 & 3 & 2 \\ 5 & 6 & 3 \\ 9 & 7 & 10 \\ \end{array}} \right ) \]




\[ \{-12,-7,-9\} \]







 

2.78 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). One reason to convert Matrix \(A\) to its row echelon form, is to find the rank of \(A\). If matrix \(A\) is a \(4\times 4\), and when converted to its row echelon form we find that one of the rows is all zeros, then the rank of \(A\) will be 3 and not full rank.



Mathematica


\[ \left ({\begin{array}{cccc} 1 & 1 & -2 & 1 \\ 3 & 2 & 4 & -4 \\ 4 & 3 & 3 & -4 \\ \end{array}} \right ) \]




\[ \left ({\begin{array}{cccc} 1 & 0 & 0 & 2 \\ 0 & 1 & 0 & -3 \\ 0 & 0 & 1 & -1 \\ \end{array}} \right ) \]



 


Matlab





 


Maple


\[ \left [{\begin{array}{cccc} 1&0&0&2\\ \noalign{\medskip }0&1&0&-3 \\ \noalign{\medskip }0&0&1&-1\end{array}} \right ] \]



 

2.79 Convert 2D matrix to show the location and values

Given \begin{align*} A =& \left ({\begin{array}{cccc} 41 & 45 & 49 & 53 \\ 42 & 46 & 50 & 54 \\ 43 & 47 & 51 & 55 \\ 44 & 48 & 52 & 56 \\ \end{array}} \right ) \end{align*}

Generate the matrix

\[ \left ({\begin{array}{ccc} 1 & 1 & 41 \\ 2 & 1 & 42 \\ 3 & 1 & 43 \\ 4 & 1 & 44 \\ 1 & 2 & 45 \\ 2 & 2 & 46 \\ 3 & 2 & 47 \\ 4 & 2 & 48 \\ 1 & 3 & 49 \\ 2 & 3 & 50 \\ 3 & 3 & 51 \\ 4 & 3 & 52 \\ 1 & 4 & 53 \\ 2 & 4 & 54 \\ 3 & 4 & 55 \\ 4 & 4 & 56 \\ \end{array}} \right ) \]

Which gives at each row, the location and the value in the original matrix.

Mathematica


Another way


But I think the simplist is to use Table


Matlab


gives


Maple


\[ A= \left [ \begin{array}{cccc} 41&42&43&44\\ \noalign{\medskip }45&46&47&48\\ \noalign{\medskip }49&50&51&52\\ \noalign{\medskip }53&54&55&56 \end{array} \right ] \]


\[ \left [ \begin{array}{ccc} 1&1&41\\ \noalign{\medskip }1&2&45 \\ \noalign{\medskip }1&3&49\\ \noalign{\medskip }1&4&53 \\ \noalign{\medskip }2&1&42\\ \noalign{\medskip }2&2&46 \\ \noalign{\medskip }2&3&50\\ \noalign{\medskip }2&4&54 \\ \noalign{\medskip }3&1&43\\ \noalign{\medskip }3&2&47 \\ \noalign{\medskip }3&3&51\\ \noalign{\medskip }3&4&55 \\ \noalign{\medskip }4&1&44\\ \noalign{\medskip }4&2&48 \\ \noalign{\medskip }4&3&52\\ \noalign{\medskip }4&4&56\end{array} \right ] \]

2.80 Find rows in matrix with zeros in them, and then remove the zeros

Given \(A = \{\{1, 0, 9\}, \{5, 0, 6\}, \{4, 1, 9\}, \{7, 0, 11\}, \{8, 1, 2\}\}\)

Find rows with zero in middle, and then remove the zeros from these found.

The result should be \(\{\{1, 9\}, \{5, 6\}, \{7, 11\}\}\)

Mathematica

Many ways to do this. See this


{{1, 9}, {5, 6}, {7, 11}}

Maple


{{1, 9}, {5, 6}, {7, 11}}

2.81 How to apply a function to two lists at the same time?

Given list \(a = \{1,2,3,4\}\) and list \(b=\{5,6,7,8\}\) how to to call function \(f(x,y)\) by taking \(x,y\) from \(a,b\) one a time so that the result gives \(f(1,5),f(2,6),f(3,7),f(4,8)\)?



Mathematica


{f[1, 5], f[2, 6], f[3, 7], f[4, 8]}

Maple


[f(1, 5), f(2, 6), f(3, 7), f(4, 8)]



 

2.82 How to apply a function to two lists are the same time, but with change to entries?

  2.82.1 example 1
  2.82.2 example 2

2.82.1 example 1

Given list \(a = \{1,2,3,4\}\) and list \(b=\{5,6,7,8\}\) how to to call function \(f(x,y)\) by taking \(\sin (x),\cos (y)\) from \(a,b\) one a time so that the result gives \[ f(\sin (1),\cos (5)),f(\sin (2),\cos (6)),f(\sin (3),\cos (7)),f(\sin (4),\cos (8)) \]

Mathematica


{f[Sin[1], Cos[5]], f[Sin[2], Cos[6]], f[Sin[3], Cos[7]], f[Sin[4], Cos[8]]}

Maple


[f(sin(1), cos(5)), f(sin(2), cos(6)), f(sin(3), cos(7)), f(sin(4), cos(8))]

2.82.2 example 2

Given list \(a = \{1,2,3,4\}\) and list \(b=\{5,6,7,8\}\) how to to call function \(f(2+x+x^2+\sin (x),2+\cos (y))\) by taking \(x,y\) from \(a,b\) one a time so that the result gives \[ \{f(3+\sin (1),5+\cos (5)),f(8+\sin (2),6+\cos (6)),f(15+\sin (3),7+\cos (7)),f(24+\sin (4),8+\cos (8))\} \]

Mathematica


{f[3+Sin[1],5+Cos[5]],f[8+Sin[2],6+Cos[6]],f[15+Sin[3],7+Cos[7]],f[24+Sin[4],8+Cos[8]]}

Maple

In Maple, since it does not have slot numbers # to use, it is better to make a function on the fly, which does the exact same thing.


[f(3 + sin(1), 2 + cos(5)), f(8 + sin(2), 2 + cos(6)), f(15 + sin(3), 2 + cos(7)), f(24 + sin(4), 2 + cos(8))]

2.83 How to select all primes numbers from a list?

Given list of numbers from 1 to 100, select the prime numbers.

Mathematica


{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97}

Maple


[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

2.84 How to collect result inside a loop when number of interation is not known in advance?

Lets say we want to collect result obtained inside loop, but do not know in advance how many iteration we need.

In Mathematica, Sow and Reap are used. In Maple, an Array can be used or a queue or a table.

Mathematica


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

Maple


[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

2.85 How flip an array around?

Given \(A=[1,2,3,4,5]\) change it to \(A=[5,4,3,2,1]\)



Mathematica


{5, 4, 3, 2, 1}



Matlab


ans =
 5     4     3     2     1



Maple


[5, 4, 3, 2, 1]



 

2.86 How to divide each element by its position in a list?

  2.86.1 Mathematica
  2.86.2 Maple
  2.86.3 Matlab

Given \(A=[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]\) change it to \[ [2,{\frac{3}{2}},{\frac{5}{3}},{\frac{7}{4}},{\frac{11}{5}},{\frac{13}{6}},{\frac{17}{7}},{\frac{19}{8}},{\frac{23}{9}},{\frac{29}{10}}] \]

2.86.1 Mathematica


\[ \left \{2,\frac{3}{2},\frac{5}{3},\frac{7}{4},\frac{11}{5},\frac{13}{6},\frac{17}{7},\frac{19}{8},\frac{23}{9},\frac{29}{10}\right \} \]

2.86.2 Maple


\[ [2,{\frac{3}{2}},{\frac{5}{3}},{\frac{7}{4}},{\frac{11}{5}},{\frac{13}{6}},{\frac{17}{7}},{\frac{19}{8}},{\frac{23}{9}},{\frac{29}{10}}] \]

2.86.3 Matlab


2.0000    1.5000    1.6667    1.7500    2.2000    2.1667    2.4286    2.3750    2.5556    2.9000

Chapter 3
signal processing, Image processing, graphics, random numbers

 3.1 Generate and plot one pulse signal of different width and amplitude
 3.2 Generate and plot train of pulses of different width and amplitude
 3.3 Find the discrete Fourier transform of a real sequence of numbers
 3.4 Generate uniform distributed random numbers
 3.5 Obtain Fourier Series coefficients for a periodic function
 3.6 Determine and plot the CTFT (continuous time Fourier Transform) for continuous time function
 3.7 Determine the DTFT (Discrete time Fourier Transform) for discrete time function
 3.8 Determine the Inverse DTFT (Discrete time Fourier Transform)
 3.9 Use FFT to display the power spectrum of the content of an audio wav file
 3.10 Plot the power spectral density of a signal
 3.11 Display spectrum of 2D image
 3.12 Obtain the statistical maximum likelihood estimates (MLE) of probability distributions
 3.13 Make a histogram of data sampled from some probability distribution
 3.14 apply a filter on 1D numerical data (a vector)
 3.15 apply an averaging Laplacian filter on 2D numerical data (a matrix)
 3.16 How to find cummulative sum
 3.17 How to find the moving average of a 1D sequence?
 3.18 How to select N random values from a set of numbers?
 3.19 How to sample a sin signal and plot it?
 3.20 How find the impulse response of a difference equation?

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

Problem: Generate one signal of some width and height.



Mathematica


pict



 

Matlab


pict

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

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



Mathematica


pict



 


Matlab


pict



 

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

Given the sequence of numbers \(x(n)=\left [{1,2,3}\right ] \), Find \(X(k) ={\displaystyle \sum \limits _{m=0}^{N-1}} x(m) e^{-j\frac{2\pi }{N}mk}\) where \(N\) is the length of the data sequence \(x(m)\) and \(k=0 \cdots N-1\)



Mathematica





 


Matlab





 


Maple

Maple need an Array for input, not list, so have to convert this is little strange





 

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.


Fortran

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


3.4 Generate uniform distributed random numbers

  3.4.1 How to generate 5 uniform distributed random numbers from 0 to 1?
  3.4.2 How to generate 5 random numbers from a to b?
  3.4.3 How to generate MATRIX of random numbers from a to b?

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



Mathematica





 


Matlab





 

Fortran


compile and run


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





 


Matlab





 

Fortran


compile and run


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



Matlab



Fortran


compile and run


3.5 Obtain Fourier Series coefficients for a periodic function

Given a continuous time function \(x\left ( t\right ) =\sin (\frac{2 \pi }{T_0} t) \), find its Fourier coefficients \(c_{n}\), where \[ x\left ( t\right ) =\sum _{n=-\infty }^{\infty }c_{n}e^{j\omega _{0}nt}\] and \[ c_{n}=\frac{1}{T_{0}}\int _{-\frac{T_{0}}{2}}^{\frac{T_{0}}{2}}x(t)e^{-j\omega _{0}nt}dt \] Where \(T_{0}\) is the fundamental period of \(x\left ( t\right ) \) and \(\omega _{0}=\frac{2\pi }{T_{0}}\)

Mathematica


\[ \frac{1}{2} i e^{-\frac{2 i \pi t}{\text{T0}}}-\frac{1}{2} i e^{\frac{2 i \pi t}{\text{T0}}} \]


pict


pict


pict


pict


\[ \frac{e^{-i t}}{2}+\frac{e^{i t}}{2}-\frac{1}{4} e^{-2 i t}-\frac{1}{4} e^{2 i t}+\frac{1}{2} \]


pict


pict

Plot phase


pict

Find Fourier series for periodic square wave


pict


\[ 0.5 +0.31831 e^{-i \pi t} + 0.31831 e^{i \pi t} - 0.106103 e^{-3 i \pi t} - 0.106103 e^{3 i \pi t} + 0.063662 e^{-5 i \pi t}+0.063662 e^{5 i \pi t} \]


pict


pict

Show phase


pict

3.6 Determine and plot the CTFT (continuous time Fourier Transform) for continuous time function

Plot the CTFT \(X(\omega )\) of \(x(t)=3 \sin (t)\). By definition \begin{align*} X(\omega ) &= \int _{t=-\infty }^{\infty } x(t) e^{-i\omega t} \mathop{dt} \end{align*}

Mathematica




pict

Matlab



Need to do the plot.

Maple


gives


3.7 Determine the DTFT (Discrete time Fourier Transform) for discrete time function

Find DTFT \(X(\Omega )\) of \(x[n]= \sin (\frac{\pi n}{8})\). By definition \begin{align*} X(\Omega ) &= \sum _{n=-\infty }^{\infty } x[n] e^{-i\Omega n} \end{align*}

Mathematica


Which gives


3.8 Determine the Inverse DTFT (Discrete time Fourier Transform)

Find \(x[n]\), give its DTFT \(X(\Omega )\) \begin{align*} x[n] &= \frac{1}{2 \pi } \int _{-\pi }^{\pi } X(\Omega ) e^{i\Omega n} \mathop{d\Omega } \end{align*}

Mathematica


Which gives


3.9 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









now load the samples and do power spectrum










pict



 

Matlab


pict

3.10 Plot the power spectral density of a signal

Problem: Given a signal \[ 3\sin \left ( 2\pi f_{1}t\right ) +4\cos \left ( 2\pi f_{2}t\right ) +5\cos \left ( 2\pi f_{3}t\right ) \] for the following frequency values (in Hz) \[ f_{1}=20 ,f_{2}=30, f_{3}=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


pict


pict

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

pict

And using InterpolationOrder -> 2

pict

Matlab


pict


pict

3.11 Display spectrum of 2D image

Mathematica


pict








pict


pict


pict

Matlab


pict


pict


pict

3.12 Obtain the statistical maximum likelihood estimates (MLE) of probability distributions

This example uses the normal distribution and Poisson as examples. The maximum likelihood estimates of the population parameters is found by solving equation(s) using the standard method of taking logs and differentiation to solve for the maximum. Mathematica and Matlab version will be added at later time.

Maple


pict

3.13 Make a histogram of data sampled from some probability distribution

This example uses the normal distribution. First random data is generated, then histogram of the data is made.

Matlab do not have an option, (unless I missed it) to make a relative histogram (this is where the total area is 1) but not hard to implement this.



Mathematica


pict



 


Matlab


pict



 


Maple


pict



 

3.14 apply a filter on 1D numerical data (a vector)

Problem: suppose we want to apply say an averaging filter on a list of numbers. Say we want to replace each value in the list by the average 3 values around each value (a smoothing filter).

In Mathematica, ListConvolve is used, in Matlab conv() is used.



Mathematica





 


Matlab





 

3.15 apply an averaging Laplacian filter on 2D numerical data (a matrix)

Problem: Apply a Laplacian filter on 2D data. In Mathematica, ListConvolve is used, in Matlab conv2() is used.



Mathematica





 


Matlab





 

3.16 How to find cummulative sum

compute \(\sum \limits _{k=1}^{10} \frac{1}{k\left ( k+1\right ) }\)



Mathematica





 


Matlab





 

3.17 How to find the moving average of a 1D sequence?

Given some sequence such as \(1,2,3,4,5,6,7\) how to find the moving average for different window sizes?



Mathematica

For window size \(k=2\)





For window size \(k=3\)





 


Matlab

For a window size \(k=2\)





For window size \(k=3\)