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

April 9, 2014

### Contents

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 ﬁle 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 deﬁned by the transfer function

$H\left(s\right)=\frac{25}{{s}^{2}+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] 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]); Maple restart:  with(DynamicSystems):  sys := TransferFunction(25/(s^2+4*s+25)):  ResponsePlot(sys, Step(),duration=4);

### 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 deﬁned by the transfer function

$H\left(s\right)=\frac{1}{{s}^{2}+0.2s+1}$

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

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

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

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

Problem: Find the response for the continuous time system deﬁned by the transfer function

$H\left(s\right)=\frac{1}{{s}^{2}+0.2s+1}$

when the input is given by

$u\left(t\right)=sin\left(t\right)$

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

### 4 Obtain the poles and zeros of a transfer function

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

$H\left(s\right)=\frac{25}{{s}^{2}+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\left(s\right)$ 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 deﬁned by the transfer function

$H\left(s\right)=\frac{5s}{{s}^{2}+4s+25}$

 Mathematica Clear["Global*"];  sys = TransferFunctionModel[(5 s)/(s^2+4 s+25),s];  ss  = MinimalStateSpaceModel[StateSpaceModel[sys]] (a=Normal[ss][[1]])//MatrixForm  (b=Normal[ss][[2]])//MatrixForm  (c=Normal[ss][[3]])//MatrixForm  (d=Normal[ss][[4]])//MatrixForm 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 ﬁnd the step response

Problem: Find the state space representation and the step response of the continuous time system deﬁned 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 deﬁned by the transfer function

$H\left(s\right)=\frac{5s}{{s}^{2}+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}}   ] Matlab clear all;  s   = tf('s');  sys = 5*s / (s^2 + 4*s + 25);  bode(sys);  grid;  set(gcf,'Position',[10,10,400,400]);

### 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 deﬁned by

$H\left(s\right)=\frac{1+s}{{s}^{2}+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"] {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] 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

### 10 Convert transfer function to state space

Problem: Given the continuous time S transfer function deﬁned by

$H\left(s\right)=\frac{1+s}{{s}^{2}+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]] TransferFunctionModel[ss] 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 deﬁned by

$H\left(s\right)=\frac{{s}^{4}+8{s}^{3}+16{s}^{3}+9s+9}{{s}^{3}+6{s}^{2}+11s+6}$

obtain the partial-fractions decomposition.

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

 Mathematica 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 deﬁned in the following ﬁgure.

Function f(t) to obtain its Laplace transform

Comment: Mathematica solution was easier than Matlab's.

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

Not ﬁnding 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\left(s\right)=\frac{{s}^{4}+5{s}^{3}+6{s}^{2}+9s+30}{{s}^{4}+6{s}^{3}+21{s}^{2}+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[%]] 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\left(s\right)=\frac{{\omega }_{n}^{2}}{{s}^{2}+2\xi {\omega }_{n}s+{\omega }_{n}^{2}}$

in order to illustrate the response when the system is over, under, and critically damped. use ${\omega }_{n}=1$ and change $\xi$ over a range of values that extends from under damped to over damped.

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

### 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\left(s\right)=\frac{{\omega }_{n}^{2}}{{s}^{2}+2\xi {\omega }_{n}s+{\omega }_{n}^{2}}$

in order to illustrate the impulse response when the system is over, under, and critically damped. use ${\omega }_{n}=1$ and change $\xi$ over a range of values that extends from under damped to over damped.

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

### 16 Convert diﬀerential equation to transfer functions and state space

Problem: Obtain the transfer and state space representation for the diﬀerential equation

$3\frac{{d}^{2}y}{d{t}^{2}}+2\frac{dy}{dt}+y\left(t\right)=u\left(t\right)$

 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] tf = Simplify[TransferFunctionModel[ss]] 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\left(s\right)=\frac{{\omega }_{n}^{2}}{{s}^{2}+2\xi {\omega }_{n}s+{\omega }_{n}^{2}}$

Display the output and input on the same plot showing how the steady state error changes as the un damped natural frequency ${\omega }_{n}$ changes. Do this for ramp and step input.

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

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

Use maximum time of $10$ seconds and $\xi =0.707$ and change $omeg{a}_{n}$ from $0.2$ to $1.2$.

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

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

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

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

Problem: Convert

$H\left(s\right)=\frac{1}{{s}^{2}+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"] 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\left(z\right)=\frac{z}{z-1}$

ﬁnd $x\left[n\right]={F}^{-1}\left(z\right)$ 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}] 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

#### 19.2 problems/example 2

Problem: Given

$F\left(z\right)=\frac{5z}{{\left(z-1\right)}^{2}}$

ﬁnd $x\left[n\right]={F}^{-1}\left(z\right)$

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

### 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\left[n\right]=u\left[n\right]$ deﬁned as $x=\left\{1,1,1,\cdots \phantom{\rule{0.3em}{0ex}}\right\}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{1em}{0ex}}$ for $n\ge 0\phantom{\rule{0.3em}{0ex}}$, ﬁnd 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\left[n\right]={\left(\frac{1}{3}\right)}^{n}u\left(n\right)+{\left(0.9\right)}^{n-3}u\left(n\right)$

 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 ﬁnd and plot the plant output

system diagram

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

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"] 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}}] Show[{inputPlot,plantPlot},PlotLabel->"input/output plot"] 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]);

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

Problem: Given

$L\left(s\right)=\frac{s}{\left(s+4\right)\left(s+5\right)}$

as the open loop transfer function, how to ﬁnd $G\left(s\right)$, the closed loop transfer function for a unity feedback?

 Mathematica Clear["Global*"];  sys = TransferFunctionModel[s/((s+4)(s+5)),s];  closedLoop = SystemsModelFeedbackConnect[sys]; The system wrapper can be removed in order to obtain the rational polynomial expression as follows First[Flatten[closedLoop[[1,1]]/closedLoop[[1,2]]]] 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 ${x}_{1}\left[n\right]$ and ${x}_{2}\left[m\right]$, determine the linear convolution of the 2 sequences. Assume ${x}_{1}=\left[1,2,3,4,5\right]$ and ${x}_{2}=\left[3,4,5\right]$.

 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 ${x}_{1}\left[n\right]$ and ${x}_{2}\left[m\right]$, 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] {a,b}=JordanDecomposition[m];  MatrixForm[b] 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\left(2\pi {f}_{1}t\right)+4cos\left(2\pi {f}_{2}t\right)+5cos\left(2\pi {f}_{3}t\right)$

for the following frequency values (in Hz)

${f}_{1}=20,{f}_{2}=30,{f}_{3}=50$

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

 Mathematica 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] 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]   Here is the same plot as above, but using InterpolationOrder -> 0 And using InterpolationOrder -> 2 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'); 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]);

### 27 Solve the continuous-time algebraic Riccati equation

Problem: Solve for $X$ in the Riccati equation

${A}^{\prime }X+XA-XB{R}^{-1}{B}^{\prime }X+{C}^{\prime }C=0$

given

$\begin{array}{llll}\hfill A& =\left(\begin{array}{cc}\hfill -3\hfill & \hfill 2\hfill \\ \hfill 1\hfill & \hfill 1\hfill \end{array}\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill B& =\left(\begin{array}{c}\hfill 0\hfill \\ \hfill 1\hfill \end{array}\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill C& =\left(\begin{array}{cc}\hfill 1\hfill & \hfill -1\hfill \end{array}\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill R& =3\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$
 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]] 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

$\frac{1}{s\left(s+0.5\right)}$

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

The Riccati equation is given by

${A}^{\prime }X+XA-XB{R}^{-1}{B}^{\prime }X+{C}^{\prime }C=0$

Let $R=\left[3\right]$.

 Mathematica Clear["Global*"];  sys  = TransferFunctionModel[1/(s (s+.5)),s];  dsys = ToDiscreteTimeModel[sys,                             1,                             z,                             Method->"ZeroOrderHold"] ss = StateSpaceModel[dsys] 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[%] 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\left(z\right)$ along with the impulse response of its continuous time approximation $H\left(s\right)$

Plot the impulse response of $H\left(z\right)=z∕\left({z}^{2}-1.4z+0.5\right)$ and using sampling period of $T=0.5$ ﬁnd 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  ] (*Plot the impulse response of the continuous system*)  Plot[samplePeriod *continouseTimeResponse,{t,0,maxSimulationTime},               PlotRange->All,               Frame->True,               ImageSize->300] (*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"}}  ] (*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"}}  ] 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');

### 30 Find the system type given an open loop transfer function

Problem: Find the system type for the following transfer functions

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

To ﬁnd the system type, the transfer function is put in the form $\frac{k{\sum }_{i}\left(s-{s}_{i}\right)}{{s}^{M}{\sum }_{j}\left(s-{s}_{j}\right)}$, then the system type is the exponent $M$. Hence it can be seen that the ﬁrst system above has type one since the denominator can be written as ${s}^{1}\left(s-1\right)$ $\phantom{\rule{1em}{0ex}}$and the second system has type 2 since the denominator can be written as ${s}^{2}\left(s-1\right)$ and the third system has type 5. The following computation determines the type

 Mathematica 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

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

Find its eigenvalues and eigenvectors.

 Mathematica Remove["Global*"]  (a = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}) // MatrixForm N[Eigenvalues[a]]//MatrixForm N[Eigenvectors[a]]//MatrixForm 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

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

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 coeﬃcients 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

$\left(\begin{array}{cc}\hfill 1\hfill & \hfill 2\hfill \\ \hfill 3\hfill & \hfill 2\hfill \end{array}\right)$

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

 Mathematica 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 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=\left(\begin{array}{cc}\hfill 1\hfill & \hfill 2\hfill \\ \hfill 1\hfill & \hfill 3\hfill \end{array}\right)$

and

$b=\left(\begin{array}{c}\hfill 5\hfill \\ \hfill 8\hfill \end{array}\right)$

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

 Mathematica 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] 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 inﬁnity 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

$\begin{array}{llll}\hfill y& =x-1\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill y& =2x+1\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$

So

$A=\left(\begin{array}{cc}\hfill 1\hfill & \hfill 1\hfill \\ \hfill -2\hfill & \hfill 1\hfill \end{array}\right)$

and

$b=\left(\begin{array}{c}\hfill 1\hfill \\ \hfill -1\hfill \end{array}\right)$

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

$\begin{array}{ccc}\hfill 4x+4y+2z\hfill & \hfill =\hfill & \hfill 12\hfill \\ \hfill 5x+6y+3z\hfill & \hfill =\hfill & \hfill 7\hfill \\ \hfill 9x+7y+10z\hfill & \hfill =\hfill & \hfill 9\hfill \end{array}$

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 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 MatrixForm[RowReduce[mat]] 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 deﬁnite. This is done by checking that all the eigenvalues of the matrix A have negative real part.

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

Use

$sys=\frac{5s}{{s}^{2}+4s+25}$

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 coeﬃcient matrix

$A=\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill -2\hfill & \hfill -3\hfill \end{array}\right]$

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

 Mathematica Remove["Global*"];  (mat = {{0,1,0},{0,0,1},{-1,-2,-3}})//MatrixForm p = -LyapunovSolve[mat, IdentityMatrix[Length[mat]]];  MatrixForm[N[p]] 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.1\phantom{\rule{1em}{0ex}}sec$, generate the closed loop transfer function, and that poles of the closed loop transfer function are inside the unit circle

System block diagram.

 Mathematica Remove["Global*"]  plant  = TransferFunctionModel[s/(1.0+s),s];  plantd = ToDiscreteTimeModel[plant,                               1,                               z,                               Method->"ZeroOrderHold"] d = ToDiscreteTimeModel[        TransferFunctionModel[(s+1)/(s+5.0),s],             1,             z,             Method->"ZeroOrderHold"] sys=SystemsModelSeriesConnect[d,plantd] loopBack=SystemsModelFeedbackConnect[sys] 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] ss=StateSpaceModel[loopBack] 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]); [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 diﬀerence and check that all values in the resulting diﬀerence 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] mat2=Conjugate[Transpose[mat]];  (diff=mat-mat2)//MatrixForm Now check that each element in the diﬀerence 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)> 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}}}] 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}] 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]); ### 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\left(t\right)$ 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}}}] 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] 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]); ### 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}}}] 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] 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]); ### 45 Draw the root locus from the open loop transfer function Problem: Given $L\left(s\right)$, the open loop transfer function, draw the root locus. Let $L\left(s\right)=\frac{{s}^{2}+2s+4}{s\left(s+4\right)\left(s+6\right)\left({s}^{2}+1.4s+1\right)}$ Root locus is the locus of the closed loop dominant pole as the gain $k$ is varied from zero to inﬁnity.  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] 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]); ### 46 Find the cross correlation between two sequences Problem: Given $\begin{array}{llll}\hfill A& =\left[0,0,2,-1,3,7,1,2,-3,0,0\right]\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill B& =\left[0,0,1,-1,2,-2,4,1,-2,5,0,0\right]\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$ 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 diﬀerent norms of a vector Problem: Given the vector say $v=1,2,4$ Find its norm for $p=1,2,\infty$  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, ﬁnd the set of orthonormal vectors that span the same space as $A$ and verify the result. Let $A=\left[\begin{array}{cccc}\hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 2\hfill \\ \hfill 1\hfill & \hfill 2\hfill & \hfill 3\hfill & \hfill 4\hfill \\ \hfill 2\hfill & \hfill 0\hfill & \hfill 2\hfill & \hfill 0\hfill \end{array}\right]$ 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 ﬁrst 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\left(y,x\right)$ Problem: Plot the function $f\left(x,y\right)={x}^{2}-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] 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]); ### 50 Find the point of intersection of 3 surfaces Problem: Given 3 surfaces ﬁnd the point where they intersect by solving the linear system and display the solution. Let the 3 surfaces be $\begin{array}{llll}\hfill z& =2x+y\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill z& =2y\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill z& =2\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$ In Mathematica, plot the surfaces using the Plot3D command, and then use LinearSolve to solve the system of equations in order to ﬁnd 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] {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)"]  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 ### 51 Obtain the LU decomposition of a matrix Problem: Given a matrix $A$, ﬁnd 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. LU decomposition.  Mathematica Remove["Global*"] mat = {{1, -3, 5}, {2, -4, 7}, {-1, -2, 1}}; MatrixForm[mat] {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_}/;j1,{3,3}]+ IdentityMatrix[3] Out[]= {{1,0,0},{2,1,0},{-1,-(5/2),1}} MatrixForm[lower] upper=lu SparseArray[{i_,j_}/;j>=i->1,{3,3}] Out[] = {{1,-3,5},{0,2,-3},{0,0,-(3/2)}} MatrixForm[upper] lower.upper Out[8]= {{1,-3,5},{2,-4,7},{-1,-2,1}} MatrixForm[%] 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 ﬁle Problem: download a wav ﬁle and display the frequency spectrum of the audio signal using FFT. The Matlab example was based on Matheworks tech note 1702.  Mathematica Remove["Global*"]; fname = "stereo.wav"; SetDirectory[ToFileName[Extract ["FileName"/.NotebookInformation[ EvaluationNotebook[]],{1},FrontEndFileName]]]; ele = Import[fname,"Elements"] 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}] 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)|'); ### 53 Solve homogeneous 1st order linear diﬀerential equation with constant coeﬃcients and initial conditions Problem: Solve ${y}^{\prime }\left(t\right)=3y\left(t\right)$ with initial condition $y\left(0\right)=1$ and plot the solution for $t=0\cdots 2$ seconds.  Mathematica Remove["Global*"]; sol = First@DSolve[{y'[t]==3y[t],y[0]==1},y[t],t]; y = y[t]/.sol Out[26]= E^(3 t) Plot[y,{t,0,2}, FrameLabel->{{"y(t)",None}, {"t","Solution of y'=3y, y(0)=1"}}, Frame->True, GridLines->Automatic, GridLinesStyle->Automatic, RotateLabel->False, ImageSize->300, AspectRatio->1] 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 ### 54 Solve homogeneous 2nd order linear diﬀerential equation with constant coeﬃcients and initial conditions Problem: Solve ${y}^{\prime \prime }\left(t\right)-1.5{y}^{\prime }\left(t\right)+5y\left(t\right)=0$ with initial conditions $y\left(0\right)=1,{y}^{\prime }\left(0\right)=0$ To use Matlab ode45, the second order ODE is ﬁrst converted to state space formulation as follows Given ${y}^{\prime \prime }\left(t\right)-1.5{y}^{\prime }\left(t\right)+5y\left(t\right)=0$ let $\begin{array}{llll}\hfill {x}_{1}& =y\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill {x}_{2}& ={y}^{\prime }\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & ={x}_{1}^{\prime }\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$ hence ${x}_{1}^{\prime }={x}_{2}$ and $\begin{array}{llll}\hfill {x}_{2}^{\prime }& ={y}^{\prime \prime }\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =1.5{y}^{\prime }-5y\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =1.5{x}_{2}-5{x}_{1}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$ Hence we can now write $\left[\begin{array}{c}\hfill {x}_{1}^{\prime }\hfill \\ \hfill {x}_{2}^{\prime }\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -5\hfill & \hfill 1.5\hfill \end{array}\right]\left[\begin{array}{c}\hfill {x}_{1}\hfill \\ \hfill {x}_{2}\hfill \end{array}\right]$ 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}] 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 ### 55 Solve non-homogeneous 2nd order linear diﬀerential equation with constant coeﬃcients and initial conditions Problem: Solve ${y}^{\prime \prime }\left(t\right)-1.5{y}^{\prime }\left(t\right)+5y\left(t\right)=4cos\left(t\right)$ with initial conditions $y\left(0\right)=1,{y}^{\prime }\left(0\right)=0$ To use Matlab ode45, the second order ODE is converted to state space as follows Given ${y}^{\prime \prime }\left(t\right)-1.5{y}^{\prime }\left(t\right)+5y\left(t\right)=4cos\left(t\right)$, let $\begin{array}{llll}\hfill {x}_{1}& =y\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill {x}_{2}& ={y}^{\prime }\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & ={x}_{1}^{\prime }\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$ hence ${x}_{1}^{\prime }={x}_{2}$ and $\begin{array}{llll}\hfill {x}_{2}^{\prime }& ={y}^{\prime \prime }\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =1.5{y}^{\prime }-5y+4cos\left(t\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =1.5{x}_{2}-5{x}_{1}+4cos\left(t\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$ Hence we can now write $\left[\begin{array}{c}\hfill {x}_{1}^{\prime }\hfill \\ \hfill {x}_{2}^{\prime }\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -5\hfill & \hfill 1.5\hfill \end{array}\right]\left[\begin{array}{c}\hfill {x}_{1}\hfill \\ \hfill {x}_{2}\hfill \end{array}\right]+\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 4\hfill \end{array}\right]cos\left(t\right)$  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}] 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 ### 56 Solve homogeneous 2nd order linear diﬀerential equation with constant coeﬃcients with boundary values (BVP) Problem: Solve ${y}^{\prime \prime }\left(t\right)+t\phantom{\rule{1em}{0ex}}y\left(t\right)=0$ with the boundary conditions $y\left(0\right)=3,y\left(20\right)=-1$ For solving with Matlab, the ODE is ﬁrst converted to state space as follows Given ${y}^{\prime \prime }\left(t\right)+t\phantom{\rule{1em}{0ex}}y\left(t\right)=0$, let $\begin{array}{llll}\hfill {x}_{1}& =y\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill {x}_{2}& ={y}^{\prime }\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & ={x}_{1}^{\prime }\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$ Therefore ${x}_{1}^{\prime }={x}_{2}$ And $\begin{array}{llll}\hfill {x}_{2}^{\prime }& ={y}^{\prime \prime }\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =-t\phantom{\rule{1em}{0ex}}y\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =-t\phantom{\rule{1em}{0ex}}{x}_{1}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$ This results in $\left[\begin{array}{c}\hfill {x}_{1}^{\prime }\hfill \\ \hfill {x}_{2}^{\prime }\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -t\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{c}\hfill {x}_{1}\hfill \\ \hfill {x}_{2}\hfill \end{array}\right]$ 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] 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]; ### 57 Solve the 1-D heat partial diﬀerential equation (PDE) The PDE is $\frac{\partial T\left(x,t\right)}{\partial t}=k\frac{{\partial }^{2}T\left(x,t\right)}{\partial {x}^{2}}$ Problem: given a bar of length $L$ and initial conditions $T\left(x,0\right)=sin\left(\pi x\right)$ and boundary conditions $T\left(0,t\right)=0,T\left(L,t\right)=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] ContourPlot[y,{x,0,barLength},{t,0,timeDuration}, PlotPoints->15, Contours->15, ContourLines->False, ColorFunction->Hue, PlotRange->All, ImageSize->300] 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; ### 58 Show the eﬀect of diﬀerent boundary/initial conditions on the solution of the 1-D heat partial diﬀerential equation (PDE) The PDE is $\frac{\partial T\left(x,t\right)}{\partial t}=k\frac{{\partial }^{2}T\left(x,t\right)}{\partial {x}^{2}}$ Problem: given a bar of length $L$ , solve the above 1-D heat PDE for 4 diﬀerent 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. ### 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

#### 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\left(t\right)$

 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}}] 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 ${\nabla }^{2}T\left(x,y\right)=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 ﬁxed temperature as shown, ﬁnd the steady state temperature distribution

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

### 62 Generate and plot one pulse signal of diﬀerent width and amplitude

Problem: Generate one signal of some width and height.

 Mathematica Remove["Global*"]  f[x_,w_,h_]:=If[0<=x{{-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] 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

### 63 Generate and plot train of pulses of diﬀerent 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] 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

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

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

 Mathematica Chop[Fourier[{1, 2, 3},               FourierParameters ->{1, -1}]]    Out[5]= {6.,           -1.5 + 0.8660254037844386*I,           -1.5 - 0.8660254037844386*I} Matlab fft([1,2,3])'    ans =      6.0000     -1.5000 - 0.8660i     -1.5000 + 0.8660i
 Ada I posted the code below on comp.lang.ada on June 8,2010. Thanks to Ludovic Brenta for making improvement to the Ada code. 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;

### 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 diﬀerentiation 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]] 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)
 Maple restart;  > with(Statistics):  > nSamples:=1000:  > data := Sample(RandomVariable(Normal(0, 1)),nSamples):  > plots[display](Histogram(data,bincount=20,  >                          frequencyscale=relative));

### 80 Solve numerically the ODE ${u}^{\prime \prime \prime \prime }+u=f$ using the point collocation method

Problem: Give the ODE

$\frac{{d}^{4}u\left(x\right)}{d{x}^{4}}+u\left(x\right)=f$

Solve numerically using the point collocation method. Use 5 points and 5 basis functions. Use the Boundary conditions $u\left(0\right)=u\left(1\right)={u}^{\prime \prime }\left(0\right)={u}^{\prime \prime }\left(1\right)=0$

Use the trial function

$g\left(x\right)=\sum _{i=1}^{N=5}{a}_{i}{\left(x\left(x-1\right)\right)}^{i}$

Use $f=1$

The solution is approximated using $u\left(x\right)\approx g\left(x\right)=\sum _{i=1}^{N=5}{a}_{i}{\left(x\left(x-1\right)\right)}^{i}$.

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

$\frac{{d}^{4}g\left(x\right)}{d{x}^{4}}+g\left(x\right)-f$

Once $g\left(x\right)$ (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] 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] 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]
 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 %"]);

### 81 Plot the dynamic response factor ${R}_{d}$ of a system as a function of $r=\frac{\omega }{{\omega }_{n}}$ for diﬀerent damping ratios

Problem: Plot the standard curves showing how the dynamic response ${R}_{d}$ changes as $r=\frac{\omega }{{\omega }_{n}}$ changes. Do this for diﬀerent damping ratio $\xi$. Also plot the phase angle.

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

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

### 82 Visualize a 2D matrix

Problem: Given a 2 dimensional matrix, say $m×n\phantom{\rule{0.3em}{0ex}}$, how to visualize its content?

These are some examples showing how to visualize a matrix as a 3D data, where the height is taken as the values of the matrix entries, and the $x,y\phantom{\rule{0.3em}{0ex}}\phantom{\rule{1em}{0ex}}$indices as the coordinates.

 Mathematica mat = Table[Table[RandomReal[{0,100}],{20}],{10}];  ListPointPlot3D[mat,Filling->Axis,ImageSize->300] mat = HilbertMatrix[{20,10}];  ListPlot3D[mat,       Mesh->None,       InterpolationOrder->0,       ColorFunction->"SouthwestColors",       Filling->Axis,       ImageSize->300] 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);
 Maple restart;  A := abs(LinearAlgebra[RandomMatrix](20,10,generator=0..100)):  plots[matrixplot](A,heights=histogram); A := LinearAlgebra[HilbertMatrix](20,10):  A := convert(A,listlist):  plots[listplot3d](A);

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

Problem: Find the general solution to $Ax=b$

$\left[\begin{array}{cccc}\hfill 2\hfill & \hfill 4\hfill & \hfill 6\hfill & \hfill 4\hfill \\ \hfill 2\hfill & \hfill 5\hfill & \hfill 7\hfill & \hfill 6\hfill \\ \hfill 2\hfill & \hfill 3\hfill & \hfill 5\hfill & \hfill 2\hfill \end{array}\right]\left[\begin{array}{c}\hfill {x}_{1}\hfill \\ \hfill {x}_{2}\hfill \\ \hfill {x}_{3}\hfill \\ \hfill {x}_{3}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill {b}_{1}\hfill \\ \hfill {b}_{2}\hfill \\ \hfill {b}_{3}\hfill \end{array}\right]\phantom{\rule{1em}{0ex}}where\left[\begin{array}{c}\hfill {b}_{1}\hfill \\ \hfill {b}_{2}\hfill \\ \hfill {b}_{3}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill 4\hfill \\ \hfill 3\hfill \\ \hfill 5\hfill \end{array}\right]$

In Maple 11, the LinearAlgebra package was used.

In Mathematica one can also get the general solution, but one must ﬁnd the Null space speciﬁcally and add it to the result from LinearSolve[] since LinearSolve[] ﬁnds 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 A$\setminus$b uses least square approach to determine a solution.  I am sure there is a way, will update this once I ﬁnd 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]) general solution=, LinearSolve(A, b, free = 'x') 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) WARNING: Maple sometimes reorders the result from solve() so we can get a diﬀerent ordering of the free variables as shown above.

### 84 Numerically integrate $f\left(x\right)$ on the real line

Problem: Integrate

${\int }_{-2}^{2}\frac{1}{5}\left(\frac{1}{100}\left(322+3x\left(98+x\left(37+x\right)\right)\right)-24\frac{x}{1+{x}^{2}}\right)dx$

The exact answer is $94∕25=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\left(x,y\right)$ in 2D

Problem: Integrate

${\int }_{7}^{10}{\int }_{11}^{14}\left({x}^{2}+4y\right)dxdy$

 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 ﬁlter on 1D numerical data (a vector)

Problem: suppose we want to apply say an averaging ﬁlter 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 ﬁlter).

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 ﬁlter on 2D numerical data (a matrix)

Problem: Apply a Laplacian ﬁlter 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 $\sum _{k=1}^{10}\frac{1}{k\left(k+1\right)}$

 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 ﬁnd closed loop step response to a plant with a PID controller?

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

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

### 90 How to make table of $x,y$ values

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

 Mathematica lst = Table[{x, Sin[x]}, {x, -Pi, Pi, 1/10}];  ListLinePlot[lst] 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\left(x\right)$ 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,:))

### 91 How to ﬁnd the moving average of a 1D sequence?

Given some sequence such as $1,2,3,4,5,6,7$ how to ﬁnd the moving average for diﬀerent 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