up

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

Nasser M. Abbasi

April 9, 2014

Contents

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

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

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

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

One day I might make a PDF file of this whole web page.

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

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

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

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

PIC

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

PIC

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

PIC

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

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

H s = 1 s2 +0.2s+1

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

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

Mathematica
Clear["Global`*"];  
SetDirectory[NotebookDirectory[]]  
sys = TransferFunctionModel[1/(s^2 + 2/10 s + 1), s];  
 
yStep = Assuming[t > 0,  
          Simplify@OutputResponse[sys, UnitStep[t], t]];  
yImpulse = Simplify@OutputResponse[sys, DiracDelta[t], t];  
p = Grid[{  
   {Plot[Evaluate@{yStep, yImpulse}, {t, 0, 50},  
     PlotRange -> {{0, 50}, {-0.8, 2.0}},  
     Frame -> True,  
     FrameLabel -> {{"y(t)", None},  
                   {t, "step and impulse reponse"}},  
     GridLines -> Automatic,  
     GridLinesStyle -> Dashed,  
     ImageSize -> {300, 300},  
     PlotStyle -> {Red, Black},  
     AspectRatio -> 1]  
    },  
   {  
    Row[{"step response=", Re@yStep}]  
    },  
   {Row[{"impulse response=", Re@yImpulse}]}  
   }, Alignment -> Left, Frame -> True]  
 
yStep  
yImpulse  
 
{1 - E^(-t/10) Cos[(3 Sqrt[11] t)/10]  
  - (  
  E^(-t/10) Sin[(3 Sqrt[11] t)/10])/(3 Sqrt[11])}  
 
{(10 E^(-t/10) HeavisideTheta[t]  
  Sin[(3 Sqrt[11] t)/10])/(3 Sqrt[11])}

PIC

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

PIC

Maple

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

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

PIC

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

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

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

when the input is given by

u(t) = sin(t)

and display the response and input on the same plot.

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

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

PIC

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

PIC

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

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

PIC

4 Obtain the poles and zeros of a transfer function

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

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

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

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

5 Obtain the continuous time transfer function given the poles, zeros and gain

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

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

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

6 Convert transfer function to state space representation

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

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

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

PIC

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

PIC

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

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

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

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

The plot itself is not shown in these examples.

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

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

8 Generate Bode plot of a transfer function

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

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

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

PIC

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

PIC

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

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

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

Use sampling period of 0.1 seconds.

Mathematica
Clear["Global`*"];  
tf = TransferFunctionModel[(s + 1)/(s^2 + 1 s + 1), s];  
ds = ToDiscreteTimeModel[tf, 0.1, z,  
                         Method -> "ZeroOrderHold"]

PIC

{gm, pm} = GainPhaseMargins[ds];  
gm  
   Out[28]={{31.4159,19.9833}}  
pm  
   Out[29]={{1.41451,1.83932},{0.,Pi}}  
 
max = 100;  
yn  = OutputResponse[ds, Table[UnitStep[n], {n, 0, max}]];  
 
ListPlot[yn,  
         Joined -> True,  
         InterpolationOrder -> 0,  
         PlotRange -> {{0, max}, {0, 1.4}},  
         Frame -> True,  
         FrameLabel -> {{"y[n]", None},  
                        {"n", "unit step response"}},  
         ImageSize -> 400,  
         AspectRatio -> 1]  

PIC

Matlab
clear all; close all;  
Ts   = 0.1; %sec, sample period  
s    = tf('s');  
sys  = ( s + 1 ) / (  s^2  + s  +  1);  
sysd = c2d(sys,Ts,'zoh');  %convert from s to z domain  
tf(sysd)  
step(sysd,10);  
[Gm,Pm,Wcg,Wcp] = margin(sysd)  
set(gcf,'Position',[10,10,410,410]);  
-------  
 
   0.09984 z - 0.09033  
  ----------------------  
  z^2 - 1.895 z + 0.9048  
 
Sample time: 0.1 seconds  
Discrete-time transfer function.  
 
Gm =  
   19.9833  
 
Pm =  
  105.3851  
 
Wcg =  
   31.4159  
 
Wcp =  
    1.4145

PIC

10 Convert transfer function to state space

Problem: Given the continuous time S transfer function defined by

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

convert to state space and back to transfer function.

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

PIC

TransferFunctionModel[ss]

PIC

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

11 Obtain partial-fraction expansion

Problem: Given the continuous time S transfer function defined by

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

obtain the partial-fractions decomposition.

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

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

Matlab
clear all;  
s=tf('s');  
tf_sys    = ( s^4+8*s^3+16*s^2+9*s+6)  /  (s^3+6*s^2+11*s+6);  
[num,den] = tfdata(tf_sys,'v');  
%this gives  r1/(s-p1) + r2/(s-p2) +..... + k  
[r,p,k]   = residue(num,den)  
------  
r =  
 
   -6.0000  
   -4.0000  
    3.0000  
 
p =  
 
   -3.0000  
   -2.0000  
   -1.0000  
 
k =  
     1     2

12 Obtain Laplace transform for a piecewise functions

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

PIC

Function f(t) to obtain its Laplace transform

Comment: Mathematica solution was easier than Matlab's.

In Matlab the definition of the Laplace transform is applied to each piece separately and the result added.

Not finding the piecewise maple function to access from inside MATLAB did not help.

Mathematica
Remove["Global`*"];  
f[t_] := Piecewise[{{0,t<0},  
                    {t,t>= 0 && t< T},  
                    {T, t>T}}]  
Simplify[LaplaceTransform[f[t],t,s] ,{T>0}]  
 
Out[]= (1 - E^((-s)*T))/s^2

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

13 Obtain Inverse Laplace transform of a transfer function

Problem: Obtain the inverse Laplace transform for the function

H(s) = s4 +5s3 +6s2 +9s+30 s4 +6s3 +21s2 +46s+30

Mathematica
Remove["Global`*"];  
f = (s^4+5 s^3+6 s^2+9 s+30)/(s^4+6 s^3+21 s^2+46 s+30);  
InverseLaplaceTransform[f,s,t];  
Expand[FullSimplify[%]]

PIC

Matlab
clear all;  
syms s t  
f = (s^4+5*s^3+6*s^2+9*s+30)/(s^4+6*s^3+21*s^2+46*s+30);  
pretty(f)  
 
    4      3      2  
   s  + 5 s  + 6 s  + 9 s + 30  
  -----------------------------  
   4      3       2  
  s  + 6 s  + 21 s  + 46 s + 30  
 
pretty(ilaplace(f))  
                                                    /            399 sin(3 t) \  
                                        253 exp(-t) | cos(3 t) + ------------ |  
  23 exp(-t)   3 exp(-3 t)                          \                253      /  
  ---------- - ----------- + dirac(t) - ---------------------------------------  
      18           26                                     117  

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

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

H s = ωn2 s2 +2ξωns+ωn2

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

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

PIC

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

PIC

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

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

H s = ωn2 s2 +2ξωns+ωn2

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

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

PIC

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

PIC

16 Convert differential equation to transfer functions and state space

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

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

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

PIC

tf = Simplify[TransferFunctionModel[ss]]

PIC

Matlab
clear all;  
 
syms y(t) Y  
 
eq  = 3*diff(y,2)+2*diff(y)+y;  
lap = laplace(eq);  
lap = subs(lap,'y(0)',0);  
lap = subs(lap,'D(y)(0)',0);  
lap = subs(lap,'laplace(y(t),t,s)',Y);  
H   = solve(lap-1,Y);  
pretty(H)  
 
[num,den] = numden(H);  
num = sym2poly(num);  
den = sym2poly(den);  
 
[A,B,C,D] = tf2ss(num,den)  
 
 
        1  
  --------------  
     2  
  3 s  + 2 s + 1  
 
A =  
   -0.6667   -0.3333  
    1.0000         0  
 
B =  
     1  
     0  
 
C =  
         0    0.3333  
 
D =  
     0

17 View the steady state error for 2nd order LTI system as the un damped natural frequency changes

Problem: Given the transfer function

H s = ωn2 s2 +2ξωns+ωn2

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

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

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

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

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

Mathematica

ramp input

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

PIC

step input

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

PIC

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

Problem: Convert

H s = 1 s2 +10s+10

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

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

PIC

Matlab
s    = tf('s');  
sys  = 1/(s^2 + 10*s + 20);  
T    = 0.01; %seconds  
sysz = c2d(sys,T,'zoh')  
 
------------->  
 
sysz =  
 
  4.837e-05 z + 4.678e-05  
  -----------------------  
  z^2 - 1.903 z + 0.9048  
 
Sample time: 0.01 seconds  
Discrete-time transfer function.

19 Show the use of the inverse Z transform

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

19.1 problems/example 1

Problem: Given

F z = z z1

find x[n] = F1 z which is the inverse Ztransform.

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

PIC

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

PIC

19.2 problems/example 2

Problem: Given

F z = 5z z12

find x[n] = F1 z

Mathematica

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

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

PIC

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

PIC

20 Find the Z transform of sequence x(n)

20.1 problems/example 1

Find the Z transform for the unit step discrete function

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

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

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

20.2 problems/example 2

Find the Z transform for x[n] = 13nu n+ 0.9n3u n

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

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

21 Sample a continuous time system

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

PIC

system diagram

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

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

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

PIC

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

PIC

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

PIC

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

PIC

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

Problem: Given

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

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

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

PIC

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

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

PIC

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

23 Do linear convolution of 2 sequences

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

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

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

24 Do circular convolution of two sequences

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

24.1 problems/example 1. The sequences are of equal length

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

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

24.2 problems/example 2. The sequences are of unequal length

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

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

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

Mathematica
Remove["Global`*"];  
m={{3,-1,1,1,0,0},  
   {1, 1,-1,-1,0,0},  
  {0,0,2,0,1,1},  
  {0,0,0,2,-1,-1},  
  {0,0,0,0,1,1},  
  {0,0,0,0,1,1}};  
MatrixForm[m]

PIC

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

PIC

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

26 Plot the power spectral density of a signal

Problem: Given a signal

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

for the following frequency values (in Hz)

f1 = 20,f2 = 30,f3 = 50

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

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

PIC

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

PIC

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

PIC

And using InterpolationOrder -> 2

PIC

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

PIC

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

PIC

27 Solve the continuous-time algebraic Riccati equation

Problem: Solve for X in the Riccati equation

AX +XAXBR1BX +CC = 0

given

A = 3 2 1 1 B = 0 1 C = 1 1 R = 3

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

PIC

Matlab

MATLAB solution requires the control system toolbox to be installed.

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

28 Solve the discrete-time algebraic Riccati equation

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

1 s(s+0.5)

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

The Riccati equation is given by

AX +XAXBR1BX +CC = 0

Let R = [3].

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

PIC

ss = StateSpaceModel[dsys]

PIC

a = ss[[1,1]];  
b = ss[[1,2]];  
c = ss[[1,3]];  
d = ss[[1,4]];  
r = {{3}};  
DiscreteRiccatiSolve[{a,b},{Transpose[c].c,r}];  
MatrixForm[%]

PIC

Matlab

MATLAB solution requires the control system toolbox.

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

29 Display the impulse response of discrete time system H(z) along with the impulse response of its continuous time approximation H(s)

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

Mathematica
maxSimulationTime = 10;  
samplePeriod      = 0.5;  
tf = TransferFunctionModel[z/(z^2-1.4 z+0.5),  
          z,  
          SamplingPeriod->samplePeriod];  
 
(*Find its impulse response*)  
discreteResponse=First@OutputResponse[tf,  
                            DiscreteDelta[k],  
                            {k,0,maxSimulationTime}]  
 
Out[4]= {0.,1.,1.4,1.46,1.344,1.1516,0.94024,0.740536,  
         0.56663,0.423015,0.308905,0.22096,0.154891,  
         0.106368,0.0714694,0.0468732,0.0298878,  
         0.0184063,0.0108249,0.00595175,0.00291999}  
 
(*approximate to continuous time, use ZeroOrderHold*)  
ctf = ToContinuousTimeModel[tf, s, Method -> "ZeroOrderHold"]  
 
Out[]= TransferFunctionModel[{  
         {{5.6099211033954965 + 1.255589621432896*s}},  
          0.5609921103395501 + 1.3862943611198901*s + s^2},  
          s,  
          SamplingPeriod -> None,  
          SystemsModelLabels -> None]  
 
(*Find the impulse response of the continuous time system*)  
 
continouseTimeResponse=Chop@First@OutputResponse[ctf,  
                    DiracDelta[t],{t,0,maxSimulationTime}]  
 
Out[6]= (19.767574172152145*HeavisideTheta[t]*  
  Sin[0.2837941092083292*t])/E^(0.6931471805599452*t) +  
  (1.255589621432896*(0.9999999999999996*Cos[0.2837941092083292*t]*  
   HeavisideTheta[t] - 2.4424297688685117*HeavisideTheta[t]*  
  Sin[0.2837941092083292*t]))/E^(0.6931471805599452*t)  
 
(*Plot the impulse response of the discrete system*)  
 
 
ListPlot[  
discreteResponse,  
DataRange->{0,maxSimulationTime},  
Filling->Axis,  
FillingStyle->{Red,  
               PlotMarkers->Graphics[{PointSize[0.03],  
               Point[{0,0}]}]},PlotRange->All,  
ImageSize->300,  
ImagePadding->{{45,10},{35,10}},  
AspectRatio->0.9,  
AxesOrigin->{0,0},  
Frame->True,  
Axes->None,  
ImageSize->300,  
Frame->True  
]

PIC

(*Plot the impulse response of the continuous system*)  
Plot[samplePeriod *continouseTimeResponse,{t,0,maxSimulationTime},  
             PlotRange->All,  
             Frame->True,  
             ImageSize->300]

PIC

(*Plot both responses on the same plot*)  
 
Show[  
  ListPlot[  
    discreteResponse,  
    Filling->Axis,  
    FillingStyle->{Red,PlotMarkers->Graphics[  
                {PointSize[0.03],Point[{0,0}]}]},PlotRange->All,  
    DataRange->{0,maxSimulationTime},  
    ImageSize->300,  
    ImagePadding->{{45,10},{35,50}},  
    AspectRatio->0.9,  
    AxesOrigin->{0,0},  
    Frame->True,  
    Axes->None  
    ],  
  Plot[samplePeriod *continouseTimeResponse,  
       {t,0,maxSimulationTime},PlotRange->All  
],  
PlotRange->All,  
FrameLabel->{{Style["response",10],None},  
            {Style["time (sec)",10],  
            "impulse response"}}  
]

PIC

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

PIC

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

PIC

30 Find the system type given an open loop transfer function

Problem: Find the system type for the following transfer functions

1.
s+1s2 s
2.
s+1s3 s2
3.
s+1 s5

To find the system type, the transfer function is put in the form kissi sMjssj, then the system type is the exponent M. Hence it can be seen that the first system above has type one since the denominator can be written as s1 s1 and the second system has type 2 since the denominator can be written as s2 s1 and the third system has type 5. The following computation determines the type

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

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

31 Find the eigenvalues and eigenvectors of a matrix

Problem, given the matrix

1 2 3 4 5 6 7 8 9

Find its eigenvalues and eigenvectors.

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

PIC

N[Eigenvalues[a]]//MatrixForm

PIC

N[Eigenvectors[a]]//MatrixForm

PIC

Matlab

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

clear all; close all;  
 
A=[1 2 3;4 5 6;7 8 9];  
 
[v,e]=eig(A)  
 
>>  
v =  
   -0.2320   -0.7858    0.4082  
   -0.5253   -0.0868   -0.8165  
   -0.8187    0.6123    0.4082  
 
e =  
   16.1168         0         0  
         0   -1.1168         0  
         0         0   -0.0000  

32 Find the characteristic polynomial of a matrix

Problem, given the matrix

1 2 3 4 5 6 7 8 0

Find its characteristic polynomial.

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

Matlab

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

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

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

Problem, given the matrix

1 2 3 2

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

Mathematica
Remove["Global`*"]  
a = {{1,2},{3,2}};  
n = Length[a];  
p = CharacteristicPolynomial[a,x]  
 
Out[15]= -4-3 x+x^2  
 
(-4 IdentityMatrix[n] - 3 a + MatrixPower[a,2])//MatrixForm

PIC

Another way to do the above is as follows

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

Matlab

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

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

34 Solve Ax = b and display the solution

Problem: Solve for x given that Ax = b where

A = 1 2 1 3

and

b = 5 8

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

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

PIC

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

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

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

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

The following algorithm summarizes all the cases

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

Let the system of equations be

y = x1 y = 2x+1

So

A = 1 1 2 1

and

b = 1 1

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

PIC

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

The above algorithm can now be run as follows

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

The output of the above is

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

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

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

Problem: Given

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

Automatically convert it to the form Ax = b and solve

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

PIC

LinearSolve[a,-b]//N  
 
Out[27]= {6.71429,-2.85714,-3.14286}

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

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

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

PIC

MatrixForm[RowReduce[mat]]

PIC

Matlab
clear all;  
A=[1 1 -2 1  
   3 2 4  -4  
   4 3 3  -4];  
rref(A)  
 
>>  
ans =  
 
     1     0     0     2  
     0     1     0    -3  
     0     0     1    -1

38 Check for stability of a continuous system represented as a transfer function and state space

Problem: Given a system Laplace transfer function, check if it is stable, then convert to state space and check stability again.

In transfer function representation, the check is that all poles of the transfer function (or the zeros of the denominator) have negative real part.

In state space, the check is that the matrix A is negative definite. This is done by checking that all the eigenvalues of the matrix A have negative real part.

The poles of the transfer function are the same as the eigenvalues of the A matrix.

Use

sys = 5s s2 +4s+25

Checking stability using transfer function poles

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

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

Checking stability using state space A matrix

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

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

39 Check continuous system stability in the Lyapunov sense

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

A = 0 1 0 0 0 1 1 2 3

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

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

PIC

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

PIC

N[Eigenvalues[p]]  
Out[39]= {5.52534,1.78449,0.190164}  
 
r = Select[Re[N[Eigenvalues[p]]],# <= 0 &]  
Out[42]= {}  
 
If[Length[r]==0,  
  Print["System is asymptotically stable"],  
  Print["System is not asymptotically stable"]  
]  
 
   System is asymptotically stable

Matlab
clear all;  
 
A=[0   1   0  
   0   0   1  
   -1 -2  -3];  
 
p=lyap(A,eye(length(A)))  
>>  
p =  
 
    5.0000   -0.5000   -1.5000  
   -0.5000    1.5000   -0.5000  
   -1.5000   -0.5000    1.0000  
 
e=eig(p)  
>>  
e =  
 
    0.1902  
    1.7845  
    5.525  
 
 
indx=find(real(e)<=0)  
>>  
indx =  
 
   Empty matrix: 0-by-1  
 
 
if isempty(indx)  
    fprintf('System is asymptotically stable\n');  
else  
    fprintf('System is not asymptotically stable\n');  
end  
>>  
    System is asymptotically stable  

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

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

PIC System block diagram.

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

PIC

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

PIC

sys=SystemsModelSeriesConnect[d,plantd]

PIC

loopBack=SystemsModelFeedbackConnect[sys]

PIC

Now generate unit step response

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

PIC

ss=StateSpaceModel[loopBack]

PIC

poles=TransferFunctionPoles[loopBack]  
Out[69]= {{{0.543991 -0.325556 I,0.543991 +0.325556 I}}}  
 
Abs[#]&/@poles  
Out[70]= {{{0.633966,0.633966}}}

Poles are inside the unit circle, hence stable.

Matlab
clear all; close all  
 
s      = tf('s');  
plant  = s/(s+1);  
T      = 1; %sampeling time;  
plantd = c2d(plant,T,'zoh');  
d      = c2d( (s+1)/(s+5), T , 'zoh');  
sys    = series(d,plantd);  % obtain the open loop  
loop   = feedback(sys,1)    % obtain the closed loop  
>>  
 
loop =  
 
   z^2 - 1.801 z + 0.8013  
  ------------------------  
  2 z^2 - 2.176 z + 0.8038  
 
Sample time: 1 seconds  
Discrete-time transfer function.  
 
 
step(loop)  
grid  
set(gcf,'Position',[10,10,320,320]);

PIC

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

41 Check if a matrix is Hermite

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

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

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

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

PIC

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

PIC

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

Map[Abs[#]<$MachineEpsilon&,diff,{2}]  
 
Out[32]= {{True,True,True},  
          {True,True,True},  
          {True,True,True}}  
 
Count[%,False]  
Out[33]= 0

Matlab
clear all;  
i=sqrt(-1);  
mat=[-1     1-3*i    0;  
    1+3*i    0      -6*i;  
      0      6*i     1];  
 
mat2=conj(transpose(mat));  
diff = mat-mat2  
>>  
diff =  
 
     0     0     0  
     0     0     0  
     0     0     0  
 
r=abs(diff)<eps('double')  
>>  
r =  
 
     1     1     1  
     1     1     1  
     1     1     1  
 
 
all(r(:))  
>>  
ans =  
 
     1  

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

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

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

PIC

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

PIC

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

PIC

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

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

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

PIC

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

PIC

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

PIC

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

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

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

PIC

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

PIC

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

PIC

45 Draw the root locus from the open loop transfer function

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

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

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

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

PIC

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

PIC

46 Find the cross correlation between two sequences

Problem: Given

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]

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

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

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

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

Matlab

In MATLAB use the xcross in the signal processing toolbox

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

47 Find the different norms of a vector

Problem: Given the vector say

v = 1,2,4

Find its norm for p = 1,2,

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

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

48 Find orthonormal vectors that span the range of matrix A

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

A = 0 1 1 2 1 2 3 4 2 0 2 0

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

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

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

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

49 Plot the surface described by f y,x

Problem: Plot the function

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

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

PIC

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

PIC

50 Find the point of intersection of 3 surfaces

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

Let the 3 surfaces be

z = 2x+y z = 2y z = 2

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

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

PIC

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

PIC

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

PIC

PIC

51 Obtain the LU decomposition of a matrix

Problem: Given a matrix A, find matrix L and U such that LU = A

The matrix L will have 1 in all the diagonal elements and zeros in the upper triangle. The matrix U will have 0 in the lower triangle as shown in this diagram.

PICLU decomposition.

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

PIC

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

PIC

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

PIC

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

PIC

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

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

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

Mathematica
Remove["Global`*"];  
fname = "stereo.wav";  
SetDirectory[ToFileName[Extract  
  ["FileName"/.NotebookInformation[  
   EvaluationNotebook[]],{1},FrontEnd`FileName]]];  
 
ele = Import[fname,"Elements"]  
 
Out[110]= {AudioChannels,AudioEncoding,  
   Data,SampledSoundList,SampleRate,Sound}  
 
(*listen to it *)  
sound = Import[fname,"Sound"];  
EmitSound[sound]  
 
fs = Import[fname,"SampleRate"]  
 
Out[113]= 22050  
 
(* now load the samples and do power spectrum *)  
data = Import[fname,"Data"];  
{nChannel,nSamples}=Dimensions[data]  
 
Out[115]= {2,29016}  
 
py = Fourier[data[[1,All]],  
             FourierParameters->{1,-1}];  
nUniquePts = Ceiling[(nSamples+1)/2]  
 
Out[117]= 14509  
 
py = py[[1;;nUniquePts]];  
py = Abs[py];  
py = py/nSamples;  
py = py^2;  
 
If[ OddQ[nSamples],  
  py[[2;;-1]]=2*py[[2;;-1]],  
  py[[2;;-2]]=2*py[[2;;-2]]  
];  
 
f = (Range[0,nUniquePts-1] fs)/nSamples;  
ListPlot[Transpose[{f,py}],  
  Joined->True,  
  FrameLabel->{{"|H(f)|",None},  
               {"hz","Magnitude spectrum"}},  
  ImageSize->400,  
  Frame->True,  
  RotateLabel->False,  
  GridLines->Automatic,  
  GridLinesStyle->Dashed,  
  PlotRange->{{0,12000},All}]

PIC

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

PIC

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

Problem: Solve

yt = 3y t

with initial condition y 0 = 1 and plot the solution for t = 02 seconds.

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

PIC

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

PIC

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

Problem: Solve

yt1.5yt+5y t = 0

with initial conditions

y 0 = 1,y0 = 0

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

Given yt1.5yt+5y t = 0

let

x1 = y x2 = y = x1 hence
x1 = x2

and

x2 = y = 1.5y5y = 1.5x2 5x1

Hence we can now write

x1 x2 = 0 1 5 1.5 x1 x2

Now Matlab ODE45 can be used.

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

PIC

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

PIC

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

Problem: Solve

yt1.5yt+5y t = 4cos t

with initial conditions

y 0 = 1,y0 = 0

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

Given yt1.5yt+5y t = 4cos t, let

x1 = y x2 = y = x1

hence

x1 = x2

and

x2 = y = 1.5y5y+4cos t = 1.5x2 5x1 +4cos t

Hence we can now write

x1 x2 = 0 1 5 1.5 x1 x2 + 0 4 cos t

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

PIC

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

PIC

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

Problem: Solve

yt+ty t = 0

with the boundary conditions

y 0 = 3,y 20 = 1

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

Given yt+ty t = 0, let

x1 = y x2 = y = x1

Therefore

x1 = x2

And

x2 = y = ty = tx1

This results in

x1 x2 = 0 1 t 0 x1 x2

Now bvp4c() can be used.

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

PIC

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

PIC

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

The PDE is

T x,t t = k2T x,t x2

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

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

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

PIC

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

PIC

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

PIC

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

The PDE is

T x,t t = k2T x,t x2

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

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

Each plot shows the boundary conditions used.

PIC

59 Generate uniform distributed random numbers

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

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

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

Fortran
program t3  
implicit none  
real :: x(5)  
 
 CALL RANDOM_SEED()  
 CALL random_number(x)  
 print *,x  
 
end program

compile and run

$ gfortran -std=f95 -Wextra -Wall -pedantic -funroll-loops  
  -ftree-vectorize -march=native  -Wsurprising -Wconversion  
  t3.f90  /usr/lib/liblapack.a /usr/lib/libblas.a  
$ ./a.exe  
  0.99755955      0.56682467      0.96591532      0.74792767      0.36739087  

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

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

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

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

Fortran
program t3_2  
implicit none  
integer ::i  
real, parameter :: a=-1.0,b=1.0  
real :: x(5)  
 
 CALL RANDOM_SEED()  
 DO i=1,2  
    CALL random_number(x)  
    x = a+(b-a)*x  
    print *,x  
 END DO  
 
end program

compile and run

$ gfortran -std=f95 -Wextra -Wall -pedantic -funroll-loops  
  -ftree-vectorize -march=native  -Wsurprising -Wconversion  
  t3_2.f90  /usr/lib/liblapack.a /usr/lib/libblas.a  
$ ./a.exe  
  0.99511909      0.13364935      0.93183064      0.49585533     -0.26521826  
 -3.87262106E-02 -0.85249150     -0.98928964     -0.30583751     -0.31551242  
$

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

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

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

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

Fortran
program t3_3  
implicit none  
integer ::i  
real, parameter :: a=-1.0,b=1.0  
real :: x(5,5)  
 
 CALL RANDOM_SEED()  
 DO i=1,2  
    CALL random_number(x)  
    x = a+(b-a)*x  
    CALL write(x)  
 END DO  
 
 !--- internal functions below ------  
     contains  
       SUBROUTINE write(A)  
       implicit none  
       REAL, DIMENSION(:,:) :: A  
       integer :: i,j  
 
       WRITE(*,*)  
       DO i = lbound(A,1), ubound(A,1)  
          WRITE(*,*) (A(i,j), j = lbound(A,2), ubound(A,2))  
       END DO  
 
  END SUBROUTINE write  
 
end program

compile and run

$ gfortran -std=f95 -Wextra -Wall -pedantic -funroll-loops  
    -ftree-vectorize -march=native  -Wsurprising -Wconversion  
    t3_3.f90  /usr/lib/liblapack.a /usr/lib/libblas.a  
 
$ ./a.exe  
 
  0.99511909     -3.87262106E-02 -0.56409657      0.32386434      0.71138477  
  0.13364935     -0.85249150     -0.73367929     -0.96778345     -0.19742620  
  0.93183064     -0.98928964      0.80104899      0.30170965     -0.58625138  
  0.49585533     -0.30583751     -0.22646809      0.29281759      0.93707883  
 -0.26521826     -0.31551242     -0.10903549     -0.35402548      0.19679904  
 
  0.34596145      0.21138644      0.22462976      0.31809497      0.45771694  
 -8.62354040E-02  0.43809581      0.95732045      0.10801017     -0.19508958  
 -0.33996975      0.79466915      0.99828446      0.95552015      0.85725522  
 -0.79923415      0.31645823     -0.48640406      0.80384660     -0.70432973  
  0.51090658     -0.69856644      0.10173070      0.31584930      0.34905851  

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

Plot the Fourier Transform of 3sin(t)

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

PIC

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

Need to do the plot.

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

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

PIC

System boundary conditions.

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

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

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

PIC

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

PIC

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

Problem: Generate one signal of some width and height.

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

PIC

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

PIC

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

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

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

PIC

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

PIC

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

Given the sequence of numbers x(n) = 1,2,3, Find X k = m=0N1x mej2π Nmk where N is the length of the data sequence x m and k = 0N 1

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

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

Ada

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

The compiler used, and the output is shown below.

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

Fortran

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

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

65 Find eAt where A is a matrix

In state space one of the most important operations is the determination of eAt where A is the system matrix.

Mathematica
ClearAll["Global`*"]  
mat={{0,1},  
     {-2,-3}}  
 
MatrixExp[mat t]  
 
Out[6]= {{problems/e^(-2 t) (-1+2 E^t),E^(-2 t) (-1+E^t)},  
         {-2 E^(-2 t) (-1+E^t),-E^(-2 t) (-2+E^t)}}  
 
 
MatrixForm[%]

PIC

Now verify the above result by solving for eAt using the method would one would do by hand, if a computer was not around.

There are a number of methods to do this by hand. The eigenvalue method, based on the Cayley Hamilton theorem will be used here.

Find the eigenvalues of |A- lambdaI|

m = mat-lambda IdentityMatrix[Length[mat]]  
 
Out[9]= {{-lambda,  1        },  
         {-2,       -3-lambda}}  
 
Det[m]  
 
Out[14]= 2+3 lambda+lambda^2  
 
sol=Solve[%==0,lambda]  
 
Out[15]= {{lambda->-2},{lambda->-1}}  
 
eig1=lambda/.sol[[1]]  
eig2=lambda/.sol[[2]]  
 
Out[16]= -2  
Out[17]= -1  
 
(*setup the equations to find b0,b1*)  
 
eq1 = Exp[eig1 t]==b0+b1 eig1;  
eq2 = Exp[eig2 t]==b0+b1 eig2;  
sol = First@Solve[{problems/eq1,eq2},{b0,b1}]  
 
Out[20]= {b0->E^(-2 t) (-1+2 E^t),  
          b1->-E^(-2 t)+E^-t}  
 
(*Now find e^At*)  
 
b0=b0/.sol[[1]]  
 
Out[24]= E^(-2 t) (-1+2 E^t)  
 
b1=b1/.sol[[2]]  
 
Out[25]= -E^(-2 t)+E^-t  
 
b0 IdentityMatrix[Length[mat]]+b1 mat;  
Simplify[%]  
 
Out[27]= {{problems/e^(-2 t) (-1+2 E^t),E^(-2 t) (-1+E^t)},  
          {-2 E^(-2 t) (-1+E^t),-E^(-2 t) (-2+E^t)}}  
 
In[28]:= MatrixForm[%]

PIC

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

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

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

Mathematica and Matlab 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 list of lists. This means rows can have different sizes. (these are the lists inside the list).

Hence each list (or row) can be different size. In Matlab a matrix must have same size rows and so a matrix in mathematica ofocurse, since that is the defintion of a matrix, but the point is that one can have rows with different sizes (The result is no longer called a matrix).

In Matlab one can also have ragged arrays, but one uses cells for that, and can't use the standard matrix. In Mathematica, there is one data structure for both, it is all the same data structure a List of lists.

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 we want to find the location where A(i,j) = 2 where i is the row number and j is the column number. Hence given this matrix

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

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

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

From help

help find  
   [I,J] = find(X,...) returns the row and column indices

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

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

which gives

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

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

Both Matlab and Mathematica give the same result, but Mathematica is row-wise oriented while Matlab is column-wise.

Mathematica use for matrix manipulate will take more time to master compared to Matlab, since Mathematica data struct is more general and little more complex (ragged arrays) compared to Matlab since it is more general and supports symbolics and not just numbers, but at the end of the day, one can do pretty much the same thing in each system as the examples below illustrate

66.1 Multiply matrix with a vector

How to perform the following matrix/vector multiplication?

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

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

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

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

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

compile and run

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

Fortran
program t1  
implicit none  
 
integer, parameter :: N=3  
real(8), parameter, DIMENSION(N, N) :: &  
           A=DBLE(reshape([ 1.0, 2.0, 3.0,&  
                       4.0, 5.0, 6.0,&  
                       7.0, 8.0, 9.0], shape(A), order=[2,1]))  
 
real(8), parameter, DIMENSION(N, 1) :: &  
           x=DBLE(reshape([1.0, 2.0, 3.0], shape(x)))  
 
real(8), DIMENSION(N, 1) :: result  
integer, DIMENSION(2)    :: siz  
 
result = MATMUL(A,x)  
 
Write(*,*) result  
siz = SHAPE(result)  
print *, "number of rows = ",siz(1)  
print *, "number of columns = ",siz(2)  
 
end program

compile and run

$ gfortran -std=f2003 -Wextra -Wall -pedantic -funroll-loops  
  -ftree-vectorize -march=native  -Wsurprising -Wconversion  
  t1.f90  /usr/lib/liblapack.a /usr/lib/libblas.a  
$ ./a.exe  
   14.000000000000000    32.000000000000000    50.000000000000000  
 number of rows =            3  
 number of columns =            1

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

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

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

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

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

66.3 Insert a row into a matrix

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

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

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

Fortran
program t3  
  implicit none  
  integer i,nRow,nCol  
  integer(4), ALLOCATABLE ::A(:,:)  
 
  A = reshape ([1,2,3,   &  
                4,5,6,   &  
                7,8,9], [3,3], order=[2,1])  
  nRow=size(A,1)  
  nCol=size(A,2)  
 
  print *, 'before'  
  do i=LBOUND(A, 1),UBOUND(A, 1)  
      print *,A(i,:)  
  end do  
 
  A = reshape ([A(1,:),    &  
                90,91,92,  &  
                A(2:,:)], [nRow+1,nCol] , order=[2,1])  
 
  print *, 'after'  
  do i=LBOUND(A, 1),UBOUND(A, 1)  
     print *,A(i,:)  
  end do  
 
end program t3  
 
 
>gfortran -std=f2008 t3.f90  
>./a.out  
 before  
           1           2           3  
           4           5           6  
           7           8           9  
 after  
           1           2           3  
          90          91          92  
           4           7           5  
           8           6           9

66.4 Insert a column into a matrix

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

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

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

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

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

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

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

Fortran

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

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

66.6 Generate an n by m zero matrix

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

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

66.7 Rotate a matrix by 900

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

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

66.8 Generate a diagonal matrix with given values on the diagonal

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

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

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

66.9 Sum elements in a matrix along the diagonal

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

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

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

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

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

66.11 Check if a Matrix is diagonal

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

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

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

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

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

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

66.13 Replicate a matrix

Given Matrix

     1     2  
     3     4

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

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

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

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

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

66.14 Find the location of the maximum value in a matrix

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

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

66.15 Swap 2 columns in a matrix

Give a matrix

     1     2  
     3     4  
     5     6

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

     2     1  
     4     3  
     6     5

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

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

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

Mathematica

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

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

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

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

Question posted on the net

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

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

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

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

66.18 problems/extract values from matrix given their index

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

Example, given

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

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

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

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

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

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

66.19 Convert N by M matrix to a row of length NM

Given

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

covert the matrix to one vector

[1 2 3 4 5 6 7 8 9]

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

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

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

Example, given

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

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

      1     9     0    10  
      3     9     2    10

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

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

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

Given

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

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

     [6 7 2]

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

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

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

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

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

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

66.22 Locate rows in a matrix with column being a string

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

For example, given

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

and given list

{'foo', 'faa'}

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

Hence the result of the above should be

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

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

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

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

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

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

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

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

Assume valid list of indices.

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

PIC

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

Mathematica

Three methods are shown.

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

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

pos = {2, 4};

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

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

method 2: (due to Mr Wizard at stackoverflow)

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

method 3: (me)

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

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

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

66.24 Convert list of separated numerical numbers to strings

Problem: given a list of numbers such as

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

convert the above to list of strings

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

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

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

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

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

66.25 Obtain elements that are common to two vectors

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

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

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

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

Given

     4     2     5  
     2     7     9  
    10     1     2

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

     2     1     2  
     4     2     5  
    10     7     9

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

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

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

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

Given

     4     2     5  
     2     7     9  
    10     1     2

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

     2     4     5  
     2     7     9  
     1     2     10

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

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

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

Given

     4     2     5  
     2     7     9  
    10     1     2

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

     2     7     9  
     4     2     5  
    10     1     2

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

Mathematica
mat={{4, 2, 5},  
     {2, 7, 9},  
     {10,1, 2}};  
Sort[mat]  
 
Out[182]= {{2,  7, 9},  
           {4,  2, 5},  
           {10, 1, 2}}  

Matlab
A=[4  2 5;  
   2  7 9;  
   10 1 2];  
sortrows(A)  
 
 
ans =  
     2     7     9  
     4     2     5  
    10     1     2  

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

Given

     4     2     5  
     2     7     9  
    10     1     2

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

    10     1     2  
     4     2     5  
     2     7     9

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
mat={{4, 2, 5},  
     {2, 7, 9},  
     {10,1, 2}};  
SortBy[mat,#[[2]]&]  (*sort by second element as key*)  
 
Out[191]= {{10, 1, 2},  
           {4,  2, 5},  
           {2,  7, 9}}  

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

66.30 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 = 50 75 0 50 0 100 0 75 100 75 100 0 0 75 100 0 75 100 will become B = 1 75 0 1 0 100 0 1 100 1 100 0 0 1 100 0 1 100

Mathematica
A={ {50, 75,  0},  
    {50, 0,   100},  
    {0,  75,  100},  
    {75, 100, 0},  
    {0,  75,  100},  
    {0,  75,  100}};  
 
ReplacePart[A,DeleteDuplicates[Position[A,_?(#!=0&)],  
      (#1[[1]]==#2[[1]]&)]->-1]  
 
Out[9]= {{-1,75,0},  
         {-1,0,100},  
         {0,-1,100},  
         {-1,100,0},  
         {0,-1,100},  
         {0,-1,100}}

Solution due to Bob Hanlon (from Mathematica newsgroup)

r=Union[Position[A,_?(# !=0&)],SameTest->(#1[[1]]==#2[[1]]&)];  
ReplacePart[A,r->-1]

Solution by Fred Simons (from Mathematica newsgroup)

r=Split[Position[A,_?(# !=0&),{2}],#1[[1]]==#2[[1]]&][[All,1]];  
ReplacePart[A,r->-1]

Solution due to Adriano Pascoletti (from Mathematica newsgroup)

r=(First /@ Split[#1,First[#1]===First[#2]&]&)[Position[A,x_/;x!=0]]  
ReplacePart[A,r->-1]

Solution due to Andrzej Kozlowski (from Mathematica newsgroup)

r=MapIndexed[Flatten[{#2,Position[#1,x_/;x!=0,1,1]}]&,A]  
ReplacePart[A,r->-1]

Solution due to Oliver Ruebenkoenig (from Mathematica newsgroup)

r=First /@ SplitBy[ Drop[ArrayRules[SparseArray[A]][[All,1]],-1],First]  
ReplacePart[A,r->-1]

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

r=MapIndexed[{First[#2],1+LengthWhile[#1,#==0&]}&,A]  
ReplacePart[A,r->-1]

Matlab
A=[50 75 0;  
   50 0 100;  
   0 75 100;  
   75 100 0;  
   0 75 100;  
   0 75 100];  
 
[I,J] = ind2sub(size(A),find(A~=0));  
[b,c] = unique(I,'first');  
A(sub2ind(size(A),b,J(c)))=-1  
 
A =  
 
    -1    75     0  
    -1     0   100  
     0    -1   100  
    -1   100     0  
     0    -1   100  
     0    -1   100  

This solution due to Bruno Luong (from matlab newsgroup)

B = cumsum(A~=0,2)>0  
B = [false(size(B,1),1) B]  
A(logical(diff(B,1,2))) = -1

This solution due to Jos from matlab newsgroup

A((cumsum(~~A,2)==1) & (A ~= 0)) = -1

66.31 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

In[43]:= Remove["Global`*"]  
 
v1={a,b,c};  
v2={problems/e,f,g};  
Outer[Times,v1,v2]  
 
Out[46]= {{a e, a f, a g},  
          {b e, b f, b g},  
          {c e, c f, c g}}

Outer sum

Outer[Plus,v1,v2]  
 
Out[48]= {{a+e, a+f, a+g},  
          {b+e, b+f, b+g},  
          {c+e, c+f, c+g}}

using numerical vectors. Outer product

v1={1,2,3};  
v2={4,5,6};  
Outer[Times,v1,v2]  
 
Out[52]= {{4,  5,  6},  
          {8,  10, 12},  
          {12, 15, 18}}

Outer sum

Outer[Plus,v1,v2]  
 
  {{5, 6, 7},  
   {6, 7, 8},  
   {7, 8, 9}}  

Matlab

Outer product

v1=[1 2 3];  
v2=[4 5 6];  
 
v1'*v2  
 
ans =  
     4     5     6  
     8    10    12  
    12    15    18

Outer sum

v1=[1 2 3];  
v2=[4 5 6];  
 
meshgrid(v1)+meshgrid(v2)'  
 
ans =  
     5     6     7  
     6     7     8  
     7     8     9

Maple

Due to Carl Love from the Maple newsgroup

restart;  
>  
> Outer:= proc(OP, L1, L2)  
>    local a,b;  
>    [seq](seq(OP(a,b), b in L2), a in L1)  
> end proc:  
>  
> L1:=[a,b,c];  
> L2:=[d,e,f];  
> Outer(`*`,L1,L2);  
 
[a*d, a*e, a*f, b*d, b*e, b*f, c*d, c*e, c*f]  
 
> Outer(`+`,L1,L2);  
 
[a+d, a+e, a+f, b+d, b+e, b+f, c+d, c+e, c+f]  

66.32 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 = 2 3 3 4 0 1 2 2 0 0 0 1

Mathematica
mat={{2,3,3,4},  
     {0,-1,-2,2},  
     {0,0,0,1}};  
 
{nRow,nCol} = Dimensions[mat]  
 
Out[14]= {3,4}  
 
Print["Rank (or dimension of the range space)=",  
       MatrixRank[mat]]  
 
   Rank (or dimension of the range space)=3  
 
 
Print["Dimension of the Null Space=",  
      nCol-MatrixRank[mat]]  
 
   Dimension of the Null Space=1  
 
Print["Basis for Null Space of mat =",NullSpace[mat]]  
 
   Basis for Null Space of mat ={{3,-4,2,0}}

Malab
A=[2 3 3 4;  
   0 -1 -2 2;  
   0 0 0 1]  
 
[nRow,nCol]=size(A);  
 
r = rank(A);  
fprintf('Dimension of range space of A=%d\n',r);  
fprintf('Dimension of the null space = %d\n',nCol-r);  
fprintf('Basic for null space of A =');  
null(A,'r')'  
 
>>>  
Dimension of range space of A=3  
Dimension of the null space = 1  
Basic for null space of A =  
ans =  
 
    1.5000   -2.0000    1.0000         0  

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

Problem: Find the SVD for the matrix 1 2 3 4 5 6

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
mat = {{1,2,3},  
       {4,5,6}}  
 
{u,s,v}=N[SingularValueDecomposition[mat]];  
u  
 
Out[175]= {{0.386318,-0.922366},  
           {0.922366,0.386318}}  
 
s  
 
Out[176]= {{9.50803,0.,0.},  
           {0.,0.77287,0.}}  
 
v  
 
Out[177]= {{0.428667,0.805964,0.408248},  
           {0.566307,0.112382,-0.816497},  
           {0.703947,-0.581199,0.408248}}  
 
(*Reconstruct A from its components *)  
 
u.s.Transpose[v]  
 
Out[178]= {{1.,2.,3.},  
           {4.,5.,6.}}  

Matlab
A=[1 2 3;  
   4 5 6];  
[u,s,v]=svd(A);  
u  
s  
v  
>>>  
 
u =  
 
   -0.3863   -0.9224  
   -0.9224    0.3863  
 
s =  
 
    9.5080         0         0  
         0    0.7729         0  
 
v =  
 
   -0.4287    0.8060    0.4082  
   -0.5663    0.1124   -0.8165  
   -0.7039   -0.5812    0.4082  
 
 
 
u*s*v'  
>>  
ans =  
 
    1.0000    2.0000    3.0000  
    4.0000    5.0000    6.0000  

Maple
> restart;  
> with(LinearAlgebra):  
> A:=Matrix([[1,2,3],[4,5,6.]]);  
                                   [1    2    3 ]  
                              A := [            ]  
                                   [4    5    6.]  
 
> m,n:=Dimensions(A):  
> u,s,v:=SingularValues(A, output=['U', 'S', 'Vt']):  
> u;  
                   [0.386317703118612    -0.922365780077058]  
                   [                                       ]  
                   [0.922365780077058    0.386317703118612 ]  
 
> s:=DiagonalMatrix(s,m,n);  
                    [9.50803200069572            0            0]  
               s := [                                          ]  
                    [       0            0.772869635673485    0]  
 
> v;  
        [0.428667133548626    0.566306918848035     0.703946704147444 ]  
        [                                                             ]  
        [0.805963908589298    0.112382414096594     -0.581199080396110]  
        [                                                             ]  
        [0.408248290463863    -0.816496580927726    0.408248290463863 ]  
 
> s2:=DiagonalMatrix(SingularValues(A, output=['S']),m,n);  
                    [9.50803200069572            0            0]  
              s2 := [                                          ]  
                    [       0            0.772869635673485    0]  
 
> evalf(u.s.v);  
                               [1.    2.    3.]  
                               [              ]  
                               [4.    5.    6.]  
 
>

66.34 Solve Ax = b

Solve for x in the following system of equations

[ 1 2 3 ][x1]   [ 1 ]  
[ 7 5 6 ][x2]=  [ 2 ]  
[ 7 8 8 ][x3]   [ 3 ]  

Mathematica
mat={{1,2,3},  
      {7,5,6},  
      {7,8,9}};  
 
b = {1,2,3}  
LinearSolve[mat,b]  
 
Out[58]= {0,0,1/3}  

Matlab
format long  
A=[1 2 3;  
   7 5 6;  
   7 8 9];  
b=[1;2; 3];  
A\b  
 
ans =  
                   0  
                   0  
   0.333333333333333

Fortran
program t2  
implicit none  
integer, parameter :: N=3  
real(8),   DIMENSION(N, N) :: &  
           A=DBLE(reshape([ 1.0, 2.0, 3.0,&  
                       7.0, 5.0, 6.0,&  
                       7.0, 8.0, 9.0], shape(A), order=[2,1]))  
 
real(8), parameter, DIMENSION(N, 1) :: &  
           b=DBLE(reshape([1.0, 2.0, 3.0], shape(b)))  
 
real(8), DIMENSION(N, 1) :: IPIV  
integer :: INFO  
 
 
 CALL DGETRF( N, N, A, N, IPIV, INFO )  
 if (INFO .eq. 0 ) then  
     CALL DGETRS( 'N', N,1, A, N, IPIV, b , N, INFO)  
     if ( INFO .eq. 0 ) then  
        print *, 'solution is',b  
     else  
        print *, 'failed DGETRS, INFO=',INFO  
     end if  
 else  
     print *, 'failed DGETRF, INFO=',INFO  
 end if  
end program

compile and run

$ gfortran -std=f2003 -Wextra -Wall -pedantic -funroll-loops  
  -ftree-vectorize -march=native  -Wsurprising -Wconversion  
  t2.f90  /usr/lib/liblapack.a /usr/lib/libblas.a  
 
$ ./a.exe  
 solution is  
 -5.28677630773884192E-018 -3.70074341541718826E-017  0.33333333333333337  

66.35 Find all nonzero elements in a matrix

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

0 0 1 10 0 2 3 0 0

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.

mat  = {{0 , 0, 1},  
        {10, 0, 2},  
        {3 , 4, 0}};  
 
sp  = SparseArray[mat];  
pos = sp["NonzeroPositions"]  
 
Out[13]= {{1,3},{2,1},{2,3},{3,1},{3,2}}  
 
values=sp["NonzeroValues"]  
 
Out[14]= {1,10,2,3,4}

Or standard list operations can be used

mat  = {{0 , 0, 1},  
        {10, 0, 2},  
        {3 , 4, 0}};  
 
(*find values not zero *)  
 
Cases[mat,x_ /; Not[PossibleZeroQ[x]] :>x ,2]  
 
Out[34]= {1,10,2,3,4}  
 
(*find the index *)  
 
Position[mat,x_/;x!=0]  
 
Out[45]= {{1,3},{2,1},{2,3},{3,1},{3,2}}

Matlab
A=[0 0 1; 10 0 2; 3 4 0];  
values= nonzeros(A)  
 
>>>  
values =  
 
    10  
     3  
     4  
     1  
     2  
 
 
%find locations  
[I,J]=find(A)  
 
>>>  
I =  
 
     2  
     3  
     3  
     1  
     2  
 
 
J =  
 
     1  
     1  
     2  
     3  
     3  

66.36 problems/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) = x2 evaluate it on (1,2,3) such that the result is (1,4,9).

Mathematica
Clear[f,x]  
f[x_]:=x^2  
v={1,2,3};  
f[v]  
 
Out[45]= {1,4,9}

Matlab
clear all;  
v=[1 2 3]  
f=@(x) x.^2  
f(v)  
 
ans =  
 
     1     4     9

66.37 generates equally spaced N points between x1 and x2

Mathematica
x1=1;  
x2=3;  
FindDivisions[{x1,x2},5]  
 
Out[48]= {1,3/2,2,5/2,3}  
 
 
N[%]  
Out[49]= {1.,1.5,2.,2.5,3.}

Matlab
clear all;  
x1 = 1;  
x2 = 3;  
linspace(x1,x2,5)  
 
ans =  
    1.0000    1.5000    2.0000    2.5000    3.0000

66.38 problems/evaluate and plot a f(x,y) on 2D grid of coordinates

Evaluate xexpx2y2 on 2D cartesian grid between x = 22 and y = 44 using h = 0.4 for grid spacing.

Mathematica
f[x_,y_]:=x Exp[-x^2-y^2];  
data = Table[f[x,y],{x,-2,2,.2},{y,-4,4.2}];  
ListPlot3D[data,  
           PlotRange->All,  
           ColorFunction->"Rainbow",  
           AxesLabel->{"x","y","z"},  
           LabelStyle->12,  
           Ticks->{Range[-2,2,1],Range[-4,4,1],Automatic},  
           DataRange->{{-2,2},{-4,4},Automatic},  
           ImageSize->400]

PIC

The above can also be done using Plot3D

f[x_,y_]:=x Exp[-x^2-y^2];  
Plot3D[f[x,y],{x,-2,2},{y,-4,4},  
       PlotRange->All,  
       ColorFunction->"Rainbow",AxesLabel->{"x","y","z"},  
       LabelStyle->12,  
       ImageSize->400]

PIC

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

Matlab
[X,Y] = meshgrid(-2:.2:2, -4:.2:4);  
 Z = X .* exp(-X.^2 - Y.^2);  
 surf(X,Y,Z)  
 xlabel('x');  
 ylabel('y');  
 zlabel('z');

PIC

66.39 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
mat=Table[RandomReal[],{3},{3}]  
Det[mat]  
 
Out[15]= -0.317792605942287

Matlab
A=rand(3,3);  
det(A)  
 
ans =  
 
  -0.276653966272994

66.40 Generate sparse matrix with n by n matrix repeated on its diagonal

Given matrix 1 2 3 4 5 6 7 8 9 , generate the following sparse matrix with this matrix on the diagonal

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

Mathematica
mat = N[{{1,2,3},{4,5,6},{7,8,9}}];  
sp  = SparseArray[Band[{1,1}]->ConstantArray[mat,{3}]]  
 
Out[69]= SparseArray[<27>,{9,9}]  
 
MatrixForm[sp]

PIC

Matlab
A  = [1 2 3;  
      4 5 6;  
      7 8 9];  
Iy = speye(3);  
full(Iy)  
>>>  
ans =  
 
     1     0     0  
     0     1     0  
     0     0     1  
>>>  
 
sp = kron(Iy,A);  
full(sp)  
 
>>>  
ans =  
 
     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

66.41 Generate sparse matrix for the tridiagonal representation of the classic second difference operator on 1D grid

The second derivative d2u dx2 is approximated by ui12ui+ui+1 h2 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
numberOfInternalGridPoints = 4;  
n  = 3*internalGridPoints;  
sp = SparseArray[{Band[{1,1}]->-2,  
                  Band[{1,2}]->1,  
                  Band[{2,1}]->1},  
                  {n,n}]  
 
Out[103]= SparseArray[<34>,{12,12}]  
 
MatrixForm[sp]

PIC

Matlab
internalGridPoints = 4;  
n = 3*internalGridPoints;  
e = ones(n,1);  
sp = spdiags([e -2*e e], -1:1, n,n);  
 
full(sp)  
 
>>>  
ans =  
 
    -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

66.42 Generate sparse matrix for the Laplacian differential operator 2u for 2D grid

2u = f in 2D is defined as 2u x2 + 2u y2 = f and the Laplacian operator using second order standard differences results in ui1,j+ui+1,j+ui,j1+ui,j+14ui,j h2 = fi,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.

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

Mathematica
n  = 3;  
sp = SparseArray[{Band[{1,1}]->-4,  
                  Band[{1,2}]->1,  
                  Band[{2,1}]->1,  
                  Band[{4,1}]->1,  
                  Band[{1,4}]->1},  
                  {n^2,n^2}];  
 
Out[2]= SparseArray[<37>,{9,9}]  
 
MatrixForm[sp]

PIC

Matlab
internalPoints=3;  
e   = ones(internalPoints,1);  
spe = spdiags([e -2*e e], -1:1,internalPoints,internalPoints);  
Iz  = speye(internalPoints);  
sp  = kron(Iz,spe)+kron(spe,Iz);  
 
full(sp)  
 
>>>  
ans =  
 
    -4     1     0     1     0     0     0     0     0  
     1    -4     1     0     1     0     0     0     0  
     0     1    -4     0     0     1     0     0     0  
     1     0     0    -4     1     0     1     0     0  
     0     1     0     1    -4     1     0     1     0  
     0     0     1     0     1    -4     0     0     1  
     0     0     0     1     0     0    -4     1     0  
     0     0     0     0     1     0     1    -4     1  
     0     0     0     0     0     1     0     1    -4  
 

66.43 Generate sparse matrix for the Laplacian differential operator 2u for 3D grid

The goal is to solve

2u x2 + 2u y2 + 2u z2 = 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

PIC

The discrete approximation to the Laplacian in 3D is

2u x2 + 2u y2 + 2u z2 = 1 h2 Ui1,j,k+Ui+1,j,k+Ui,j1,k+Ui,j+1,k+Ui,k,k1 +Ui,j,k+1 6Ui,j,k

For the direct solver, the A matrix needs to be formulated. From

1 h2 Ui1,j,k+Ui+1,j,k+Ui,j1,k+Ui,j+1,k+Ui,k,k1 +Ui,j,k+1 6Ui,j,k = fi,j,k

Solving for Ui,j,k results in

Ui,j,k = 1 6 Ui1,j,k+Ui+1,j,k+Ui,j1,k+Ui,j+1,k+Ui,k,k1 +Ui,j,k+1 h2f i,j,k

To help make the A matrix, a small example with n = 2, is made. The following diagram uses the standard numbering on each node

PIC

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

PIC

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,

PIC

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

PIC

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 n2 higher than the grid point itself. n2 is the size of the grid at each surface. Hence, the general A matrix, for the above example, can now be written as

PIC

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 nx = 2, ny = 2, nz = 2. These are the number of grid points in the x,y,z directions respectively.

Mathematica
In[23]:= nx=2;  
ny=2;  
nz=2;  
sp=SparseArray[{Band[{1,1}]->-6,  
                Band[{1,2}]->1,  
                Band[{2,1}]->1,  
                Band[{nx+1,1}]->1,  
                Band[{1,ny+1}]->1,  
                Band[{nx*ny+1,1}]->1,  
                Band[{1,nx*ny+1}]->1},  
                {nx*ny*nz,nx*ny*nz}]  
 
Out[26]= SparseArray[<42>,{8,8}]  
 
Normal[sp]  
 
Out[28]= {{-6,1,1,0,1,0,0,0},  
          {1,-6,1,1,0,1,0,0},  
          {1,1,-6,1,1,0,1,0},  
          {0,1,1,-6,1,1,0,1},  
          {1,0,1,1,-6,1,1,0},  
          {0,1,0,1,1,-6,1,1},  
          {0,0,1,0,1,1,-6,1},  
          {0,0,0,1,0,1,1,-6}}  
 
MatrixForm[sp]

PIC

Matlab
nx=2;  
ny=2;  
ex=ones(nx,1);  
Lx=spdiags([ex -3*ex ex],[-1 0 1],nx,nx);  
ey=ones(ny,1);  
Ly=spdiags([ey -3*ey ey],[-1 0 1],ny,ny);  
 
Ix=speye(nx);  
Iy=speye(ny);  
L2=kron(Iy,Lx)+kron(Ly,Ix);  
 
nz=2;  
N=nx*ny*nz;  
e=ones(N,1);  
L=spdiags([e e],[-nx*ny nx*ny],N,N);  
Iz=speye(nz);  
 
A=kron(Iz,L2)+L;  
full(A)  
 
ans =  
 
    -6     1     1     0     1     0     0     0  
     1    -6     0     1     0     1     0     0  
     1     0    -6     1     0     0     1     0  
     0     1     1    -6     0     0     0     1  
     1     0     0     0    -6     1     1     0  
     0     1     0     0     1    -6     0     1  
     0     0     1     0     1     0    -6     1  
     0     0     0     1     0     1     1    -6  

Example 2: Using nx = 2, ny = 2, nz = 3.

Mathematica
nx=2;  
ny=2;  
nz=3;  
sp=SparseArray[{Band[{1,1}]->-6,  
                Band[{1,2}]->1,  
                Band[{2,1}]->1,  
                Band[{nx+1,1}]->1,  
                Band[{1,ny+1}]->1,  
                Band[{nx*ny+1,1}]->1,  
                Band[{1,nx*ny+1}]->1},  
                {nx*ny*nz,nx*ny*nz}]  
 
Out[38]= SparseArray[<70>,{12,12}]  
 
Normal[sp]  
 
Out[39]= {{-6,1,1,0,1,0,0,0,0,0,0,0},  
          {1,-6,1,1,0,1,0,0,0,0,0,0},  
          {1,1,-6,1,1,0,1,0,0,0,0,0},  
          {0,1,1,-6,1,1,0,1,0,0,0,0},  
          {1,0,1,1,-6,1,1,0,1,0,0,0},  
          {0,1,0,1,1,-6,1,1,0,1,0,0},  
          {0,0,1,0,1,1,-6,1,1,0,1,0},  
          {0,0,0,1,0,1,1,-6,1,1,0,1},  
          {0,0,0,0,1,0,1,1,-6,1,1,0},  
          {0,0,0,0,0,1,0,1,1,-6,1,1},  
          {0,0,0,0,0,0,1,0,1,1,-6,1},  
          {0,0,0,0,0,0,0,1,0,1,1,-6}}  
 
MatrixForm[sp]

PIC

Matlab
nx=2;  
ny=2;  
ex=ones(nx,1);  
Lx=spdiags([ex -3*ex ex],[-1 0 1],nx,nx);  
ey=ones(ny,1);  
Ly=spdiags([ey -3*ey ey],[-1 0 1],ny,ny);  
 
Ix=speye(nx);  
Iy=speye(ny);  
L2=kron(Iy,Lx)+kron(Ly,Ix);  
 
nz=3;  
N=nx*ny*nz;  
e=ones(N,1);  
L=spdiags([e e],[-nx*ny nx*ny],N,N);  
Iz=speye(nz);  
 
A=kron(Iz,L2)+L;  
full(A)  
 
ans =  
 
    -6     1     1     0     1     0     0     0     0     0     0     0  
     1    -6     0     1     0     1     0     0     0     0     0     0  
     1     0    -6     1     0     0     1     0     0     0     0     0  
     0     1     1    -6     0     0     0     1     0     0     0     0  
     1     0     0     0    -6     1     1     0     1     0     0     0  
     0     1     0     0     1    -6     0     1     0     1     0     0  
     0     0     1     0     1     0    -6     1     0     0     1     0  
     0     0     0     1     0     1     1    -6     0     0     0     1  
     0     0     0     0     1     0     0     0    -6     1     1     0  
     0     0     0     0     0     1     0     0     1    -6     0     1  
     0     0     0     0     0     0     1     0     1     0    -6     1  
     0     0     0     0     0     0     0     1     0     1     1    -6

Example 3: Using nx = 3, ny = 3, nz = 3.

Mathematica

nx = 3;  
ny = 3;  
nz = 3;  
sp = SparseArray[{Band[{1, 1}] -> -6,  
                  Band[{1, 2}] -> 1,  
                  Band[{2, 1}] -> 1,  
                  Band[{nx + 1, 1}] -> 1,  
                  Band[{1, ny + 1}] -> 1,  
                  Band[{nx*ny + 1, 1}] -> 1,  
                  Band[{1, nx*ny + 1}] -> 1},  
                  {nx*ny*nz,nx*ny*nz}];  
 
MatrixForm[sp]

PIC

Matlab

nx=3;  
ny=3;  
ex=ones(nx,1);  
Lx=spdiags([ex -3*ex ex],[-1 0 1],nx,nx);  
ey=ones(ny,1);  
Ly=spdiags([ey -3*ey ey],[-1 0 1],ny,ny);  
 
Ix=speye(nx);  
Iy=speye(ny);  
L2=kron(Iy,Lx)+kron(Ly,Ix);  
 
nz=3;  
N=nx*ny*nz;  
e=ones(N,1);  
L=spdiags([e e],[-nx*ny nx*ny],N,N);  
Iz=speye(nz);  
 
A=kron(Iz,L2)+L;  
full(A)  
 
<display not shown, same as Mathematica output above>  

66.44 Generate the adjugate matrix for square matrix

Given square matrix such as

1 2 3 4 5 6 7 8 9

find the adjugate matrix which is

3 6 3 6 12 6 3 6 3

Mathematica
<<"Combinatorica`"  
mat      = {{1,2,3},{4,5,6},{7,8,9}};  
cof      = Table[Cofactor[mat,{i,j}],{i,3},{j,3}];  
adjugate = Transpose[cof]  
 
Out[22]= { {-3,   6, -3},  
           { 6, -12,  6},  
           {-3,   6, -3}}

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.

A=[1 2 3;4 5 6;7 8 9]  
det(A)*inv(A)  
 
>>>  
Warning: Matrix is close to singular or badly scaled.  
Results may be inaccurate. RCOND = 1.541976e-18.  
 
ans =  
 
   -3.0000    6.0000   -3.0000  
    6.0000  -12.0000    6.0000  
   -3.0000    6.0000   -3.0000

The following is due to Matt J from the matlab newsgroup

A=[1 2 3;  
   4 5 6;  
   7 8 9]  
 
c=poly(A);  
adjunct_A=polyvalm(c(1:end-1),A)  
 
>>>  
adjunct_A =  
 
   -3.0000    6.0000   -3.0000  
    6.0000  -12.0000    6.0000  
   -3.0000    6.0000   -3.0000

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.

!-- Find the Adjugate of a square matrix  
!-- Version 6/22/2012 by Nasser M. Abbasi  
!-- gfortran 4.6.3  
program t44  
implicit none  
 
integer, parameter :: n=3  
integer :: INFO,i,j,ii,IPIV(n-1),number_row_exchange  
real (kind=kind(0.0d0)) :: A(n,n),B(n-1,n-1),C(n,n)  
logical :: M(n,n)  
 
A(1,:) = [1, 2, 3];  
A(2,:) = [4, 5, 6];  
A(3,:) = [7, 8, 9];  
 
DO j=1,n  
   DO i=1,n  
 
      M = .true.  
      M(:,j) = .false.  
      M(i,:) = .false.  
      B = reshape(pack(A, M),[n-1,n-1])  
 
      !-- LU decomposition in order to obtain determinant  
      CALL DGETRF( n-1, n-1, B, n-1, IPIV, INFO )  
 
      !-- determinant of U is the product of diagonal  
      C(i,j)= PRODUCT( [(B(ii,ii),ii=1,n-1)] )  
 
      !-- Adjust sign based on number of row exchanges and (-1)^(i+j)  
      number_row_exchange = count(IPIV /= [(ii,ii=1,n-1)])  
      C(i,j) = C(i,j)*(1-2*MODULO(number_row_exchange+i+j,2))  
 
   END DO  
END DO  
write(*,'(3F15.4)') C  
end program t44

compile, link and run

export LD_LIBRARY_PATH=/usr/lib/atlas-base/:/usr/lib/atlas-base/atlas/  
gfortran -fcheck=all -Wall t44.f90 -L/usr/lib/atlas-base/atlas/ -lblas -llapack  
 
>./a.out  
        -3.0000         6.0000        -3.0000  
         6.0000       -12.0000         6.0000  
        -3.0000         6.0000        -3.0000  
 
!>dpkg -l|grep -E "(blas|lapack)"  
!ii  libblas3gf        1.2.20110419-2ubuntu1    Basic Linear Algebra Reference implementations, shared library  
!ii  liblapack-dev     3.3.1-1                  library of linear algebra routines 3 - static version  
!ii  liblapack-doc     3.3.1-1                  library of linear algebra routines 3 - documentation  
!ii  liblapack3gf      3.3.1-1                  library of linear algebra routines 3 - shared version  
!ii  libopenblas-base  0.1alpha2.2-3            Optimized BLAS (linear algebra) library based on GotoBLAS2  
!ii  libopenblas-dev   0.1alpha2.2-3            Optimized BLAS (linear algebra) library based on GotoBLAS2  
!>

Ada

with Ada.Text_Io; use Ada.Text_Io;  
with Ada.Float_Text_Io; use Ada.Float_Text_Io;  
with Ada.integer_Text_Io; use Ada.integer_Text_Io;  
with Ada.Numerics.Real_Arrays;  use Ada.Numerics.Real_Arrays;  
 
 
procedure t44 is  
 
    A : constant real_matrix :=  
             (( 1.0,  2.0,  3.0),  
              ( 4.0,  5.0,  6.0),  
              ( 7.0,  8.0,  9.0));  
           --(( -3.0,  2.0,  -5.0),  
           --   ( -1.0,  0.0,  -2.0),  
           --   ( 3.0,  -4.0,  1.0));  
    C : real_matrix(A'range(1), A'range(2));  
    B : real_matrix(1 .. 2, 1 .. 2);  
 
    -- Thanks goes to Dmitry A. Kazakov this function  
    function Exclude (A : Real_Matrix; I, J : Integer)  
       return Real_Matrix is  
       AI, AJ : Integer := A'First (1);  
    begin  
       return B : Real_Matrix (1..A'Length (1)-1, 1..A'Length (2)-1) do  
          AI := A'First (1);  
          for BI in B'Range (1) loop  
             if AI = I then  
                AI := AI + 1;  
             end if;  
             AJ := A'First (2);  
             for BJ in B'Range (2) loop  
                if AJ = J then  
                   AJ := AJ + 1;  
                end if;  
                B (BI, BJ) := A (AI, AJ);  
                AJ := AJ + 1;  
             end loop;  
             AI := AI + 1;  
          end loop;  
       end return;  
    end Exclude;  
 
    -- Function to print a 2D matrix  
    procedure Put (X : Real_Matrix) is  
       begin  
          for I in X'Range (1) loop  
             for J in X'Range (2) loop  
                Put (X (I, J) ); put("    ");  
             end loop;  
             New_Line;  
          end loop;  
   end Put;  
 
begin  
 
    FOR I in A'range(1) LOOP  
      FOR J in A'range(2) LOOP  
          B := Exclude (A, I, J);  
          C(I,J) := float((-1)**( ((I+J) rem 2) ))*determinant(B);  
      END LOOP;  
    END LOOP;  
    c := transpose(c);  
    put(c);  
 
end t44;

compile and run

-- export LD_LIBRARY_PATH=/usr/lib/atlas-base/:/usr/lib/atlas-base/atlas/  
>gnatmake t44.adb  
gcc-4.6 -c t44.adb  
gnatbind -x t44.ali  
gnatlink t44.ali  
 
>./t44  
-3.00000E+00     6.00000E+00    -3.00000E+00  
 6.00000E+00    -1.20000E+01     6.00000E+00  
-3.00000E+00     6.00000E+00    -3.00000E+00  

66.45 Multiply each column by values taken from a row

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

PIC

In Matlab, the bsxfun is used.

Mathematica

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

r   = {2,3};  
mat = {{1,2},{3,4},{5,6}};  
r #&/@mat  
 
Out[12]= {{2,  6},  
          {6,  12},  
          {10, 18}}  
 
Or  
 
Map[v*#&, mat]  

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

r={2,3};  
mat={{1,2},{3,4},{5,6}};  
Inner[Times,mat,r,List]  
 
Out[66]= {{2,  6},  
          {6,  12},  
          {10, 18}}

Matlab
r=[2 3];  
A=[1 2;  
   3 4;  
   5 6]  
bsxfun(@times,r,A)  
 
 
ans =  
 
     2     6  
     6    12  
    10    18  

Fortran
program t45  
 implicit none  
 
 integer, parameter :: nRow=3, nCol=2  
 
 character(len=*), parameter :: FMT = "(2I5)"  
 integer :: v(nCol),A(nRow,nCol),B(nRow,nCol),i,j;  
 
 !-- initialization of data  
 v = [2,3]  
 A(1,:) = [1,2];  
 A(2,:) = [3,4];  
 A(3,:) = [5,6];  
 
 !F2008  
 !do concurrent (j = 1:nCol)  
 !   B(:,j) = v(j)*A(:,j)  
 !end do  
 
 do j=1,nCol  
    B(:,j) = v(j)*A(:,j)  
 end do  
 
 write(*,FMT) ( (B(i,j),j=1,nCol), i=1,nRow)  
 
end program t45

compile and run

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

Octave

Use automatic broadcasting

r .* A

66.46 problems/extract submatrix from a larger matrix by removing row/column

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

For example, given the matrix 1 2 3 4 5 6 7 8 9 then the submatrix for 1,1 is 5 6 8 9 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.

Mathematica
mat={{1,2,3},  
     {4,5,6},  
     {7,8,9}};  
{nRow, nCol} = Dimensions[mat];  
Table[ReplacePart[mat,{{i},{i,_},{_,j}}  
      :>Sequence[],Heads-> False],{i,nRow},{j,nCol}];  
 
r=Flatten[%,1]  
 
 
Out[182]= {{{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}}}

Matlab
A=[1 2 3;  
   4 5 6;  
   7 8 9]  
 
for i=1:size(A,1)  
    for j=1:size(A,1)  
        B=A;  
        B(i,:)=[];  
        B(:,j)=[]  
    end  
end  
 
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

Fortran
program t46  
implicit none  
 
integer, parameter:: dp=kind(0.d0), n=3  
!integer, parameter :: n=3  
integer :: i,j,ii  
real(dp) :: A(n,n) = dble(reshape([ 1.0, 2.0, 3.0,  &  
                                    4.0, 5.0, 6.0,  &  
                                    7.0, 8.0, 9.0], &  
                                    shape(A),       &  
                                    order=[2,1]))  
real(dp) :: B(n-1,n-1)  
logical :: mask(n,n)  
 
DO i=1,n  
   DO j=1,n  
 
      mask = .true.  
      mask(:,j) = .false.  
      mask(i,:) = .false.  
 
      !-- extract submatrix  
      B = reshape(pack(A, mask),[n-1,n-1])  
 
      DO ii=1,n-1  
         write(*,'(2F6.1)') B(ii,:)  
      END DO  
      write(*,'(/)')  
   END DO  
END DO  
end program t46

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

Ada
with Ada.Text_Io; use Ada.Text_Io;  
with Ada.Float_Text_Io; use Ada.Float_Text_Io;  
with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;  
 
procedure foo is  
    A : constant Real_Matrix :=  
             (( 1.0,  2.0,  3.0),  
              ( 4.0,  5.0,  6.0),  
              ( 7.0,  8.0,  9.0));  
    B : real_matrix(1 .. 2, 1 .. 2);  
 
    -- Thanks goes to Dmitry A. Kazakov  
    -- for providing this function  
    function Exclude (A : Real_Matrix; I, J : Integer)  
       return Real_Matrix is  
       AI, AJ : Integer := A'First (1);  
    begin  
       return B : Real_Matrix(1..A'Length(1)-1,1..A'Length(2)-1) do  
          AI := A'First (1);  
          for BI in B'Range (1) loop  
             if AI = I then  
                AI := AI + 1;  
             end if;  
             AJ := A'First (2);  
             for BJ in B'Range (2) loop  
                if AJ = J then  
                   AJ := AJ + 1;  
                end if;  
                B (BI, BJ) := A (AI, AJ);  
                AJ := AJ + 1;  
             end loop;  
             AI := AI + 1;  
          end loop;  
       end return;  
    end Exclude;  
 
    -- Function to print a 2D matrix  
    procedure put(A: Real_Matrix) is  
    begin  
     for I in A'range(1) loop  
      for J in A'range(2) loop  
          put(A(I,J));  
      end loop;  
      new_line;  
     end loop;  
    end put;  
 
begin  
 
    FOR I in A'range(1) LOOP  
      FOR J in A'range(2) LOOP  
          B := Exclude (A, I, J);  
          put(B);  
          new_line(1);  
      END LOOP;  
    END LOOP;  
 
end foo;

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

66.47 delete one row from a matrix

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

1 2 3 4 5 6 7 8 9

Mathematica
mat={{1,2,3},  
     {4,5,6},  
     {7,8,9}};  
 
mat[[2]]=Sequence[];  
mat  
 
Out[81]= {{1,2,3},  
          {7,8,9}}

or a little longer solution using Pick

{nRow,nCol}=Dimensions[mat];  
mask=Table[True,{nRow},{nCol}];  
mask[[2]]=False;  
Pick[mat,mask]  
 
Out[70]= {{1,2,3},  
          {7,8,9}}

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

Maple
A:=<1,2,3;  
    4,5,6;  
    7,8,9>;  
                                    [1  2  3]  
                                    [       ]  
                               A := [4  5  6]  
                                    [       ]  
                                    [7  8  9]  
 
A:=LinearAlgebra:-DeleteRow(A,2);  
                                    [1  2  3]  
                               A := [       ]  
                                    [7  8  9]

66.48 delete one column from a matrix

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

1 2 3 4 5 6 7 8 9

Mathematica
mat={{1,2,3},  
     {4,5,6},  
     {7,8,9}};  
 
mat[[All,2]]=Sequence[];  
mat  
 
>>>>  
Out[93]= {{1,3},  
          {4,6},  
          {7,9}}  

or a little longer solution using Pick

{nRow,nCol}=Dimensions[mat];  
mask=Table[True,{nRow},{nCol}];  
mask[[All,2]]=False;  
mat=Pick[mat,mask];  
mat  
 
>>>>>  
Out[98]= {{1,3},  
          {4,6},  
          {7,9}}  

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

Maple
A:=<1,2,3;  
    4,5,6;  
    7,8,9>;  
                                    [1  2  3]  
                                    [       ]  
                               A := [4  5  6]  
                                    [       ]  
                                    [7  8  9]  
 
A:=LinearAlgebra:-DeleteColumn(A,2);  
                                      [1  3]  
                                      [    ]  
                                 A := [4  6]  
                                      [    ]  
                                      [7  9]

66.49 generate random matrix so that each row adds to 1

Generate the random matrix and divide each row by its total

Mathematica
mat = Table[RandomReal[],{3},{4}]  
s   = Total[mat,{2}]  
b   = mat/s  
Total[b,{2}]  
 
>>>>>  
Out[21]= {{0.57575,0.64264,0.293957,0.870987},  
          {0.534529,0.0291316,0.23062,0.180575},  
          {0.891017,0.812128,0.751813,0.417362}}  
 
Out[22]= {2.38333,0.974855,2.87232}  
 
Out[23]= {{0.241573,0.269639,0.123339,0.365449},  
          {0.548316,0.029883,0.236568,0.185233},  
          {0.310208,0.282743,0.261744,0.145305}}  
 
Out[24]= {1.,1.,1.}  

Matlab
A = rand(3,4)  
s = sum(A,2);  
B = bsxfun(@rdivide,A,s)  
sum(B,2)  
 
>>>>  
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

66.50 generate random matrix so that each column adds to 1

Generate the random matrix and divide each column by its total

Mathematica

This method due to Bob Hanlon from the Math group

mat =  Table[RandomReal[],{3},{4}]  
s   = Total[mat,{1}]  
b   = #/s& /@ mat  
Total[b,{1}]  
 
>>>>  
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

b   = Transpose[Transpose[mat]/s];  
Total[b,{1}]  
Out[37]= {1.,1.,1.,1.}

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

b=Inner[Divide,mat,s,List];  
Total[b,{1}]  
 
Out[71]= {1.,1.,1.,1.}

Matlab
A=rand(3,4)  
s=sum(A,1);  
B=bsxfun(@rdivide,A,s)  
sum(B,1)  
 
>>>>>>  
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

Maple
A:=LinearAlgebra:-RandomMatrix(3,4);  
                                [ 9  -95   51  24]  
                                [                ]  
                           A := [99  -20   76  65]  
                                [                ]  
                                [60  -25  -44  86]  
 
b:=MTM:-sum(A,1);  #sum columns  
                          b := [168, -140, 83, 175]  
 
nCol:=LinearAlgebra:-ColumnDimension(A):  
B:=convert([seq(A[..,i]/b[i],i=1..nCol)],Matrix);  
 
                                [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]  
MTM:-sum(B,1);  
                                [1, 1, 1, 1]

66.51 sum all rows in a matrix

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

1 2 3 4 5 6 7 8 9

Mathematica
mat={{1,2,3},  
      {4,5,6},  
      {7,8,9}}  
Total[mat,{2}]  
 
 
Out[2]= {6,  
        15,  
        24}  

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

66.52 sum all columns in a matrix

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

1 2 3 4 5 6 7 8 9

Mathematica
In[3]:= mat={{1,2,3},  
             {4,5,6},  
             {7,8,9}}  
Total[mat,{1}]  
 
Out[4]= {12,15,18}  

Matlab
A=[1 2 3;4 5 6;7 8 9];  
sum(A,1)  
 
 
ans =  
    12    15    18  

66.53 find in which columns values that are not zero

given matrix 0 39 0 55 100 0 34 0 0 0 9 0 0 0 50 find the column index which contain values > 0 but do it from top row to the bottom row. Hence the result should be 21strow,1,22ndrow,1,2,3 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
mat={{0, 39,  0},  
     {55,100, 0},  
     {34,0,   0},  
     {0, 9,   0},  
     {0, 0,   50}};  
Position[#,x_/;x>0]&/@mat;  
Flatten[%]  
 
 
Out[18]= {2,1,2,1,2,3}

Matlab
A=[0  39   0  
   55 100  0  
   34 0    0  
   0  9    0  
   0  0   50];  
[I,~]=find(A'>0)  
I'  
 
ans =  
     2     1     2     1     2     3  

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

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.

Mathematica
1.
first method
v1={1,3,4,5,8,9,20.2,30,-44};  
v2={1,7,4};  
v1=DeleteCases[v1,Alternatives@@v2]  
 
Out[61]= {3,5,8,9,20.2,30,-44}

2.
second method
v1={1,3,4,5,8,9,20.2,30,-44};  
v2={1,7,4};  
v1=DeleteCases[v1,x_/;MemberQ[v2,x]]  
 
Out[67]= {3,5,8,9,20.2,30,-44}

Matlab
v1=[1,3,4,5,8,9,20.2,30,-44];  
v2=[1,7,4];  
v1(ismember(v1,v2))=[]  
 
 
v1 =  
 3.0000  5.0000 8.0000 9.0000  20.2000 30.0000 -44.0000

Fortran
program f  
  implicit none  
  REAL, ALLOCATABLE :: A(:),B(:)  
  A =[1.0,3.0,4.0,5.0,8.0,9.0,20.2,30.0,-44.0];  
  B =[1.0,7.0,4.0];  
  A = pack(A, A .NE. B)  
  Print *,A  
end program f  
 
>gfortran f2.f90  
>./a.out  
   3.0000000 5.0000000 8.0000000  9.0000000  20.200001  30.000000 -44.000000

66.55 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
len = 52; n = 4;  
v   = RandomReal[{0, 1}, len];  
Mean /@ Partition[v, n]  
 
{0.750653,0.410073,0.401005,  
0.138907,0.466247,0.568257,  
0.535362,0.517755,0.555368,  
0.705857,0.502319,0.453571,0.357949}

Matlab
N = 52;  
n = 4;  
V = rand(N,1);  
mean(reshape(V,n,N/n))'  
 
ans =  
 
    0.3026  
    0.3589  
    0.4745  
    0.6249  
    0.3042  
    0.5428  
    0.2387  
    0.3200  
    0.6224  
    0.4408  
    0.5657  
    0.5469  
    0.5696  

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

Given matrix

1 5 2 3 4 8 7 2

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

1 5 2 3

Mathematica
a = {{1, 5}, {2, 3}, {4, 8}, {7, 2}};  
Min@Position[a, {x_, y_} /; y > 6]  
(* 3 *)  
 
a[[1 ;; % - 1]]  
(* {{1, 5}, {2, 3}} *)

Matlab
A = [1,5;2,3;4,8;7,2];  
A(1:find(A(:,2)>6)-1,:)  
 
ans =  
      1     5  
      2     3

67 Display spectrum of 2D image

Mathematica
img=Import["ExampleData/lena.tif"];  
Image[img,ImageSize->300]

PIC

ImageDimensions[img]  
 
Out[38]= {150,116}  
 
ImageChannels[img]  (*see how many channels*)  
 
Out[39]= 3  
 
data=ImageData[img];  (*get data*)  
{nRow,nCol,nChannel}=Dimensions[data]  
 
Out[41]= {116,150,3}  
 
(*look at each channel*)  
Map[Image[data[[All,All,#]],ImageSize->100]&,  
         Range[1,nChannel]]

PIC

(*get channel 2 to FFT but center it first*)  
d = data[[All,All,2]];  
d = d*(-1)^Table[i+j,{i,nRow},{j,nCol}];  
 
(*make FFT,center,  and look at spectrum and phase*)  
 
fw = Fourier[d,FourierParameters->{1,1}];  
 
(*adjust for better viewing as needed*)  
fudgeFactor = 100;  
abs = fudgeFactor * Log[1+Abs@fw];  
Labeled[Image[abs/Max[abs],  
       ImageSize->300],  
       Style["Magnitude spectrum", 18]]

PIC

arg = Arg@fw;  
Labeled[Image[arg/Max[arg],  
        ImageSize->300],  
        Style["Phase spectrum", 18]]

PIC

Matlab
close all; clear all;  
%image in same folder  
img = imread('lena.tif','tif');  
imagesc(img)

PIC

img = fftshift(img(:,:,2));  
F = fft2(img);  
figure;  
imagesc(100*log(1+abs(fftshift(F)))); colormap(gray);  
title('magnitude spectrum');

PIC

figure;  
imagesc(angle(F));  colormap(gray);  
title('phase spectrum');

PIC

68 Draw root locus for a discrete system

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

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

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

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

PIC

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

PIC

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

PIC

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

PIC

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

PIC

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

PIC

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

PIC

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

PIC

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

In[168]:= sys=sys/.{g->9.8,m->1,M->20,L->1};  
y = Chop[OutputResponse[sys,UnitStep[t],t]]  
 
Out[169]= {problems/e^(-6.41561 t) (0.000115693 E^(3.2078 t) UnitStep[t]-  
           0.000231385 E^(6.41561 t) UnitStep[t]+  
           0.000115693 E^(9.62341 t) UnitStep[t]+  
           0.0238095 E^(6.41561 t) t^2 UnitStep[t]),  
 
           E^(-6.41561 t) (-0.00242954 E^(3.2078 t) UnitStep[t]+  
           0.00485909 E^(6.41561 t) UnitStep[t]-  
           0.00242954 E^(9.62341 t) UnitStep[t])}

The following obtained the response as interpolatingFunction since the full range of time domain is now specified

sys = sys/.{g->9.8,m->1,M->20,L->1};  
y   = OutputResponse[sys,UnitStep[t],{t,0,3}]  
 
Out[148]= {InterpolatingFunction[{{0.,3.}},<>][t],  
           InterpolatingFunction[{{0.,3.}},<>][t]}  
 
Plot[ y[[1]],{t,0,3},PlotLabel->"y(t) response",  
      FrameLabel->{"time (sec)","y (meter)"},  
      Frame->True,  
      PlotRange->All,  
      ImageSize->300,  
      GridLines->Automatic,  
      GridLinesStyle->Dashed,  
      RotateLabel->False]

PIC

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

PIC

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

PIC

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

PIC

70 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

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

71 Solve 2nd order ODE (Van Der Pol) and generate phase plot

Problem: Solve

ytμ 1y2yt+y t = 0

for different values of μ to compare effect of changing μ on the solution trajectories. The initial conditions are

y 0 = 0.1 yt = 0

In both Mathematica and Matlab, numerical ODE solver was used.

For Matlab, The 2nd order ODE is first converted to 2 first order ODE's, and then solve the coupled ODE's using ode45. In Mathematica NDSolve was used. This does not require the conversion.

Starting by writing the 2nd order ODE as 2 first order ODE's

x1 = y x2 = y derivatives x1 = y x2 = y x1 = x2 x2 = μ 1x12x2 +x1

Below is the solution of the differential equation for different values of μ = 1

Mathematica
process[mu_]:=Module[{sol,p1,p2,p3,data,system,z,x1,  
                      x2,t,ode1,ode2,timeScale},  
ode1 =x1'[t]==x2[t];  
ode2 =x2'[t]==z(1-x1[t]^2)x2[t]-x1[t];  
timeScale ={t,0,80};  
 
sol=First@NDSolve[{ode1,ode2/.z->mu,x1[0]==0.01,x2[0]==0},  
                  {x1[t],x2[t]},  
                  timeScale,Method->{"Extrapolation",  
                    Method->"LinearlyImplicitEuler"}];  
 
sol = {x1[t],x2[t]}/.sol;  
 
p1 = Plot[Evaluate[sol[[1]]],Evaluate[timeScale],  
          Frame->True,  
          FrameLabel->{"time","y(t)"},  
          PlotLabel->"\[Mu]="<>ToString[mu] ];  
 
p2 = Plot[Evaluate[sol[[2]]],Evaluate[timeScale],  
          Frame->True,FrameLabel->{"time","y'(t)"},  
          PlotLabel->"\[Mu]="<>ToString[mu],  
          PlotStyle->RGBColor[1,0,0]];  
 
data = Table[{problems/evaluate[sol[[1]]],  
              Evaluate[sol[[2]]]},  
              Evaluate[timeScale]];  
 
p3=ListPlot[data,Joined->True,Frame->True,  
            PlotLabel->"\[Mu]="<>ToString[mu],  
            FrameLabel->{"y(t)","y'(t)"},  
            PlotRange->All];  
 
{p1,p2,p3}  
];  
mu={0.1,.2,1,3,5,10};  
r = process/@mu;  
Grid[r]

PIC

Matlab
function e71  
mu=[0.1,0.2,1,3,5,10];  
 
for n=1:length(mu)  
    process(mu(n),n-1);  
end  
 
function process(mu,n)  
timeScale=[0 80];  
ic=[0.1;0];  
 
[t,x]=ode45(@ode,timeScale,ic,[],mu);  
subplot(6,3,(3*n)+1);  
plot(t,x(:,1));  
title(sprintf('mu=%1.2f',mu));  
xlabel('time'); ylabel('y(t)');  
 
subplot(6,3,(3*n)+2);  
plot(t,x(:,2),'r');  
title(sprintf('mu=%1.2f',mu));  
xlabel('time'); ylabel('y''(t)');  
 
subplot(6,3,(3*n)+3);  
plot(x(:,2),x(:,1));  
title(sprintf('mu=%1.2f',mu));  
xlabel('y(t)'); ylabel('y''(t)');  
 
 
function xdot=ode(t,x,mu)  
xdot=[x(2) ; mu*(1-x(1)^2)*x(2)-x(1)];

PIC

72 Generate a plot using x and y data

Problem: Generate x and y data for sin(x) and plot the data.

Mathematica
Remove["Global`*"];  
data = Table[{x,Sin[x]}, {x, -2Pi,2 Pi,Pi/100}];  
ListPlot[data, Frame -> True,  
         GridLines->Automatic,  
         GridLinesStyle->Dashed,  
         ImageSize->300,  
         AspectRatio->1,  
         PlotStyle->Red,  
         FrameLabel->{{"sin(x)",None},  
                      {"x","plot of sin(x)"}},  
         FrameTicks->{{-2Pi,-Pi,0,Pi,2Pi},  
                       None,Automatic,None}]

PIC

Matlab
clear all;  
x = -2*pi:pi/100:2*pi;  
plot(x,sin(x),'r');  
grid on;  
title('plot of sin(x)');  
xlabel('x');  
ylabel('sin(x)');  
set(gca,'XTick',-2*pi:pi:2*pi)  
set(gca,'XTickLabel',{'-2*pi','-pi','0','pi','2*pi'})  
set(gcf,'Position',[10,10,340,340]);

PIC

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

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

G s = ωn2 s2 +2ζωns+ωn2

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

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

The transfer function is written as

Y s U s = ωn2 s2 +2ζωns+ωn2

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

Hence the differential equation, assuming zero initial conditions becomes

yt+2ζω nyt+ω n2y t = ω n2u t

To solve the above in Matlab using ode45, the differential equation is converted to 2 first order ODE's as was done before resulting in

x1 x2 = 0 1 ωn 2 2ζωn x1 x2 + 0 ωn2 u t

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

x1 x2 = 0 1 ωn 2 2ζωn x1 x2 + 0 ωn2

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

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

PIC

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

PIC

PIC

74 Find the linear convolution of 2 sequences where the origin is located at an arbitrary position in the sequence

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.

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

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.

(*Write down the h  and take its DTFT *)  
Clear[n,w]  
h = 1/4DiscreteDelta[n+1] + 1/2DiscreteDelta[n]+1/4 DiscreteDelta[n-1];  
hTransformed = FourierSequenceTransform[h,n,w]  
 
Out[16]= ((1/4)*(1 + E^(I*w))^2)/E^(I*w)  
 
(*Write down the x sequence and take its DTFT*)  
 
x = 1 DiscreteDelta[n]+2 DiscreteDelta[n-1]-  
    2 DiscreteDelta[n-2]+DiscreteDelta[n-3]+  
    4 DiscreteDelta[n-4]-DiscreteDelta[n-5]+  
    DiscreteDelta[n-6]+2 DiscreteDelta[n-7];  
 
xTransformed = FourierSequenceTransform[x,n,w]  
 
Out[13]= (2 + 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))/E^(7*I*w)  
 
(*Now multiply the DTFT's and take the inverse*)  
z = InverseFourierSequenceTransform [xTransformed hTransformed,w,n]

PIC

(*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*)  
 
InverseFourierSequenceTransform[xTransformed hTransformed hTransformed,w,n]//N

PIC

75 Generate direction field plot of a first order differential equation

Problem: Given the following non autonomous differential equation, plot the line fields which represents the solutions of the ODE.

dy x dx = x2 y

Direction field plot (or slope plot) shows solutions for the ODE without actually solving the ODE.

The following are the steps to generate direction field plot for dy dx = f(x,y)

1.
generate 2 lists of numbers. The y list and the x list. These 2 lists are the coordinates at which a small slop line will be plotted.
2.
At each of the above coordinates (y,x) evaluate f(x,y).
3.
Plot small line starting at the point (x,y) and having slop of f(x,y). The length of the line is kept small. Normally it has an arrow at the end of it.

Using Matlab, the points are first generated (the (x,y) coordinates) then the slope f(x,y) evaluated at each of these points, then the command quiver() is used. Next contour() command is used to plot few lines of constant slope.

In Mathematica, the command VectorPlot is used.

Mathematica
f[x_,y_]:= x^2-y  
ListVectorPlot[Table[{1,f[x,y]},{x,-2,2,0.1},{y,-2,2,0.1}],  
      FrameLabel->{{"y(x)",None},  
                   {"x","direction fields for y'(x)=x^2-y"}},  
      ImageSize->350,LabelStyle->Directive[14]]

PIC

StreamPlot[{1,f[x,y]},{x,-2,2},{y,-2,2},  
   FrameLabel->{{"y(x)",None},  
     {"x","direction fields for y'(x)=x^2-y"}},  
   ImageSize->350,  
   LabelStyle->Directive[14]]

PIC

Matlab
clear all; close all;  
 
%direction fields for dy/dx=xy  
f     = @(X,Y) X.^2 - Y;  
[X,Y] = meshgrid(-2:.4:2,-2:.4:2);  
YlengthOfVector = f(X,Y);  
XlengthOfVector = ones(size(YlengthOfVector));  
 
quiver(X,Y,XlengthOfVector,YlengthOfVector);  
xlabel('x'); ylabel('y');  
hold on;  
 
contour(X,Y,YlengthOfVector,[-5 -4 -3 -2 -1 0 1 2 3 4 5]);  
 
title(...  
 'direction fields for $\frac{dy}{dx}=x^2 - y$','interpreter',...  
 'latex','fontsize',12)

PIC

76 Obtain Fourier Series coefficients for a periodic function?

Given a continuous time function x t, find its Fourier coefficients cn, where

x t =n=c nejω0nt

and

cn = 1 T 0T 0 2 T 0 2 x(t)ejω0ntdt

Where T 0 is the fundamental period of x t and ω0 = 2π T 0

Mathematica
Clear[T,w0]  
sol = FourierSeries[Sin[2 Pi/T t],t,3,  
          FourierParameters->{1,2 Pi/T}]  
 
Out[199]= ((1/2)*I)/E^((2*I*Pi*t)/T) -  
          (1/2)*I*E^((2*I*Pi*t)/T)  
 
data=Table[{i,FourierCoefficient[Sin[2 Pi/T t],  
           t,i,FourierParameters->{1,2 Pi/T}]},{i,-5,5}];  
head = {"n","Subscript[c, n]"};  
Grid[Insert[data,head,1],Frame->All]

PIC

mag = data;  
mag[[All,2]] = Map[Abs[#]&,data[[All,2]]];  
ListPlot[mag,AxesOrigin->{0,0},  
         Filling->Axis,  
         FillingStyle->{{Thick,Red}},  
         AxesLabel->{"n","|Subscript[c, n]|"}]

PIC

phase = data;  
phase[[All,2]]=Map[Arg[#]&,data[[All,2]]]* 180/Pi;  
ListPlot[phase,AxesOrigin->{0,0},  
         Filling->Axis,  
         FillingStyle->{{Thick,Red}},  
         AxesLabel->{"n","Phase Subscript[c, n] degrees"},  
         ImageSize->300]

PIC

(*Find Fourier Series for Cos[w0 t] + Sin[w0 t]^2*)  
 
w0 = (2 Pi)/T;  
T  = 2 Pi;  
f[t_] := Cos[w0 t]+Sin[w0 t]^2;  
Plot[f[t],{t,-10 Pi,10 Pi},PlotRange->All,ImageSize -> 300]

PIC

sol = FourierSeries[f[t],t,3,  
             FourierParameters->{1,(2 Pi)/T}]  
 
Out[216]= 1/2 + 1/(E^(I*t)*2) + E^(I*t)/2 -  
          1/4/E^(2*I*t) - (1/4)*E^(2*I*t)  
 
data = Table[{i,  
          FourierCoefficient[f[t],t,i,  
              FourierParameters->{1,2 Pi/T}]},{i,-5,5}];  
 
head = {"n","Subscript[c, n]"};  
Grid[Insert[data,head,1],Frame->All]

PIC

mag = data;  
mag[[All,2]] = Map[Abs[#]&,data[[All,2]]];  
ListPlot[mag,AxesOrigin->{0,0},  
     Filling->Axis,  
     FillingStyle->{{Thick,Red}},  
     AxesLabel->{"n","|Subscript[c, n]|"},  
     ImageSize->300]

PIC

(*Find Fourier series for periodic square wave*)  
f[x_] := SquareWave[{0,1},(x+.5)/2] ;  
Plot[f[x],{x,-6,6},Filling->Axis,ImageSize->300]

PIC

T   = 2;  
sol = Chop[FourierSeries[f[t],t,6,  
           FourierParameters->{1,(2 Pi)/T}]]  
 
Out[223]= 0.5 + 0.3183098861837907/E^(I*Pi*t) +  
  0.3183098861837907*E^(I*Pi*t) -  
  0.1061032953945969/E^(3*I*Pi*t) -  
  0.1061032953945969*E^(3*I*Pi*t) +  
  0.06366197723675814/E^(5*I*Pi*t) +  
  0.06366197723675814*E^(5*I*Pi*t)  
 
data = Chop[Table[{i,FourierCoefficient[f[t],t,i,  
          FourierParameters->{1,2 Pi/T}]},{i,-8,8}]];  
 
head = {"n","Subscript[c, n]"};  
Grid[Insert[data,head,1],Frame->All]

PIC

data=Table[{i,FourierCoefficient[f[t],t,i,  
            FourierParameters->{1,2 Pi/T}]},{i,-20,20}];  
mag = data;  
mag[[All,2]]=Map[Abs[#]&,data[[All,2]]];  
ListPlot[mag,AxesOrigin->{0,0},  
         Filling->Axis,  
         FillingStyle->{{Thick,Red}},  
         PlotRange->All,  
         AxesLabel->{"n","|Subscript[c, n]|"},  
         ImageSize->300]

PIC

phase = data;  
phase[[All,2]]=Map[Arg[#]&,data[[All,2]]]* 180/Pi;  
ListPlot[phase,AxesOrigin->{0,0},  
         Filling->Axis,  
         FillingStyle->{{Thick,Red}},  
         AxesLabel->{"n","Phase Subscript[c, n] degrees"},  
         ImageSize->300]

PIC

77 Plot the constant energy levels for a nonlinear pendulum

Problem:

Plot the constant energy levels for the nonlinear pendulum in 𝜃,𝜃̇

PIC

Assume that m = 1,g = 9.8ms2,L = 10m

Answer:

The constant energy curves are curves in the Y-X plane where energy is constant. The Y-axis represents 𝜃̇, and the X-axis represents 𝜃

We assume the pendulum is given an initial force when in the initial position (𝜃 = 00) that will cause it to swing anticlock wise. The pendulum will from this point obtain an energy which will remain constant since there is no damping.

The higher the energy the pendulum posses, the larger the angle 𝜃 it will swing by will be.

If the energy is large enough to cause the pendulum to reach 𝜃 = π (the upright position) with an non zero angular velocity, then the pendulum will continue to rotate in the same direction and will not swing back and forth.

The expression for the energy E for the pendulum is first derived as a function of 𝜃,𝜃̇

E = PE +KE = mgL 1cos𝜃+ 1 2mL2𝜃̇2

The school PDF report contains more information on this topic.

This was solved in Mathematica using the ListContourPlot[] command after generating the energy values.

The original Matlab implementation is left below as well as the Maple implementation. However, these were not done using the contour method, which is a much simpler method. These will be updated later.

The following is the resulting plot for m = 1,g = 9.8ms2,L = 10m

Mathematica
energy[theta_,speed_,m_,g_,len_]:=  
    m g len(1-Cos[theta])+(1/2) m len^2 speed^2;  
 
m=1; g=9.8; len=10;  
 
data = Table[energy[i,j,m,g,len],  
             {j,-4,4,.01},  
             {i,-3 Pi,3 Pi,Pi/20.}];  
 
ListContourPlot[data,  
      InterpolationOrder->2,  
      Contours->20,  
      MaxPlotPoints->30,  
      DataRange->{{-3 Pi,3Pi},{-4,4}},  
      FrameTicks->{ {{-4,-3,-2,-0,1,2,3,4},None},  
                    {{-3 Pi,-2 Pi,-Pi,0,Pi,2Pi,3Pi},None}},  
      FrameLabel->{{"\[Theta]'(t)",None},  
        {"\[Theta](t)",  
        "constant energy levels for nonlinear pendulum model"}},  
      ImageSize->400,  
      LabelStyle->14,  
      RotateLabel->False]

PIC

Matlab
function HW2  
 
close all; clear all;  
m=1; g=9.8; L=10;  
 
nPoints=40;  %number of data points  
 
nEnergyLevels = 8; %number of energy levels  
lowAnglesToVisit = linspace(0,pi,nEnergyLevels);  
lowEnergyLevels(1:nEnergyLevels) = m*g*L*(1-cos(lowAnglesToVisit));  
highAnglesToVisit = linspace(pi,2*pi,2*nEnergyLevels);  
highEnergyLevels=zeros(1,2*nEnergyLevels);  
 
initialHighEnergy=2*m*g*L;  
Q=0.2;  
for i=1:2*nEnergyLevels  
    highEnergyLevels(i) = initialHighEnergy+(Q*i*initialHighEnergy);  
end  
 
A = zeros(length(lowAnglesToVisit)+length(highAnglesToVisit),2);  
A(:,1) = [lowAnglesToVisit highAnglesToVisit];  
A(:,2) = [lowEnergyLevels highEnergyLevels];  
 
[nAngles,~]=size(A);  
data=zeros(nPoints,2);  
 
for j=1:nAngles  
 
    currentAngle=A(j,1);  
    currentEnergyLevel =A(j,2) ;  
    angles=linspace(0,currentAngle,nPoints);  
    data(1:nPoints,1)=angles(:);  
 
    for m=1:nPoints  
        data(m,2)=speed(currentEnergyLevel,angles(m));  
    end  
    doPlots(data);  
end  
 
title(sprintf('constant energy curves, quantum=%1.3f',Q));  
xlabel('angle (position)');  
ylabel('speed');  
set(gca,'xtick',[-4*pi,-3*pi,-2*pi,-pi,0,pi,2*pi,3*pi,4*pi]);  
set(gca,'XTickLabel',...  
  {'-4pi','-3pi','-2pi','-pi','0','pi','2pi','3pi','4pi'})  
 
function s=speed(energy,angle)  
m=1; g=9.8; L=10;  
if angle<pi  
    s=sqrt(2/(m*L^2)*(energy-m*g*L*(1-cos(angle))));  
else  
    s=sqrt(2/(m*L^2)*(energy-m*g*L*(1+cos(angle-pi))));  
end  
 
function doPlots(data)  
plotCurves(data,0);  
plotCurves(data,2*pi);  
plotCurves(data,-2*pi);  
 
function plotCurves(data,k)  
plot(data(:,1)+k,data(:,2));  
hold on;  
plot(data(:,1)+k,-data(:,2));  
plot(-data(:,1)+k,-data(:,2));  
plot(-data(:,1)+k,data(:,2));

PIC

PIC

Maple

This solution below is due to Matrin Eisenberg

> restart;  
> plotEnergyLevels:= proc(m::positive, g::positive, L::positive, Erange::positive, numContours::posint)  
>   local plotOptions, maxTheta, thetaDot, plotContour, Emax, energies, contours;  
>   plotOptions := args[6..-1];  
>   maxTheta := E-> arccos(max(1-E/m/g/L, -1));  
>   thetaDot := E-> theta-> sqrt(max(2/m/L^2 * (E - m*g*L*(1-cos(theta))), 0));  
>   plotContour := E-> plot(thetaDot(E), 0..maxTheta(E), numpoints=5, plotOptions);  
>  
>   # Create first quadrant of middle region of the graph.  
>   Emax := Erange*m*g*L;  
>   energies := {problems/emax * n/numContours $ n=1..numContours};  
>   if Erange > 2 then energies := energies union {2*m*g*L}; fi;  
>   contours := map(plotContour, energies);  
>  
>   # Create other quadrants.  
>   map2(rcurry, plottools[reflect], {[[0,0],[0,1]], [0,0], [[0,0],[1,0]]});  
>   contours := contours union map(f-> map(f, contours)[], %);  
>  
>   # Create left and right regions of the graph.  
>   map2(rcurry, plottools[translate], {-2*Pi,2*Pi}, 0);  
>   contours := contours union map(f-> map(f, contours)[], %);  
>  
>   # Complete the graph.  
>   plots[display]([%[], plot([[-2*Pi,0], [0,0], [2*Pi,0]], plotOptions, style=point)],  
>     view=[-3*Pi..3*Pi, -thetaDot(Emax)(0)..thetaDot(Emax)(0)],  
>     title="Constant energy levels, undamped pendulum",  
>     labels=["angle", "angular speed"],  
>     xtickmarks=[seq](evalf(k*Pi)=sprintf("%a pi", k), k=-3..3));  
> end:  
>plotEnergyLevels(1, 9.8, 10, 5, 15);

PIC

78 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
> restart;   with(Statistics):  
> X := RandomVariable(Normal(mu,sigma));  
                                    X := _R  
 
> lik:=product(PDF(X,x[i]),i=1..n);  
                                /                             2 \  
                                |     1/2          (x[i] - mu)  |  
                         n      |    2    exp(-1/2 ------------)|  
                      --------' |                          2    |  
                     '  |  |    |                     sigma     |  
              lik :=    |  |    |1/2 ---------------------------|  
                        |  |    |              1/2              |  
                        |  |    \            Pi    sigma        /  
                       i = 1  
 
> lik:=expand(ln(lik)) assuming positive;  
                                                            2  
                                    1/2          (x[i] - mu)  
                        n          2    exp(-1/2 ------------)  
                      -----                              2  
                       \                            sigma  
               lik :=   )   ln(1/2 ---------------------------)  
                       /                     1/2  
                      -----                Pi    sigma  
                      i = 1  
 
>  
 
> eq1:=diff(lik,mu)=0;  
                                       /  n         \  
                                       |-----       |  
                               n mu    | \     x[i] |  
                     eq1 := - ------ + |  )   ------| = 0  
                                   2   | /         2|  
                              sigma    |----- sigma |  
                                       \i = 1       /  
 
> eq2:=diff(lik,sigma)=0;  
                                   /  n                         \  
                              2    |----- /                  2 \|  
                    n     n mu     | \    |  2 x[i] mu   x[i]  ||  
         eq2 := - ----- + ------ + |  )   |- --------- + ------|| = 0  
                  sigma        3   | /    |        3          3||  
                          sigma    |----- \   sigma      sigma /|  
                                   \i = 1                       /  
 
>  
> solve({problems/eq1,eq2},{mu,sigma});  
        n                        /  n       \2   /  n        \  
      -----                      |-----     |    |-----      |  
       \                         | \        |    | \        2|       2  
        )   x[i]          RootOf(|  )   x[i]|  - |  )   x[i] | n + _Z )  
       /                         | /        |    | /         |  
      -----                      |-----     |    |-----      |  
      i = 1                      \i = 1     /    \i = 1      /  
{mu = ----------, sigma = ---------------------------------------------},  
          n                                     n  
 
            n  
          -----  
           \  
            )   x[i]  
           /  
          -----  
          i = 1  
    {mu = ----------, sigma = 0}  
              n  
 
>

79 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
data = RandomReal[NormalDistribution[],1000];  
Histogram[data,30,"ProbabilityDensity",  
          ImageSize -> 300]]

PIC

Matlab
data = randn(1000,1);  
numberOfBins = 20;  
[freq,bin]   = hist(data,numberOfBins);  
binWidth     = (bin(end)-bin(1))/numberOfBins;  
currentArea  = sum(binWidth*freq);  
freq         = freq/currentArea;  
 
bar(bin,freq)

PIC

Maple
restart;  
> with(Statistics):  
> nSamples:=1000:  
> data := Sample(RandomVariable(Normal(0, 1)),nSamples):  
> plots[display](Histogram(data,bincount=20,  
>                          frequencyscale=relative));

PIC

80 Solve numerically the ODE u+u = f using the point collocation method

Problem: Give the ODE  

d4u x dx4 +u x = f

Solve numerically using the point collocation method. Use 5 points and 5 basis functions. Use the Boundary conditions u 0 = u 1 = u0 = u1 = 0

Use the trial function

g x = i=1N=5a i x x1i

Use f = 1

The solution is approximated using u x g x = i=1N=5a i x x1i.

N equations are generated for N unknowns (the unknowns being the undetermined coefficients of the basis functions). This is done by setting the error to zero at those points. The error being

d4g x dx4 +g xf

Once g x (the trial function is found) the analytical solution is used to compare with the numerical solution.

Mathematica
Clear["Global`*"];  
nBasis  = 5;  
nPoints = 5;  
a = Array[v,{nBasis}]  
 
Out[392]= {v[1],v[2],v[3],v[4],v[5]}  
 
trial    = Sum[a[[n]] (x(x-1))^n,{n,1,nBasis}];  
residual = D[trial,{x,4}]+trial-1;  
mat      = Flatten[Table[residual==0/.x->n/(2*nPoints-1),  
             {n,1,nPoints-1}]];  
 
mat = Join[mat,{2 a[[1]]+2a[[2]]==0}];  
sol = N@First@Solve[mat,a]  
 
Out[397]= {v[1.]->-0.0412493,  
           v[2.]->0.0412493,  
           v[3.]->0.000147289,  
           v[4.]->-0.0000245233,  
           v[5.]->-3.28259*10^-8}  
 
 
numericalSolution=Chop@FullSimplify@Sum[sol[[n,2]](x(x-1))^n,  
                   {n,1,nBasis}]  
 
numericalSolutionDiff4 = D[numericalSolution,{x,4}];  
numericalMoment = D[numericalSolution,{x,2}];  
 
Out[398]= x (0.0412493 +x^2 (-0.0826459+x (0.0416667 +  
               x (-0.00034374+x (-1.51283*10^-8+x (0.0000984213 +  
               x (-0.0000248515+  
               (1.64129*10^-7-3.28259*10^-8 x) x)))))))

now the analytical solution is obtained using DSolve and compared to the above numerical solution

eq = u''''[x]+u[x]==1;  
bc = {u[0]==0,u[1]==0,u''[0]==0,u''[1]==0};  
analyticalSolution = First@DSolve[{problems/eq,bc},u[x],x];  
analyticalSolution = u[x]/.analyticalSolution;  
analyticalSolution = FullSimplify[analyticalSolution]  
 
Out[346]= (2 (Cos[1/Sqrt[2]]-Cosh[1/Sqrt[2]])  
     (Cos[1/Sqrt[2]]+Cosh[1/Sqrt[2]]-  
     Cos[x/Sqrt[2]] Cosh[(-1+x)/Sqrt[2]]-  
     Cos[(-1+x)/Sqrt[2]] Cosh[x/Sqrt[2]]))/(Cos[Sqrt[2]]-Cosh[Sqrt[2]])  
 
analyticalMoment = D[analyticalSolution,{x,2}];  
dispError[pt_] := ((analyticalSolution-numericalSolution)/.x->pt)/  
                       (analyticalSolution/.x->.5)  
 
momentError[pt_]:=((analyticalMoment-numericalMoment)/.x->pt)/  
                          (analyticalMoment/.x->.5)

Now the percentage displacement and moment errors are plotted

Plot[dispError[z],{z,0,1},  
  ImagePadding->70,  
  ImageSize->400,  
  Frame->True,  
  GridLines->Automatic,  
  GridLinesStyle->Dashed,  
  PlotStyle->Red,  
  FrameLabel->{{"error",None},  
    {"x","percentage displacement error"}},  
    LabelStyle->14]

PIC

Plot[momentError[z],{z,0,1},  
  ImagePadding->70,  
  ImageSize->400,  
  Frame->True,  
  GridLines->Automatic,  
  GridLinesStyle->Dashed,  
  PlotStyle->Red,  
  FrameLabel->{{"error",None},  
    {"x","percentage moment error"}},  
    LabelStyle->14,PlotRange->All]

PIC

Now the Residual error distribution over x is plotted

Plot[Abs[numericalSolutionDiff4+numericalSolution-1],{x,0,1},  
  ImagePadding->80,  
  ImageSize->400,  
  Frame->True,  
  GridLines->Automatic,  
  GridLinesStyle->Dashed,  
  PlotStyle->Red,  
  FrameLabel->{{"error",None},  
   {"x","Residual error distribution over x"}},  
   LabelStyle->14,  
   PlotRange->All]

PIC

Maple
> #Solve u''''+u=1 using point collocation. Zero B.C.  
> #written by Nasser Abbasi for fun on April 25,2006.  
> restart;  
> Digits:=15;  
> nBasis  := 5:  
> nPoints := 5:  
> coef    := [seq(a[i],i=1..nBasis)]:  
> g := (x,n)->coef[n]*(x*(x-1))^n:  #basis function  
> f := x->sum(g(x,k),k=1..nBasis):  #trial function  
> moment := x->diff(f(x),x$2):  
> residual := x->diff(f(x),x$4)+f(x)-1:  
> A := seq(subs(x=N/(2*nPoints-1),residual(x)=0),N=1..nPoints-1):  
> A := A,2*(coef[1]+coef[2]=0):  
> sol   := solve([A],coef):  
> coef  := map(rhs,sol[1]):  
> evalf(coef):  
> numericalSolution:=x->sum(coef[n]*(x*(x-1))^n  ,n=1..nBasis):  
> `The numerical solution is `, numericalSolution(x);  
> numericalSolutionMoment := x->diff(numericalSolution(x),x$2):  
>  
> #Now obtain exact solution from Maple  
> eq    := diff(W(x),x$4)+W(x)=1:  
> bc    := W(0)=0,W(1)=0,D(D(W))(0)=0,D(D(W))(1)=0:  
>  
> exact := unapply(rhs(dsolve({problems/eq,bc})),x):  
>  
> exactMoment := x->diff(exact(x),x$2):  
>  
> displacmentError := x->(exact(x)-numericalSolution(x))/exact(.5):  
> momentError := x-> (exactMoment(x)-numericalSolutionMoment(x))/  
>                        evalf(subs(x=.5,exactMoment(x))):  
>  
> plot(displacmentError(x),x=0..1,labels=['x',"Disp. error %"]);  
> plot(momentError(x),x=0..1,labels=['x',"Moment error %"]);

PIC

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

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

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

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

PIC

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

PIC

82 Visualize a 2D matrix

Problem: Given a 2 dimensional matrix, say m×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,yindices as the coordinates.

Mathematica
mat = Table[Table[RandomReal[{0,100}],{20}],{10}];  
ListPointPlot3D[mat,Filling->Axis,ImageSize->300]

PIC

mat = HilbertMatrix[{20,10}];  
ListPlot3D[mat,  
     Mesh->None,  
     InterpolationOrder->0,  
     ColorFunction->"SouthwestColors",  
     Filling->Axis,  
     ImageSize->300]

PIC

Matlab
clear all; close all;  
 
A     = (100).*rand(20,20);  
[X,Y] = meshgrid(1:20,1:20);  
 
surfl(X,Y,A);  
shading interp  
colormap(gray);  
 
figure;  
A     = hilb(20);  
[X,Y] = meshgrid(1:20,1:20);  
surfl(X,Y,A);

PIC

PIC

Maple
restart;  
A := abs(LinearAlgebra[RandomMatrix](20,10,generator=0..100)):  
plots[matrixplot](A,heights=histogram);

PIC

A := LinearAlgebra[HilbertMatrix](20,10):  
A := convert(A,listlist):  
plots[listplot3d](A);

PIC

83 Find the particular and homogenous solution to undetermined system of equations

Problem: Find the general solution to Ax = b

2 4 6 4 2 5 7 6 2 3 5 2 x1 x2 x3 x3 = b1 b2 b3 where b1 b2 b3 = 4 3 5

In Maple 11, the LinearAlgebra package was used.

In Mathematica one can also get the general solution, but one must find the Null space specifically and add it to the result from LinearSolve[] since LinearSolve[] finds particular solutions only.

In Matlab the same thing needs to be done. I am not sure now how to make Matlab give me the same particular solution as Maple and Mathematica since Matlab Ab uses least square approach to determine a solution.  I am sure there is a way, will update this once I find out.

Mathematica
vars = {x1,x2,x3,x4};  
eqs  ={2 x1+4 x2+6 x3+4 x4==4,  
       2 x1+5 x2+7 x3+6 x4==3,  
       2 x1+3 x2+5 x3+2 x4==5};  
 
{b,mat} = CoefficientArrays[eqs,vars];  
Normal[mat]  
 
Out[162]= {{2,4,6,4},  
           {2,5,7,6},  
           {2,3,5,2}}  
 
Normal[b]  
 
Out[163]= {-4,-3,-5}  
 
 
(*Mathematica LinearSolve gives one  
  solution (particular solution)*)  
 
particularSolution = LinearSolve[mat,b]  
 
Out[164]= {-4,1,0,0}  
 
(*find the homogenous solution (nullspace) and add it  
   to the above particular solution*)  
 
homogenousSolution = (Transpose[NullSpace[mat]])  
 
Out[165]= {{2,-1},  
           {-2,-1},  
           {0,1},  
           {1,0}}  
 
(*Add the particular solution to the homogenous solution  
  to get the general solution*)  
 
fullSolution = particularSolution+homogenousSolution  
 
Out[166]= {{-2,-5},  
           {-1,0},  
           {0,1},  
           {1,0}}  
 
(*To obtain the general solution right away Solve  
  can be used instead*)  
 
sol = Flatten[Solve[eqs,vars]]  
 
During evaluation of In[167]:= Solve::svars: Equations may  
not give solutions for all "solve" variables. >>  
 
Out[167]= {x3->3/2-x1/2-x2/2,  
           x4->-(5/4)+x1/4-x2/4}  
 
generalSolution = vars/.sol  
 
Out[168]= {x1,  
           x2,  
           3/2-x1/2-x2/2,  
           -(5/4)+x1/4-x2/4  
           }  

Matlab
clear all;  
A=[2 4 6 4;  
   2 5 7 6;  
   2 3 5 2]  
 
b=[4 3 5]'  
particularSolution=A\b  
>>  
Warning: Rank deficient, rank = 2, tol =  4.033641e-15.  
 
particularSolution =  
 
         0  
         0  
    1.5000  
   -1.2500  
 
 
nullSolution=null(A,'r')  
>>  
nullSolution =  
 
    -1     2  
    -1    -2  
     1     0  
     0     1  

Maple
restart;  
with(LinearAlgebra);  
vars := x[1], x[2], x[3], x[4];  
sys := 2*x[1]+4*x[2]+6*x[3]+4*x[4] = 4,  
       2*x[1]+5*x[2]+7*x[3]+6*x[4] = 3,  
       2*x[1]+3*x[2]+5*x[3]+2*x[4] = 5;  
 
`eqs=`, convert([sys], Vector);  
A, b := GenerateMatrix([sys], [vars])

PIC

`general solution=`, LinearSolve(A, b, free = 'x')

PIC

Can solve this system to get the general solution using the solve command as follows

s := solve([sys], [vars]);  
r := subs(op(s), [vars]);  
`general solution=`, convert(%, Vector)

PIC

WARNING: Maple sometimes reorders the result from solve() so we can get a different ordering of the free variables as shown above.

84 Numerically integrate f(x) on the real line

Problem: Integrate

221 5 1 100 322+3x 98+x 37+x24 x 1+x2dx

The exact answer is 9425 = 3.76

Mathematica
f[x_] := (1/5)(1/100(322+3*x(98+x(37+x)))-24(x/(1+x^2)))  
r = Integrate[f[x],{x,-2,2}]  
 
Out[15]= 94/25  
 
N[r]  
 
Out[17]= 3.76

To compare with Matlab, replace 1 by 1.0 in the expression.

f[x_]:=(1.0/5)(1/100(322+3*x(98+x(37+x)))-24(x/(1+x^2)))  
r = Integrate[f[x],{x,-2,2}];  
InputForm[r]  
 
Out[62]=  3.7600000000000007  

Matlab
clear all;  
format long  
f=@(x) (1/5)*(1/100*(322+3*x.*(98+x.*(37+x)))-24*(x/(1+x.^2)));  
integral(f,-2,2)  
 
>>>  
ans =  
 
   3.760000000000001  
 
 
integral(f,-2,2,'AbsTol',1e-6,'RelTol',1e-6)  
>>>  
ans =  
 
   3.760000000000001  

85 Numerically integrate f(x,y) in 2D

Problem: Integrate

7101114 x2 +4ydxdy

The exact answer is 1719

Mathematica
f[x_,y_] := x^2+4 y  
r = Integrate[f[x,y],{x,11,14},{y,7,10}]  
 
Out[69]= 1719

Matlab
clear all;  
format long  
f=@(x,y) x.^2+4*y;  
integral2(f,11,14,7,10)  
 
ans =  
     1.719000000000004e+03  

86 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
data={1,2,7,9,3,4,10,12};  
filter=(1/3){1,1,1};  
data[[2;;-2]]=ListConvolve[filter,data];  
N[data]  
 
Out[66]= {1.,  
          3.33333,  
          6.,  
          6.33333,  
          5.33333,  
          5.66667,  
          8.66667,  
          12.}

Matlab
data   = [1 2 7 9 3 4 10 12];  
filter = (1/3)*[1 1 1];  
data(2:end-1) = conv(data,filter,'valid');  
data'  
 
ans =  
 
    1.0000  
    3.3333  
    6.0000  
    6.3333  
    5.3333  
    5.6667  
    8.6667  
   12.0000

87 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
data={{0,4,10,5,3},  
      {4,4,1,8,5},  
      {5,1,2,3,8},  
      {8,6,8,8,10},  
      {10,3,7,7,8}};  
 
filter= (1/4){{0,1,0},  
              {1,0,1},  
              {0,1,0}};  
 
data[[2;;-2,2;;-2]]=ListConvolve[filter,data];  
N[%]  
 
Out[58]= {{0., 4.,   10.,  5.,  3.},  
          {4., 2.5,  6.,   3.5, 5.},  
          {5., 4.25, 3.25, 6.5, 8.},  
          {8., 5.,   5.75, 7.,  10.},  
          {10.,3.,   7.,   7.,  8.}}  

Matlab
data=[0 4 10 5 3;  
      4,4,1,8,5;  
      5,1,2,3,8;  
      8,6,8,8,10;  
      10,3,7,7,8];  
 
filter= (1/4)*[0,1,0;  
              1,0,1;  
              0,1,0];  
 
data(2:end-1,2:end-1)=conv2(data,filter,'valid')  
 
data =  
         0    4.0000   10.0000    5.0000    3.0000  
    4.0000    2.5000    6.0000    3.5000    5.0000  
    5.0000    4.2500    3.2500    6.5000    8.0000  
    8.0000    5.0000    5.7500    7.0000   10.0000  
   10.0000    3.0000    7.0000    7.0000    8.0000

88 How to compute k=110 1 kk+1

Mathematica
In[3]:= N[Sum[1/(k*(k + 1)), {k, 10}]]  
 
Out[3]= 0.9090909090909091

Matlab
format long  
k=1:10;  
s=cumsum(1./(k.*(k+1)));  
s(end)  
 
 
ans =  
   0.909090909090909

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

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

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

PIC

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

PIC

90 How to make table of x,y values

Generate a list or table of (x,sin(x)) values and plot them.

Mathematica
lst = Table[{x, Sin[x]}, {x, -Pi, Pi, 1/10}];  
ListLinePlot[lst]  

PIC

Matlab

Matlab does not have a nice function like Table so one can either make the x variable and then use it to make sin(x) or use arrayfun

close all;  
x=-pi:.1:pi;  
plot(x,sin(x));

Using arrayfun

A=cell2mat(arrayfun(@(x) [x;sin(x)],-pi:.1:pi, 'UniformOutput',false));  
plot(A(1,:),A(2,:))

PIC

91 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

  v = {1, 2, 3, 4, 5, 6, 7, 8};
  f = {1/2, 1/2};
  ListConvolve[f, v] // N
  
  Out[34]= {1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}

For window size k = 3

  v = {1, 2, 3, 4, 5, 6, 7, 8};
  f = Table[1/3, {3}];
  ListConvolve[f, v] // N
  
  Out[40]= {2., 3., 4., 5., 6., 7.}

Matlab

For a window size k = 2

  V=[1 2 3 4 5 6 7 8];
  f=[1/2 1/2];
  conv(V,f,'valid')
  
  ans =
      1.5000    2.5000    3.5000    4.5000    5.5000    6.5000    7.5000

For window size k = 3

  V = [1 2 3 4 5 6 7 8];
  k = 3;
  f = ones(k,1)/k;
  conv(V,f,'valid')
  
  ans =
      2.0000    3.0000    4.0000    5.0000    6.0000    7.0000

92 How to select N random values from a set of numbers?

Given the set v1,2,3,5,6,7,11,12,20,21 how to select say m = 5 random numbers from it?

Mathematica
  method 1
  --------
  
  a = {1, 2, 3, 5, 6, 7, 11, 12, 20, 21};
  m = 5;
  b = Table[0, {m}];
  Do[k = RandomInteger[{1, Length[a]}];
    b[[i]] = a[[k]];
    a = Delete[a, k],
    {i, 1, m}
    ];
  b
  
  Out[47]= {6, 21, 3, 5, 11}
  
  method 2 (Version 9)
  --------------------
  RandomSample[a, m]
  (*  {1, 6, 11, 7, 20}  *)
  

Matlab
  A = [1,2,3,5,6,7,11,12,20,21];
  m = 5;
  B = zeros(m,1);
  for i = 1:m
      k = randi(length(A),1);
      B(i) = A(k);
      A(k) = [];
  end
  B
  ------------------
  
  B =
       2
      20
       7
      11
       1