home
PDF letter size
PDF legal size

General report

How to solve Basic Engineering and Mathematics Problems

using Mathematica, Matlab and Maple

November 30, 2019   Compiled on November 30, 2019 at 2:43pm

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

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

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

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

### 1.1 Creating tf and state space and diﬀerent Conversion of forms

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} $\text{tf} = \left [{\begin{array}{c}{\frac{5\,{s}^{2}+15\,s+10}{{s}^{3}+10\,{s}^{2}+24\,s}}\end{array}} \right ]$ .syntaxhighlighter textarea {font-size: 14px !important;} $5*s^2+15*s+10$ .syntaxhighlighter textarea {font-size: 14px !important;} $s*(s^2+10*s+24)$

#### 1.1.2 Convert transfer function to state space representation

1.1.2.1 problem 1
1.1.2.2 problem 2

##### 1.1.2.1 problem 1

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} $\left [{\begin{array}{cc} 0 & 1 \\ -25 & -4 \end{array}} \right ]$ .syntaxhighlighter textarea {font-size: 14px !important;} $\left [{\begin{array}{c} 0\\ 1\end{array}} \right ]$ .syntaxhighlighter textarea {font-size: 14px !important;} $\left [{\begin{array}{cc} 0&5\end{array}} \right ]$ .syntaxhighlighter textarea {font-size: 14px !important;} $\left [{\begin{array}{c} 0\end{array}} \right ]$

##### 1.1.2.2 problem 2

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cc|c|c} 0 & 1 & 0 \\ -1 & -1 & 1 \\ \hline 1 & 1 & 0 \\ \end{array}} \right )_{}$ .syntaxhighlighter textarea {font-size: 14px !important;} $\frac{s+1}{s^2+s+1}$

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

#### 1.1.3 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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} $\frac{1}{3{s}^{2}+2\,s+1}$ .syntaxhighlighter textarea {font-size: 14px !important;} \left \{ \left [{\begin{array}{cc} 0&1\\ \noalign{\medskip }-1/3&-2/3 \end{array}} \right ] , \left [{\begin{array}{c} 0\\ \noalign{\medskip } 1\end{array}} \right ] , \left [{\begin{array}{cc} 1/3&0\end{array}} \right ] , \left [{\begin{array}{c} 0\end{array}} \right ] \right \}

#### 1.1.6 Convert from continuous transfer function to discrete time transfer function

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

#### 1.1.7 Convert a Laplace transfer function to an ordinary diﬀerential equation

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.2 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(s) = \frac{25}{s^{2}+4s+25}$

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;}

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

Problem: Find the impulse and step responses for the continuous time system 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 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple Using Maple DynamicSystems .syntaxhighlighter textarea {font-size: 14px !important;} Using Laplace transform method: .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.4 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(s)=\frac{1}{s^{2}+0.2s+1}$ when the input is given by $u(t)=\sin (t)$ and display the response and input on the same plot.

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} Using DynamicSystem package .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.5 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(s)=\frac{25}{s^{2}+4s+25}$

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.6 Generate Bode plot of a transfer function

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} or can plot the the two bode ﬁgures on top of each others as follows .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.7 How to check that state space system $$x'=Ax+Bu$$ is controllable?

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cccc} 0 & 1 & 0 & 2 \\ 1 & 0 & 2 & 0 \\ 0 & -2 & 0 & -10 \\ -2 & 0 & -10 & 0 \\ \end{array}} \right )$ .syntaxhighlighter textarea {font-size: 14px !important;} True .syntaxhighlighter textarea {font-size: 14px !important;} 4

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} 4

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} \left [{\begin{array}{cccc} 0&1&0&2\\ \noalign{\medskip }1&0&2&0 \\ \noalign{\medskip }0&-2&0&-10\\ \noalign{\medskip }-2&0&-10&0 \end{array}} \right ] .syntaxhighlighter textarea {font-size: 14px !important;} true .syntaxhighlighter textarea {font-size: 14px !important;} true .syntaxhighlighter textarea {font-size: 14px !important;} 4

### 1.8 Obtain partial-fraction expansion

Problem: Given the continuous time S transfer function deﬁned by $H(s)=\frac{s^{4}+8s^{3}+16s^{2}+9s+9}{s^{3}+6s^{2}+11s+6}$ obtain the partial-fractions decomposition.

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} $s+2- \frac{7}{s+2}-{\frac{9}{2\,s+6}}+{\frac{9}{2\,s+2 }}$ .syntaxhighlighter textarea {font-size: 14px !important;} $[s,2,\frac{7}{s+2},-{\frac{9}{2\,s+6}},{\frac{9}{2\,s+2}}]$

### 1.9 Obtain Laplace transform for a piecewise functions

Problem: Obtain the Laplace transform for the function 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 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple With Maple, had to use Heaviside else Laplace will not obtain the transform of a piecewise function. .syntaxhighlighter textarea {font-size: 14px !important;} ${\frac{1-{{e}^{-sT}}}{{s}^{2}}}$

### 1.10 Obtain Inverse Laplace transform of a transfer function

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

Mathematica

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

Matlab

Maple

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

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

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

Mathematica

Matlab

Maple

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

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

1.12.1 Mathematica

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

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

The steady state error is the 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 $$\omega _{n}$$ from $$0.2$$ to $$1.2$$.

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

#### 1.12.1 Mathematica

1.12.1.1 ramp input
1.12.1.2 step input

### 1.13 Show the use of the inverse Z transform

1.13.1 example 1
1.13.2 example 2

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

#### 1.13.1 example 1

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

#### 1.13.2 example 2

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} Inverse Z is 5 n .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.14 Find the Z transform of sequence x(n)

1.14.1 example 1
1.14.2 example 2

#### 1.14.1 example 1

Find the Z transform for the unit step discrete function

Given the unit step function $$x[n]=u[n]$$ deﬁned as $$x=\{1,1,1,\cdots \}\,\$$ for $$n\geq 0\,$$, ﬁnd its Z transform.

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

#### 1.14.2 example 2

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.15 Sample a continuous time system

Given the following system, sample the input and ﬁnd and plot the plant output

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

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

Mathematica

Matlab

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} The system wrapper can be removed in order to obtain the rational polynomial expression as follows .syntaxhighlighter textarea {font-size: 14px !important;} $\frac{s}{s^2+10 s+20}$

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cccccc} 3 & -1 & 1 & 1 & 0 & 0 \\ 1 & 1 & -1 & -1 & 0 & 0 \\ 0 & 0 & 2 & 0 & 1 & 1 \\ 0 & 0 & 0 & 2 & -1 & -1 \\ 0 & 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 & 1 \\ \end{array}} \right )$ .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cccccc} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & 1 & 0 & 0 & 0 \\ 0 & 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 2 \\ \end{array}} \right )$

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cc} 0.589517 & 1.82157 \\ 1.82157 & 8.81884 \\ \end{array}} \right )$

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.19 Solve the discrete-time algebraic Riccati equation

Problem: Given a continuous-time system represented by a transfer function $\frac{1}{s(s+0.5)}$ convert this representation to state space and sample the system at sampling period of $$1$$ second, and then solve the discrete-time Riccati equation.

The Riccati equation is given by $A^{\prime }X+XA-XBR^{-1}B^{\prime }X+C^{\prime }C=0$

Let $$R=[3]$$.

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cc} 0.671414 & -0.977632 \\ -0.977632 & 2.88699 \\ \end{array}} \right )$

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

Plot the impulse response of $$H(z)=z/(z^2-1.4 z+0.5)$$ and using sampling period of $$T=0.5$$ ﬁnd continuous time approximation using zero order hold and show the impulse response of the system and compare both responses.

Mathematica

Find its impulse response

approximate to continuous time, use ZeroOrderHold

Find the impulse response of the continuous time system

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

Plot the impulse response of the discrete system

Plot the impulse response of the continuous system

Plot both responses on the same plot

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

Matlab

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

Problem: Find the system type for the following transfer functions

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

To ﬁ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 )$$ 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 .syntaxhighlighter textarea {font-size: 14px !important;} Out[171]= 1 .syntaxhighlighter textarea {font-size: 14px !important;} Out[173]= 2 .syntaxhighlighter textarea {font-size: 14px !important;} Out[175]= 5

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.22 Find the eigenvalues and eigenvectors of a matrix

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ \end{array}} \right )$ .syntaxhighlighter textarea {font-size: 14px !important;} $\left \{\frac{3}{2} \left (5+\sqrt{33}\right ),\frac{3}{2} \left (5-\sqrt{33}\right ),0\right \}$ $\{16.1168,-1.11684,0.\}$ .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{ccc} -\frac{-15-\sqrt{33}}{33+7 \sqrt{33}} & \frac{4 \left (6+\sqrt{33}\right )}{33+7 \sqrt{33}} & 1 \\ -\frac{15-\sqrt{33}}{-33+7 \sqrt{33}} & \frac{4 \left (-6+\sqrt{33}\right )}{-33+7 \sqrt{33}} & 1 \\ 1 & -2 & 1 \\ \end{array}} \right )$ $\left ({\begin{array}{ccc} 0.283349 & 0.641675 & 1. \\ -1.28335 & -0.141675 & 1. \\ 1. & -2. & 1. \\ \end{array}} \right )$

 Matlab Matlab generated eigenvectors are such that the sum of the squares of the eigenvector elements add to one. .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.23 Find the characteristic polynomial of a matrix

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $-x^3+6 x^2+72 x+27$

 Matlab Note: Matlab generated characteristic polynomial coeﬃcients are negative to what Mathematica generated. .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

Problem, given the matrix $\left ({\begin{array} [c]{cc}1 & 2\\ 3 & 2 \end{array}} \right )$ Verify that matrix is a zero of its characteristic polynomial. The Characteristic polynomial of the matrix is found, then evaluated for the matrix. The result should be the zero matrix.

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $$x^2-3 x-4$$ .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cc} 0 & 0 \\ 0 & 0 \\ \end{array}} \right )$ Another way is as follows .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cc} 0 & 0 \\ 0 & 0 \\ \end{array}} \right )$

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

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.25 How to check for stability of 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}$

#### 1.25.1 Checking stability using transfer function poles

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

#### 1.25.2 Checking stability using state space A matrix

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.26 Check continuous system stability in the Lyapunov sense

Problem: Check the stability (in Lyapunov sense) for the state coeﬃcient matrix $A=\begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ -1 & -2 & -3 \end{bmatrix}$

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

We must transpose the matrix $$A$$ when calling these functions, since the Lyapunov equation is deﬁned as $$A^TP+PA=-Q$$ and this is not how the software above deﬁnes them. By simply transposing the $$A$$ matrix when calling them, then the result will be correct.

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ -1 & -2 & -3 \\ \end{array}} \right )$ .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{ccc} 2.3 & 2.1 & 0.5 \\ 2.1 & 4.6 & 1.3 \\ 0.5 & 1.3 & 0.6 \\ \end{array}} \right )$ .syntaxhighlighter textarea {font-size: 14px !important;} $$\{6.18272,1.1149,0.202375\}$$

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} \left [{\begin{array}{c} 6.18272045921436+ 0.0\,i \\ \noalign{\medskip } 1.11490451203192+ 0.0\,i\\ \noalign{\medskip } 0.202375028753723+ 0.0\,i\end{array}} \right ]

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

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

System block diagram.

Mathematica

Now generate unit step response

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

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

Poles are inside the unit circle, hence stable.

Matlab

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.31 Draw the root locus from the open loop transfer function

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

Root locus is the locus of the closed loop dominant pole as the gain $$k$$ is varied from zero to inﬁnity.

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.32 Find $$e^{A t}$$ where $$A$$ is a matrix

Mathematica

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

Now verify the result by solving for $$e^{At}$$ using the method would one would do by hand, if a computer was not around. There are a number of methods to do this by hand. The eigenvalue method, based on the Cayley Hamilton theorem will be used here. Find the eigenvalues of $$|A-\lambda I|$$

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

$\lambda ^2+3 \lambda +2$

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

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

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

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

Matlab

### 1.33 Draw root locus for a discrete system

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

Mathematica

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

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

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

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

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

Here is example of obtaining the analytical expression of the response

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

Matlab

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

#### 1.35.1 example 1, single input, single output closed loop

Given the following simple closed loop system, show the step response. Let mass $$M=1\text{kg}$$, damping coeﬃcient $$c=1 \text{newton-seconds per meter}$$ and let the stiﬀness coeﬃcient be $$k=20\text{newton per meter}$$.

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

Mathematica

Matlab

Another way to do the above is

### 1.36 Compare the eﬀect on the step response of a standard second order system as $$\zeta$$ changes

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

Change $$\zeta$$ and plot the step response for each to see the eﬀect of changing $$\zeta$$ (the damping coeﬃcient).

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

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

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

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

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

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

Mathematica

Matlab

### 1.37 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 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.38 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 .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

### 1.39 How to make Nyquist plot?

1.39.1 Example 1
1.39.2 Example 2

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

Here are two examples.

#### 1.39.1 Example 1

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

#### 1.39.2 Example 2

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

 .syntaxhighlighter textarea {font-size: 14px !important;}

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

Mathematica and Matlab allow one to do pretty much the same operations in the area of linear algebra and matrix manipulation. Two things to keep in mind is that Mathematica uses a more general way to store data.

Mathematica uses ragged arrays or a list of lists. This means rows can have diﬀerent sizes. (these are the lists inside the list). So a Mathematica matrix is stored in a list of lists. This is similar in a way to Matlab cell data structure, since each raw can have diﬀerent length. In a standard matrix each row must have the same length.

In Matlab one can also have ragged arrays, these are called cells. In Mathematica, there is one data structure for both.

Another thing to keep in mind is that Matlab, due to its Fortran background is column major when it comes to operations on matrices. This simple example illustrate this diﬀerence. Suppose we have a matrix $$A$$ of 3 rows, and want to ﬁnd the location where $$A(i,j)=2$$ where $$i$$ is the row number and $$j$$ is the column number. Given this matrix $A = \left ({\begin{array}{cccc} 1 & 2 & 3 & 4 \\ 2 & 3 & 1 & 5 \\ 5 & 6 & 7 & 2 \end{array}} \right )$ Then the result of using the find() command in Matlab is

The Matlab result gives the order of the rows it found the element at based on searching column wise since it lists the second row ﬁrst in its result. Compare this to Mathematica Position[] command

which gives

Mathematica searched row-wise.

Mathematica use for matrix manipulate takes more time to master compared to Matlab, since Mathematica data structure is more general and little more complex (ragged arrays) compared to Matlab’s since Mathematica also has to support symbolics in its commands and not just numbers

In Maple the following short cuts can be used enter vectors and matrices: For row vector:  v:=<1|2|3|4> and for column vector  v:=<1,2,3,4> and for matrix of say 2 rows and 3 columns  A:=<1,2|3,4|5,6> so — acts as column separator. There are other ways to do this (as typical in Maple), but I ﬁnd the above the least confusing. For transpose  A^%T can be used.

### 2.1 Multiply matrix with a vector

How to perform the following matrix/vector multiplication?

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

compile and run

Fortran

compile and run

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} $\left [{\begin{array}{c} 14\\ 32\\ 50 \end{array}} \right ]$

 Python .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.2 Insert a number at speciﬁc position in a vector or list

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Fortran .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} Using <> notation .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Python Python uses zero index. .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.3 Insert a row into a matrix

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple Using <<>> notation .syntaxhighlighter textarea {font-size: 14px !important;} Using Matrix/Vector .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Python .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Fortran .syntaxhighlighter textarea {font-size: 14px !important;} Compile and run .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.4 Insert a column into a matrix

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Fortran .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} Using Matrix/Vector .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Python .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.5 Build matrix from other matrices and vectors

2.5.1 First example
2.5.2 second example

#### 2.5.1 First example

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

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

The command ArrayFlatten is essential for this in Mathematica.

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $$\left ({\begin{array}{cc} 1 & 4 \\ 2 & 5 \\ 3 & 6 \\ 7 & 10 \\ 8 & 11 \\ 9 & 12 \\ \end{array}} \right )$$

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} \left [{\begin{array}{cc} 1&4\\ \noalign{\medskip }2&5 \\ \noalign{\medskip }3&6\\ \noalign{\medskip }7&10\\ \noalign{\medskip } 8&11\\ \noalign{\medskip }9&12\end{array}} \right ]

 Python .syntaxhighlighter textarea {font-size: 14px !important;} Another way .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Fortran .syntaxhighlighter textarea {font-size: 14px !important;} Using the RESHAPE command .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

#### 2.5.2 second example

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

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

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

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

Mathematica

Maple

Matlab

Fortran

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Or

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

Fortran

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

### 2.7 Generate an n by m zero matrix

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} $\left [{\begin{array}{cccc} 0&0&0&0\\ 0&0&0&0 \\ 0&0&0&0\end{array}} \right ]$

 Fortran .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Ada .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.8 Rotate a matrix by 90 degrees

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Fortran Using additional copy for the matrix .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.9 Generate a diagonal matrix with given values on the diagonal

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.10 Sum elements in a matrix along the diagonal

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} Another ways .syntaxhighlighter textarea {font-size: 14px !important;} 15

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} Out[49]= 45

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} ans = 45

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} 45

### 2.12 Check if a Matrix is diagonal

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} true

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.14 Replicate a matrix

Given Matrix

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.15 Find the location of the maximum value in a matrix

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple This below ﬁnds position of ﬁrst max. .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Maple support for such operations seems to be not as strong as Matlab. One way to ﬁnd locations of all elements is by using explicit loop .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.16 Swap 2 columns in a matrix

Give a matrix

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

 Mathematica In Mathematica, to join 2 matrices side-by-side, use Join with '2' as the third argument. To join them one on top of the other, use '1' as the third argument .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;}

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

Question posted on the net

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.19 extract values from matrix given their index

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

Example, given

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Another method (same really as above, but using Part explicit) .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.20 Convert $$N$$ by $$M$$ matrix to a row of length $$N M$$

Given

covert the matrix to one vector

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple Maple reshapes along columns, like Matlab. To get same result as Mathematica, we can transpose the matrix ﬁrst. To get same result as Matlab, do not transpose. .syntaxhighlighter textarea {font-size: 14px !important;} Notice the result is a row matrix and not a vector. To get a vector .syntaxhighlighter textarea {font-size: 14px !important;} They look the same on the screen, but using whattype we can ﬁnd the type. .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.21 ﬁnd rows in a matrix based on values in diﬀerent columns

Example, given Matrix

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

Given

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} another way is to ﬁnd the index using Position and then use Extract .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} another way is to use Cases[]. This is the shortest way .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.23 Locate rows in a matrix with column being a string

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

and given list

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} But notice that in Matlab, the answer is a cellarray. To access the numbers above .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

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

Assume valid list of indices.

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

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

 Mathematica Three methods are shown. method 1: (credit for the idea to Mike Honeychurch at stackoverﬂow). It turns out it is easier to work with what we want to keep instead of what we want to delete so that Part[] can be used directly. Hence, given a list of row numbers to remove, such as .syntaxhighlighter textarea {font-size: 14px !important;} Start by generating list of the rows and columns to keep by using the command Complement[], followed by using Part[] .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} method 2: (due to Mr Wizard at stackoverﬂow) .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} method 3: (me) use Pick. This works similar to Fortran pack(). Using a mask matrix, we set the entry in the mask to False for those elements we want removed. Hence this method is just a matter of making a mask matrix and then using it in the Pick[] command. .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.25 Convert list of separated numerical numbers to strings

Problem: given a list of numbers such as

convert the above to list of strings

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} answer below is due to Bruno Luong at Matlab newsgroup .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.26 Obtain elements that are common to two vectors

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

Given

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

Given

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.29 Sort a matrix row-wise using ﬁrst column as key

Given

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.30 Sort a matrix row-wise using non-ﬁrst column as key

Given

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

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.31 Replace the ﬁrst nonzero element in each row in a matrix by some value

Problem: Given a matrix, replace the ﬁrst nonzero element in each row by some a speciﬁc value. This is an example. Given matrix $$A$$ below, replace the ﬁrst non-zero element in each row by $$-1$$, then

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} Solution due to Bob Hanlon (from Mathematica newsgroup) .syntaxhighlighter textarea {font-size: 14px !important;} Solution by Fred Simons (from Mathematica newsgroup) .syntaxhighlighter textarea {font-size: 14px !important;} Solution due to Adriano Pascoletti (from Mathematica newsgroup) .syntaxhighlighter textarea {font-size: 14px !important;} Solution due to Oliver Ruebenkoenig (from Mathematica newsgroup) .syntaxhighlighter textarea {font-size: 14px !important;} Solution due to Szabolcs Horvát (from Mathematica newsgroup) .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} This solution due to Bruno Luong (from matlab newsgroup) .syntaxhighlighter textarea {font-size: 14px !important;} This solution due to Jos from matlab newsgroup .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.32 Perform outer product and outer sum between two vector

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

 Mathematica using symbolic vectors. Outer product .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Outer sum .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} using numerical vectors. Outer product .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Outer sum .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab Outer product .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Outer sum .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple Due to Carl Love from the Maple newsgroup .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} 3,4 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.35 Solve $$Ax=b$$

Solve for x in the following system of equations

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Fortran

compile and run

### 2.36 Find all nonzero elements in a matrix

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

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

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

 Mathematica In Mathematica, standard Mathematica matrix operations can be used, or the matrix can be converted to SparseArray and special named operation can be used on it. .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Or standard list operations can be used .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.37 evaluate f(x) on a vector of values

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.38 generates equally spaced N points between $$x_1$$ and $$x_2$$

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

Mathematica

The above can also be done using Plot3D

I need to sort out the orientation diﬀerence between the two plots above.

Matlab

### 2.40 Find determinant of matrix

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.41 Generate sparse matrix with $$n$$ by $$n$$ matrix repeated on its diagonal

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{ccccccccc} 1. & 2. & 3. & 0 & 0 & 0 & 0 & 0 & 0 \\ 4. & 5. & 6. & 0 & 0 & 0 & 0 & 0 & 0 \\ 7. & 8. & 9. & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1. & 2. & 3. & 0 & 0 & 0 \\ 0 & 0 & 0 & 4. & 5. & 6. & 0 & 0 & 0 \\ 0 & 0 & 0 & 7. & 8. & 9. & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1. & 2. & 3. \\ 0 & 0 & 0 & 0 & 0 & 0 & 4. & 5. & 6. \\ 0 & 0 & 0 & 0 & 0 & 0 & 7. & 8. & 9. \\ \end{array}} \right )$

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.42 Generate sparse matrix for the tridiagonal representation of second diﬀerence operator in 1D

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cccccccccccc} -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 \\ \end{array}} \right )$

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.43 Generate sparse matrix for the Laplacian diﬀerential operator $$\nabla ^{2}u$$ for 2D grid

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

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.44 Generate sparse matrix for the Laplacian diﬀerential operator $$\nabla ^{2}u$$ for 3D grid

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

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

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

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

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

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

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

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

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

Now that the $$A$$ structure is understood, the A matrix can be generated.

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

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Example 2: Using $$n_x=2$$, $$n_y=2$$, $$n_z=3$$.

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Example 3: Using $$n_x=3$$, $$n_y=3$$, $$n_z=3$$.

Matlab

Matlab result is below.

### 2.45 Generate the adjugate matrix for square matrix

Given square matrix such as

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

ﬁnd the adjugate matrix which is

$\begin{pmatrix} -3 & 6 & -3 \\ 6 & -12 & 6 \\ -3 & 6 & -3 \\ \end{pmatrix}$

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab Will try to ﬁnd function in Matlab to do this. But for non-singular matrices only the direct method of ﬁnding the inverse and then multiplying by the determinant to recover the adjunct can be used. .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} The following is due to Matt J from the matlab newsgroup .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Fortran

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

compile and run

### 2.46 Multiply each column by values taken from a row

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

In Matlab, bsxfun is used.

 Mathematica credit for this solution goes to Bob Hanlon, Adriano Pascoletti, Kurt Tekolste, David Park, and Peter. J. C. Moses from the Math group .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Another way is to use Inner[] command. Credit for this solution goes to Sswziwa Mukasa and Peter. J. C. Moses from the Math group .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Fortran

Octave

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

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

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Fortran

compile and run

compile and run

### 2.48 delete one row from a matrix

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} or a little longer solution using Pick .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.49 delete one column from a matrix

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} or a little longer solution using Pick .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.50 generate random matrix so that each row adds to 1

Generate the random matrix and divide each row by its total

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.51 generate random matrix so that each column adds to 1

Generate the random matrix and divide each column by its total

 Mathematica This method due to Bob Hanlon from the Math group .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Or can use Transpose .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Another way of doing the above, without the need to transpose 2 times is the following .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Maple

### 2.52 sum all rows in a matrix

Given the folowing 2D matrix $$A$$ ﬁnd the sum of each row

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.53 sum all columns in a matrix

Given the folowing 2D matrix $$A$$ ﬁnd the sum of each column

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.54 ﬁnd in which columns values that are not zero

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} 3,5,8,9,20.2,30,-44 .syntaxhighlighter textarea {font-size: 14px !important;} 3,5,8,9,20.2,30,-44

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Fortran

### 2.56 How to ﬁnd mean of equal sized segments of a vector

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.57 ﬁnd ﬁrst value in column larger than some value and cut matrix from there

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

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} Another method that can be used in this, but I prefer the kron method above: .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.59 repeat each column of matrix number of times

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} Another method that can be used in this, but the above is better .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

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

generate the matrix

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

 Mathematica The trick is to use Map with 2 at end. .syntaxhighlighter textarea {font-size: 14px !important;} Matlab to do

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

Given a vector or list $$1,2,3,4,5,6,7,8,9,10$$ how to ﬁnd the sum of all its elements?

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} 55

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} 55

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} 55

### 2.62 How to ﬁnd maximum of each row of a matrix?

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.63 How to ﬁnd maximum of each column of a matrix?

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.66 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 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.67 Check if a matrix is Hermite

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

A Matrix is Hermite if it is the same as the conjugate of its transpose. One way to check is to take the 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 .syntaxhighlighter textarea {font-size: 14px !important;} Now check that each element in the diﬀerence matrix is less than MachineEpsilon .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} 1

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.69 Linear convolution of 2 sequences

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.70 Circular convolution of two sequences

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

#### 2.70.1 example 1. The sequences are of equal length

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

#### 2.70.2 example 2. The sequences are of unequal length

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.71 Linear convolution of 2 sequences with origin at arbitrary position

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

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

Mathematica

case 1

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

case 2

Convolve 2 sequences where the origin is located at diﬀerent location from the ﬁrst element of the sequence, use DiracDelta function, and the DTFT approach.

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

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

Write down the x sequence and take its DTFT

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

Now multiply the DTFT’s and take the inverse

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

### 2.72 Visualize a 2D matrix

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

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

Mathematica

Matlab

Maple

### 2.73 Find the cross correlation between two sequences

Problem: Given \begin{align*} A & =[0,0,2,-1,3,7,1,2,-3,0,0]\\ B & =[0,0,1,-1,2,-2,4,1,-2,5,0,0] \end{align*}

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

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Matlab In MATLAB use the xcross in the signal processing toolbox .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

With MATLAB use the orth(A) function, With Mathematica, use {u,s,v}=SingularValueDecomposition[A] , and since the rank is 2, then the ﬁ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 .syntaxhighlighter textarea {font-size: 14px !important;} 2 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.75 Solve A x= b and display the solution

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.76 Determine if a set of linear equations A x= b, has a solution and what type of solution

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

The following algorithm summarizes all the cases

Let the system of equations be \begin{align*} y & =x-1\\ y & =2x+1 \end{align*}

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

Mathematica

2

2

The above algorithm can now be run as follows

The output of the above is

Matlab

Output is

### 2.77 Given a set of linear equations automatically generate the matrix $$A$$ and vector $$b$$ and solve $$A x=b$$

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

Automatically convert it to the form $$Ax=b$$ and solve

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{ccc} 4 & 3 & 2 \\ 5 & 6 & 3 \\ 9 & 7 & 10 \\ \end{array}} \right )$ .syntaxhighlighter textarea {font-size: 14px !important;} $\{-12,-7,-9\}$ .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 2.78 Convert a matrix to row echelon form and to reduced row echelon form

Problem: Given a matrix A, convert it to REF and RREF. Below shows how to

convert the matrix A to RREF. To convert to REF (TODO). One reason to convert Matrix $$A$$ to its row echelon form, is to ﬁnd the rank of $$A$$. If matrix $$A$$ is a $$4\times 4$$, and when converted to its row echelon form we ﬁnd that one of the rows is all zeros, then the rank of $$A$$ will be 3 and not full rank.

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cccc} 1 & 1 & -2 & 1 \\ 3 & 2 & 4 & -4 \\ 4 & 3 & 3 & -4 \\ \end{array}} \right )$ .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cccc} 1 & 0 & 0 & 2 \\ 0 & 1 & 0 & -3 \\ 0 & 0 & 1 & -1 \\ \end{array}} \right )$

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} \left [{\begin{array}{cccc} 1&0&0&2\\ \noalign{\medskip }0&1&0&-3 \\ \noalign{\medskip }0&0&1&-1\end{array}} \right ]

### 2.79 Convert 2D matrix to show the location and values

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

Generate the matrix

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

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

Mathematica

Another way

But I think the simplist is to use Table

Matlab

gives

### 2.80 How to ﬁnd equation of line which is the intersection of 2 planes in 3D?

Given plane $$2 x+3 y+z-6=0$$ and plane $$x-y-2 z+8=0$$, we want to ﬁnd equation of the line, which is the intersection of these two planes.

The idea is to solve the above two equations for $$x,y$$ and parametrise the solution in terms of $$z=t$$.

Mathematica

Now ParametricPlot3D is used to plot the line

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

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

Problem: Generate one signal of some width and height.

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;}

Matlab

### 3.2 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 .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

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

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple Maple need an Array for input, not list, so have to convert this is little strange .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

Fortran

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

### 3.4 Generate uniform distributed random numbers

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Fortran

compile and run

#### 3.4.2 How to generate 5 random numbers from a to b?

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Fortran

compile and run

#### 3.4.3 How to generate MATRIX of random numbers from a to b?

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

Mathematica

Matlab

Fortran

compile and run

### 3.5 Determine and plot the Fourier Transform for continuous time function

Plot the Fourier Transform of $$3 \sin (t)$$

Mathematica

Matlab

Need to do the plot.

### 3.6 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 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} now load the samples and do power spectrum .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

Matlab

### 3.7 Plot the power spectral density of a signal

Problem: Given a signal $3\sin \left ( 2\pi f_{1}t\right ) +4\cos \left ( 2\pi f_{2}t\right ) +5\cos \left ( 2\pi f_{3}t\right )$ for the following frequency values (in Hz) $f_{1}=20 ,f_{2}=30, f_{3}=50$ and plot the signal power spectral density so that 3 spikes will show at the above 3 frequencies. Use large enough sampling frequency to obtain a sharper spectrum

Mathematica

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

And using InterpolationOrder -> 2

Matlab

Mathematica

Matlab

### 3.9 Obtain Fourier Series coeﬃcients for a periodic function

Given a continuous time function $$x\left ( t\right )$$, ﬁnd its Fourier coeﬃcients $$c_{n}$$, where $x\left ( t\right ) =\sum _{n=-\infty }^{\infty }c_{n}e^{j\omega _{0}nt}$ and $c_{n}=\frac{1}{T_{0}}\int _{-\frac{T_{0}}{2}}^{\frac{T_{0}}{2}}x(t)e^{-j\omega _{0}nt}dt$ Where $$T_{0}$$ is the fundamental period of $$x\left ( t\right )$$ and $$\omega _{0}=\frac{2\pi }{T_{0}}$$

Mathematica

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

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

Plot phase

Find Fourier series for periodic square wave

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

Show phase

### 3.10 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

### 3.11 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 .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;}

### 3.12 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 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 3.13 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 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 3.14 How to ﬁnd cummulative sum

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

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 3.15 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$$ .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} For window size $$k=3$$ .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab For a window size $$k=2$$ .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} For window size $$k=3$$ .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 3.16 How to select N random values from a set of numbers?

Given the set $$v{1,2,3,5,6,7,11,12,20,21}$$ how to select say $$m=5$$ random numbers from it?

 Mathematica method 1 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} method 2 (Version 9) .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

### 3.17 How to sample a sin signal and plot it?

Sample a sin signal of one second period at 8 times its frequency.

Mathematica

Matlab

### 3.18 How ﬁnd the impulse response of a diﬀerence equation?

Find the impulse response $$h[n]$$ for the discrete system given by the diﬀerence equation $$y[n]-\frac{1}{2} y[n-1]=x[n]-\frac{1}{4}x[n-1]$$

analytical solution:

1.
The ﬁrst step is to replace $$y[n]$$ by $$h[n]$$ and $$x[n]$$ by $$\delta [n]$$ giving $h[n]-\frac{1}{2} h[n-1]=\delta [n]-\frac{1}{4}\delta [n-1]$
2.
Taking the Fourier transform and assuming $$H(e^{i \omega })$$ is the Fourier transform of $$h[n]$$ results in $H(e^{i \omega })-\frac{1}{2} H(e^{i \omega }) e^{- i \omega } = 1 -\frac{1}{4} e^{- i \omega }$
3.
Solving for $$H(e^{i \omega })$$ gives \begin{align*} H(e^{i \omega }) &= \frac{ 1 -\frac{1}{4} e^{- i \omega } }{ 1- \frac{1}{2} e^{-i \omega } }\\ &= \frac{1}{1- \frac{1}{2} e^{-i \omega }} - \frac{1}{4} \frac{ e^{-i \omega }}{1- \frac{1}{2} e^{-i \omega }} \end{align*}
4.
Taking the inverse discrete Fourier transform given by $h[n]= \frac{1}{2\pi } \int _{-\pi }^{\pi } H(e^{i \omega }) e^{i \omega }$ which results in $h[n]= \left ( \frac{1}{2} \right )^n u[n] - \frac{1}{4} \left ( \frac{1}{2} \right )^{n-1} u[n-1]$

Mathematica

${\begin{array}{cc} -2^{-n-1} & n>0 \\ 0 & \text{True} \\ \end{array}} +{\begin{array}{cc} 2^{-n} & n\geq 0 \\ 0 & \text{True} \\ \end{array}}$

And some values of $$h[n]$$ starting from $$n=0$$ are

$\left \{1,\frac{1}{4},\frac{1}{8},\frac{1}{16},\frac{1}{32},\frac{1}{64},\frac{1}{128}, \frac{1}{256},\frac{1}{512}, \frac{1}{1024},\frac{1}{2048}\right \}$

## Chapter 4Diﬀerential, PDE solving, integration, numerical and analytical solving of equations

### 4.1 Generate direction ﬁeld plot of a ﬁrst order diﬀerential equation

Problem: Given the following non autonomous diﬀerential equation, plot the line ﬁelds which represents the solutions of the ODE.

$\frac{dy\left ( x\right ) }{dx}=x^2 - y$

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

The following are the steps to generate direction ﬁeld plot for $$\frac{dy}{dx}=f(x,y)$$

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

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

In Mathematica, the command VectorPlot is used. In Maple dfieldplot is used.

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} $y \left ( x \right ) =-{\frac{\sqrt{3}{{\rm Ai}^{(1)}\left (x\right )}+{{\rm Bi}^{(1)}\left (x\right )}}{\sqrt{3}{{\rm Ai}\left (x\right )}+{{\rm Bi}\left (x\right )}}}$ .syntaxhighlighter textarea {font-size: 14px !important;}

### 4.2 Solve the 2-D Laplace PDE for a rectangular plate with Dirichlet boundary conditions

Problem: Solve $$\nabla ^{2}T\left (x,y\right )=0$$ on 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 .syntaxhighlighter textarea {font-size: 14px !important;} Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

### 4.3 Solve homogeneous 1st order ODE, 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 .syntaxhighlighter textarea {font-size: 14px !important;} Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

### 4.4 Solve homogeneous 2nd order ODE with constant coeﬃcients

Problem: Solve $y^{\prime \prime }\left ( t\right ) -1.5y^{\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.5y^{\prime }\left ( t\right ) +5y\left ( t\right ) =0$$ let \begin{align*} x_{1} & =y\\ x_{2} & =y^{\prime }\\ & =x_{1}^{\prime } \end{align*}

hence $x_{1}^{\prime }=x_{2}$ and \begin{align*} x_{2}^{\prime } & =y^{\prime \prime }\\ & =1.5y^{\prime }-5y\\ & =1.5x_{2}-5x_{1} \end{align*}

Hence we can now write $\begin{bmatrix} x_{1}^{\prime }\\ x_{2}^{\prime }\end{bmatrix} =\begin{bmatrix} 0 & 1\\ -5 & 1.5 \end{bmatrix}\begin{bmatrix} x_{1}\\ x_{2}\end{bmatrix}$ Now Matlab ODE45 can be used.

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $$1. \left (1. e^{0.75 t} \cos (2.10654 t)-0.356034 e^{0.75 t} \sin (2.10654 t)\right )$$ .syntaxhighlighter textarea {font-size: 14px !important;} Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

### 4.5 Solve non-homogeneous 2nd order ODE, constant coeﬃcients

Problem: Solve $y^{\prime \prime }\left ( t\right ) -1.5y^{\prime }\left ( t\right ) +5y\left ( t\right ) =4\cos \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.5y^{\prime }\left ( t\right ) +5y\left ( t\right ) =4\cos \left ( t\right )$$, let \begin{align*} x_{1} & =y\\ x_{2} & =y^{\prime }\\ & =x_{1}^{\prime } \end{align*}

hence $x_{1}^{\prime }=x_{2}$ and \begin{align*} x_{2}^{\prime } & =y^{\prime \prime }\\ & =1.5y^{\prime }-5y+4\cos \left ( t\right ) \\ & =1.5x_{2}-5x_{1}+4\cos \left ( t\right ) \end{align*}

Hence we can now write $\begin{bmatrix} x_{1}^{\prime }\\ x_{2}^{\prime }\end{bmatrix} =\begin{bmatrix} 0 & 1\\ -5 & 1.5 \end{bmatrix}\begin{bmatrix} x_{1}\\ x_{2}\end{bmatrix} +\begin{bmatrix} 0\\ 4 \end{bmatrix} \cos \left ( t\right )$

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $$\frac{1}{5183}( -1704 \sin (t) \sin ^2\left (\frac{\sqrt{71} t}{4}\right )+69 \sqrt{71} e^{3 t/4} \sin \left (\frac{\sqrt{71} t}{4}\right )+4544 \cos (t) \cos ^2\left (\frac{\sqrt{71} t}{4}\right )+639 e^{3 t/4} \cos \left (\frac{\sqrt{71} t}{4}\right )-1704 \sin (t) \cos ^2\left (\frac{\sqrt{71} t}{4}\right )+4544 \sin ^2\left (\frac{\sqrt{71} t}{4}\right ) \cos (t)))$$ .syntaxhighlighter textarea {font-size: 14px !important;} Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

### 4.6 Solve homogeneous 2nd order ODE, constant coeﬃcients (BVP)

Problem: Solve $y^{\prime \prime }\left ( t\right ) +t\ 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\ y\left ( t\right ) =0$$, let \begin{align*} x_{1} & =y\\ x_{2} & =y^{\prime }\\ & =x_{1}^{\prime } \end{align*}

Therefore $x_{1}^{\prime }=x_{2}$ And \begin{align*} x_{2}^{\prime } & =y^{\prime \prime }\\ & =-t\ y\\ & =-t\ x_{1} \end{align*}

This results in $\begin{bmatrix} x_{1}^{\prime }\\ x_{2}^{\prime }\end{bmatrix} =\begin{bmatrix} 0 & 1\\ -t & 0 \end{bmatrix}\begin{bmatrix} x_{1}\\ x_{2}\end{bmatrix}$ Now bvp4c() can be used.

Mathematica

$\frac{-\sqrt{3} \text{Ai}\left (\sqrt [3]{-1} t\right )+\text{Bi}\left (\sqrt [3]{-1} t\right )-3\ 3^{2/3} \Gamma \left (\frac{2}{3}\right ) \text{Bi}\left (20 \sqrt [3]{-1}\right ) \text{Ai}\left (\sqrt [3]{-1} t\right )+3\ 3^{2/3} \Gamma \left (\frac{2}{3}\right ) \text{Ai}\left (20 \sqrt [3]{-1}\right ) \text{Bi}\left (\sqrt [3]{-1} t\right )}{\sqrt{3} \text{Ai}\left (20 \sqrt [3]{-1}\right )-\text{Bi}\left (20 \sqrt [3]{-1}\right )}$

Matlab

### 4.7 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 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;}

### 4.8 Show the eﬀect of boundary/initial conditions on 1-D heat 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

Each plot shows the boundary conditions used.

### 4.9 Find particular and homogenous solution to undetermined system of equations

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

$\begin{bmatrix} 2 & 4 & 6 & 4\\ 2 & 5 & 7 & 6\\ 2 & 3 & 5 & 2 \end{bmatrix}\begin{bmatrix} x_{1}\\ x_{2}\\ x_{3}\\ x_{3}\end{bmatrix} =\begin{bmatrix} b_{1}\\ b_{2}\\ b_{3}\end{bmatrix} \ where\begin{bmatrix} b_{1}\\ b_{2}\\ b_{3}\end{bmatrix} =\begin{bmatrix} 4\\ 3\\ 5 \end{bmatrix}$

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$$\backslash$$b uses least square approach to determine a solution.  I am sure there is a way, will update this once I ﬁnd out.

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cccc} 2 & 4 & 6 & 4 \\ 2 & 5 & 7 & 6 \\ 2 & 3 & 5 & 2 \\ \end{array}} \right )$ .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Mathematica LinearSolve gives one solution (particular solution) .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} ﬁnd the homogenous solution (nullspace) and add it to the above particular solution .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cc} 2 & -1 \\ -2 & -1 \\ 0 & 1 \\ 1 & 0 \\ \end{array}} \right )$ Add the particular solution to the homogenous solution to get the general solution .syntaxhighlighter textarea {font-size: 14px !important;} $\left ({\begin{array}{cc} -2 & -5 \\ -1 & 0 \\ 0 & 1 \\ 1 & 0 \\ \end{array}} \right )$ To obtain the general solution right away Solve can be used instead .syntaxhighlighter textarea {font-size: 14px !important;} $$\left \{\text{x3}\to -\frac{\text{x1}}{2}-\frac{\text{x2}}{2}+\frac{3}{2},\text{x4}\to \frac{\text{x1}}{4}-\frac{\text{x2}}{4}-\frac{5}{4}\right \}$$ .syntaxhighlighter textarea {font-size: 14px !important;} $$\left \{\text{x1},\text{x2},-\frac{\text{x1}}{2}-\frac{\text{x2}}{2}+\frac{3}{2},\frac{\text{x1}}{4}-\frac{\text{x2}}{4}-\frac{5}{4}\right \}$$

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Maple .syntaxhighlighter textarea {font-size: 14px !important;} A=\left [{\begin{array}{cccc} 2&4&6&4\\ \noalign{\medskip }2&5&7&6\\ \noalign{\medskip }2&3&5&2\end{array}} \right ] b=\left [{\begin{array}{c} 4\\ \noalign{\medskip }3\\ \noalign{\medskip }5\end{array}} \right ] .syntaxhighlighter textarea {font-size: 14px !important;} \left [{\begin{array}{c} 4-x_{{3}}+2\,x_{{4}}\\ \noalign{\medskip }-1-x_{{3}}-2\,x_{{4}}\\ \noalign{\medskip }x_{{3}}\\ \noalign{\medskip }x_{{4}}\end{array}} \right ] Can solve this system to get the general solution using the solve command as follows .syntaxhighlighter textarea {font-size: 14px !important;} WARNING: Maple sometimes reorders the result from solve() so we can get a diﬀerent ordering of the free variables as shown above. \left [{\begin{array}{c} 4-x_{{3}}+2\,x_{{4}}\\ \noalign{\medskip }-1-x_{{3}}-2\,x_{{4}}\\ \noalign{\medskip }x_{{3}}\\ \noalign{\medskip }x_{{4}}\end{array}} \right ]

### 4.10 Plot the constant energy levels for a nonlinear pendulum

Problem:

Plot the constant energy levels for the nonlinear pendulum in $$\theta ,\dot{\theta }$$

Assume that $$m=1,g=9.8m/s^{2},L=10m$$

The constant energy curves are curves in the Y-X plane where energy is constant. The Y-axis represents $$\dot{\theta }$$, and the X-axis represents $$\theta$$

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

The higher the energy the pendulum posses, the larger the angle $$\theta$$ it will swing by will be.

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

The expression for the energy $$E$$ for the pendulum is ﬁrst derived as a function of $$\theta ,\dot{\theta }$$\begin{align*} E & =PE+KE\\ & =mgL\left ( 1-\cos \theta \right ) +\frac{1}{2}mL^{2}\dot{\theta }^{2} \end{align*}

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

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

The following is the resulting plot for $$m=1,g=9.8m/s^{2} ,L=10m$$

Mathematica

Matlab

Maple

The Maple solution was contributed by Matrin Eisenberg

### 4.11 Solve numerically the ODE $$u^{''''}+u=f$$ using point collocation method

Problem: Give the ODE $\frac{d^{4}u\left ( x\right ) }{dx^{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 ) ={\displaystyle \sum \limits _{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 ) ={\displaystyle \sum \limits _{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 ) }{dx^{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

\begin{equation*} \begin{split} 0.0412493 x-0.0826459 x^3+0.0416667 x^4 -0.00034374 x^5-1.51283*10^-8 x^6\\ +0.0000984213 x^7-0.0000248515 x^8 +1.64129*10^-7 x^9-3.28259*10^-8 x^10 \end{split} \end{equation*}

$-\frac{\cos \left (\frac{1-(1+i) x}{\sqrt{2}}\right )+\cos \left (\frac{1-(1-i) x}{\sqrt{2}}\right )+\cosh \left (\frac{1-(1+i) x}{\sqrt{2}}\right )+\cosh \left (\frac{1-(1-i) x}{\sqrt{2}}\right )-2 \left (\cos \left (\frac{1}{\sqrt{2}}\right )+\cosh \left (\frac{1}{\sqrt{2}}\right )\right )}{2 \left (\cos \left (\frac{1}{\sqrt{2}}\right )+\cosh \left (\frac{1}{\sqrt{2}}\right )\right )}$

$-\frac{-i \cos \left (\frac{1-(1+i) x}{\sqrt{2}}\right )+i \cos \left (\frac{1-(1-i) x}{\sqrt{2}}\right )+i \cosh \left (\frac{1-(1+i) x}{\sqrt{2}}\right )-i \cosh \left (\frac{1-(1-i) x}{\sqrt{2}}\right )}{2 \left (\cos \left (\frac{1}{\sqrt{2}}\right )+\cosh \left (\frac{1}{\sqrt{2}}\right )\right )}$

Maple

### 4.12 How to numerically solve a set of non-linear equations?

Numerically solve the following three equations for $$e,a,f$$ \begin{align*} eq1 &= r1 - a(e-1) \\ eq2 &= delT - \sqrt{\frac{a^3}{\mu }}(e*\sinh (f)-f) \\ eq3 &= r2 - a(e \cosh (f)-1) \end{align*}

 Mathematica .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

 Matlab .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} Another option is solve but slower .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

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

Problem: Solve $y^{\prime \prime }\left ( t\right ) -\mu \left ( 1-y^{2}\right ) y^{\prime }\left ( t\right ) +y\left ( t\right ) =0$ for diﬀerent values of $$\mu$$ to compare eﬀect of changing $$\mu$$ on the solution trajectories. The initial conditions are \begin{align*} y\left ( 0\right ) & =0.1\\ y^{\prime }\left ( t\right ) & =0 \end{align*}

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

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

Starting by writing the 2nd order ODE as 2 ﬁrst order ODE’s

$$\left .{\begin{array}[c]{l}x_{1}=y\\ x_{2}=y^{\prime }\end{array}} \right \}$$derivatives$$\Rightarrow \left .{\begin{array} [c]{l}x_{1}^{^{\prime }}=y^{\prime }\\ x_{2}^{^{\prime }}=y^{^{\prime \prime }}\ \end{array}} \right \} \Rightarrow \left .{\begin{array} [c]{l}x_{1}^{^{\prime }}=x_{2}\\ x_{2}^{^{\prime }}=\mu \left ( 1-x_{1}^{2}\right ) x_{2}+x_{1}\end{array}} \right \}$$

Below is the solution of the diﬀerential equation for diﬀerent values of $$\mu =1$$

Mathematica

Matlab

### 4.14 How to numerically solve Poisson PDE on 2D using Jacobi iteration method?

Problem: Solve $$\bigtriangledown ^2 u= f(x,y)$$ on 2D using Jacobi method. Assume $$f(x,y)=-\mathrm{e}^{-(x - 0.25)^2 - (y - 0.6)^2}$$. Use mesh grid norm and relative residual to stop the iterations.

Mathematica

Matlab

### 4.15 How to solve BVP second order ODE using ﬁnite elements with linear shape functions and using weak form formulation?

Solve $$y''(x)-3 y(x) = -x^2$$ over $$x=0\ldots 1$$ with boundary conditions $$x(0)=0$$ and $$x(1)=0$$ using piecewise linear trial functions.

solution

Let the trial function be the linear function $$\hat{y}(x)=c_1 x+c_2$$. The residual is $$R=\hat{y}''(x)-3\hat{y}(x)+x^2$$. Let the test function be $$w(x)$$. We need to solve for each element $$i$$ the following equation $I_i = \int w_i \, R_i\,dx$ Using the weak form, we apply integration by parts to obtain \begin{align*} I_i &= \int _{x_i}^{x_{i+1}} -\frac{dw}{dx}\frac{d\hat{y}}{dx} -3 w(x)\hat{y}(x) + w(x) x^2 \,dx + \left [ w(x) \frac{d\hat{y}}{dx} \right ]_{x_i}^{x_{i+1}} \\ &=0 \tag{96.1}\\ \end{align*}

or \begin{align*} I &= \sum _{i=1}^{i=M}{I_i} \\ &= \sum _{i=1}^{i=M}{ \int _{x_i}^{x_{i+1}} \left ( -\frac{dw}{dx}\frac{d\hat{y}}{dx} -3 w(x)\hat{y}(x) + w(x) x^2 \,dx + \left [ w(x) \frac{d\hat{y}}{dx} \right ]_{x_i}^{x_{i+1}} \right ) }\\ &=0 \end{align*}

Where $$M$$ is the number of elements. The above is the $$M$$ equations we need to solve for the unknowns $$y_i$$ at the internal nodes. In Galerkin method, the test function is $$w_1(x)=H_1(x)$$ and $$w_2(x)=H_2(x)$$ which comes from writing $$w_i(x) = \frac{d\hat{y}(x)}{dy_i}$$

Rewriting the trial function $$\hat{y}(x)$$ in terms of unknown values at nodes results in $\hat{y}(x) = H_1(x) y_i + H_2(x) y_{i+1}$ Where $$H_1(x)=\frac{x_{i+1}-x}{h}$$ and $$H_2(x)=\frac{x-x_i}{h}$$ where $$h$$ is the length of each element. Hence (96.1) becomes \begin{align} I_i &= \int _{x_i}^{x_{i+1}} -\begin{Bmatrix}H_1'(x)\\H_2'(x)\end{Bmatrix} \begin{Bmatrix}H_1'(x) H_2'(x)\end{Bmatrix} \begin{Bmatrix}y_i \\y_{i+1}\end{Bmatrix} -3 \begin{Bmatrix}H_1(x)\\H_2(x)\end{Bmatrix} \begin{Bmatrix}H_1(x) H_2(x)\end{Bmatrix} \begin{Bmatrix}y_i \\y_{i+1}\end{Bmatrix} + \begin{Bmatrix}H_1(x)\\H_2(x)\end{Bmatrix} x^2 \,dx + \left [ w(x) \frac{d\hat{y}}{dx} \right ]_{x_i}^{x_{i+1}} \tag{96.2} \end{align}

The terms $$\left [ w(x) \frac{d\hat{y}}{dx} \right ]_{x_i}^{x_{i+1}}$$ reduce to $$\left [ w(x) \frac{d\hat{y}}{dx} \right ]_{0}^{1}$$ since intermediate values cancel each other leaving the two edge values over the whole domain.

But $$H_1'(x)=\frac{-x}{h}$$ and $$H_2'(x)=\frac{x}{h}$$. Plugging these into (96.2) and integrating gives the following \begin{align} I_i &= \frac{1}{h} \begin{bmatrix} -1 & 1\\ 1 & -1 \end{bmatrix} \begin{Bmatrix}y_i \\y_{i+1}\end{Bmatrix} -h \begin{bmatrix} 1 & \frac{1}{2} \\ \frac{1}{2} & 1 \end{bmatrix} \begin{Bmatrix}y_i \\y_{i+1}\end{Bmatrix} + \frac{1}{12 h} \begin{bmatrix} 3 x_i^4-4 x_i^3 x_{i+1} + x_{i+1}^4 \\ x_i^4-4 x_i x_{i+1}^3 + 3 x_{i+1}^4 \end{bmatrix} + \left [ w(x) \frac{d\hat{y}}{dx} \right ]_{0}^{1} \tag{96.3} \end{align}

The algebra above was done with Mathematica. In the following code, $$x_1$$ is $$x_i$$ and $$x_2$$ is $$x_{i+1}$$

 .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;} .syntaxhighlighter textarea {font-size: 14px !important;}

The above equation $$I_i$$ is calculated for each element $$i$$, resulting in $$M$$ equations where $$M$$ is the number of elements. Collecting all these local equations into a global stiﬀness matrix and then adjusting the global stiﬀness matrix for boundary conditions by eliminating corresponding rows and columns, it is then solved for the unknowns $$y_i$$ as a system $$Ax=f$$ using direct solver. The following code illustrates the method.

Matlab

The analytical solution to the ODE is below. We can see it is very close to the FEM solution above with 11 elements. More elements leads to more accurate solution.

Mathematica

Analtical solution

Numerical

### 4.16 How to solve Poisson PDE in 2D using ﬁnite elements methods using rectangular element?

4.16.2 Case 1
4.16.3 $$f(x,y)=1$$
4.16.4 $$f(x,y)=xy$$
4.16.5 case 3

4.16.7 Case 4
4.16.8 Case 5
4.16.9 case 6

Solve $$\nabla ^{2}u=f\left ( x,y\right )$$ on square using FEM. Assume $$u=0$$ on boundaries and solve using $$f\left ( x,y\right ) =xy$$ and also $$f\left ( x,y\right ) =1$$ and $$f\left ( x,y\right ) =3\left ( \cos (4 \pi x)+\sin (3 \pi y)\right )$$

Use Galerkin method and weak form, Using a bilinear trial function. Let width of square be 1.

Solution

Using as an example with 9 elements to illustrate the method. The program below can be called with diﬀerent number of elements.

The trial function is bilinear

$$\tilde{u}=c_{1}+c_{2}x+c_{3}y+c_{4}xy \tag{1}$$

Looking at one element, and using local coordinates systems with element having width $$2a$$ and height $$2b$$ gives

Evaluating the trial function (1) at each corner node of the above element gives

\begin{align*} \tilde{u}_{1} & =c_{1}-ac_{2}-bc_{3}+abc_{4}\\ \tilde{u}_{2} & =c_{1}+ac_{2}-bc_{3}-abc_{4}\\ \tilde{u}_{3} & =c_{1}+ac_{2}+by_{3}+abc_{4}\\ \tilde{u}_{4} & =c_{1}-ac_{2}+bc_{3}+abc_{4} \end{align*}

Or

$\begin{Bmatrix} \tilde{u}_{1}\\ \tilde{u}_{2}\\ \tilde{u}_{3}\\ \tilde{u}_{4}\end{Bmatrix} =\begin{pmatrix} 1 & -a & -b & ab\\ 1 & a & -b & -ab\\ 1 & a & b & ab\\ 1 & -a & b & ab \end{pmatrix}\begin{Bmatrix} c_{1}\\ c_{2}\\ c_{3}\\ c_{4}\end{Bmatrix}$

Hence

\begin{align} \begin{Bmatrix} c_{1}\\ c_{2}\\ c_{3}\\ c_{4}\end{Bmatrix} & =\begin{pmatrix} 1 & -a & -b & ab\\ 1 & a & -b & -ab\\ 1 & a & b & ab\\ 1 & -a & b & ab \end{pmatrix} ^{-1}\begin{Bmatrix} \tilde{u}_{1}\\ \tilde{u}_{2}\\ \tilde{u}_{3}\\ \tilde{u}_{4}\end{Bmatrix} \nonumber \\ & =\frac{1}{4ab}\begin{pmatrix} ab & ab & ab & ab\\ -b & b & b & -b\\ -a & -a & a & a\\ 1 & -1 & 1 & -1 \end{pmatrix}\begin{Bmatrix} \tilde{u}_{1}\\ \tilde{u}_{2}\\ \tilde{u}_{3}\\ \tilde{u}_{4}\end{Bmatrix} \tag{2} \end{align}

Substituting (2) back into (1) and rearranging terms results in

\begin{align*} \tilde{u} & =c_{1}+c_{2}x+c_{3}y+c_{4}xy\\ & =\frac{1}{4ab}\left ( \left ( ab-bx-ay+xy\right ) \tilde{u}_{1}+\left ( ab+bx-ay-xy\right ) \tilde{u}_{2}+\left ( ab+bx+ay+xy\right ) +\left ( ab-bx+ay-xy\right ) \right ) \end{align*}

Since $$ab$$ is the $$\frac{1}{4}$$ of the area of element, then the above becomes

$\tilde{u}=\frac{1}{A}\left ( \left ( ab-bx-ay+xy\right ) \tilde{u}_{1}+\left ( ab+bx-ay-xy\right ) \tilde{u}_{2}+\left ( ab+bx+ay+xy\right ) +\left ( ab-bx+ay-xy\right ) \right )$

The above can now be written in term of what is called shape functions

$\tilde{u}(x,y)=N_{1}(x,y)\tilde{u_{1}}+N_{2}(x,y)\tilde{u_{2}}N_{3}(x,y)\tilde{u_{3}}+N_{4}(x,y)\tilde{u_{4}}$

Where

\begin{align*} N_{1} & =\frac{1}{A}\left ( ab-bx-ay+xy\right ) =\frac{1}{A}\left ( a-x\right ) \left ( b-y\right ) \\ N_{2} & =\frac{1}{A}\left ( ab+bx-ay-xy\right ) =\frac{1}{A}\left ( a+x\right ) \left ( b-y\right ) \\ N_{3} & =\frac{1}{A}\left ( ab+bx+ay+xy\right ) =\frac{1}{A}\left ( a+x\right ) \left ( b+y\right ) \\ N_{4} & =\frac{1}{A}\left ( ab-bx+ay-xy\right ) =\frac{1}{A}\left ( a-x\right ) \left ( b+y\right ) \end{align*}

Now that the shape functions are found, the next step is to determine the local element stiﬀness matrix. This can be found from the weak form integral over the area of the element.

$$I_{i}=\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}w\left ( \nabla ^{2}u-f\left ( x,y\right ) \right ) \,dxdy \tag{3}$$

Where $$w\left ( x,y\right )$$ is the test function. For Galerkin method, the test function $$w_{i}=\frac{d\tilde{u}}{du_{i}}=N_{i}\left ( x,y\right )$$. Hence

$w\left ( x,y\right ) =\begin{Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end{Bmatrix} =\begin{Bmatrix} \frac{1}{A}\left ( a-x\right ) \left ( b-y\right ) \\ \frac{1}{A}\left ( a+x\right ) \left ( b-y\right ) \\ \frac{1}{A}\left ( a+x\right ) \left ( b+y\right ) \\ \frac{1}{A}\left ( a-x\right ) \left ( b+y\right ) \end{Bmatrix}$

Using weak form, integration by parts is applied to (2)

$I_{i}=\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\left ( -\frac{\partial w}{\partial x_{i}}\frac{\partial \tilde{u}}{\partial x_{i}}-wf\left ( x,y\right ) \right ) \,dxdy+\overbrace{\int \limits _{\Gamma }w\frac{\partial \tilde{u}}{\partial n}}^{\text{goes to zero}}$

The term $$\int \limits _{\Gamma }w\frac{\partial \tilde{u}}{\partial n}$$ is the integration over the boundary of the element. Since there is only an essential boundary condition over all the boundaries (this is the given Dirchilet boundary condition), $$w=0$$ on the boundary and this integral vanishes. There is no natural boundary conditions for this example.

For those elements not on the external edge of the overall grid (i.e. internal elements), each contribution to this integral will cancel from the adjacent internal element. What this means is that the above reduces to just the ﬁrst integral

\begin{align*} I_{i} & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\left ( -\frac{\partial w}{\partial x_{i}}\frac{\partial \tilde{u}}{\partial x_{i}}-wf\left ( x,y\right ) \right ) \,dxdy\\ & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}-\begin{Bmatrix} \frac{\partial N_{1}}{\partial x}\\ \frac{\partial N_{2}}{\partial x}\\ \frac{\partial N_{3}}{\partial x}\\ \frac{\partial N_{4}}{\partial x}\end{Bmatrix}\begin{Bmatrix} \frac{\partial N_{1}}{\partial x} & \frac{\partial N_{2}}{\partial x} & \frac{\partial N_{3}}{\partial x} & \frac{\partial N_{4}}{\partial x}\end{Bmatrix}\begin{Bmatrix} \tilde{u}_{1}\\ \tilde{u}_{2}\\ \tilde{u}_{3}\\ \tilde{u}_{4}\end{Bmatrix} -\begin{Bmatrix} \frac{\partial N_{1}}{\partial y}\\ \frac{\partial N_{2}}{\partial y}\\ \frac{\partial N_{3}}{\partial y}\\ \frac{\partial N_{4}}{\partial y}\end{Bmatrix}\begin{Bmatrix} \frac{\partial N_{1}}{\partial y} & \frac{\partial N_{2}}{\partial y} & \frac{\partial N_{3}}{\partial y} & \frac{\partial N_{4}}{\partial y}\end{Bmatrix}\begin{Bmatrix} \tilde{u}_{1}\\ \tilde{u}_{2}\\ \tilde{u}_{3}\\ \tilde{u}_{4}\end{Bmatrix} \\ & -\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin{Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end{Bmatrix} f\left ( x,y\right ) \,dxdy \end{align*}

Hence

\begin{align*} I_{i} & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}-\begin{pmatrix} \frac{\partial N_{1}}{\partial x}\frac{\partial N_{1}}{\partial x}+\frac{\partial N_{1}}{\partial y}\frac{\partial N_{1}}{\partial y} & \frac{\partial N_{1}}{\partial x}\frac{\partial N_{2}}{\partial x}+\frac{\partial N_{1}}{\partial y}\frac{\partial N_{2}}{\partial y} & \frac{\partial N_{1}}{\partial x}\frac{\partial N_{3}}{\partial x}+\frac{\partial N_{1}}{\partial y}\frac{\partial N_{3}}{\partial y} & \frac{\partial N_{1}}{\partial x}\frac{\partial N_{4}}{\partial x}+\frac{\partial N_{1}}{\partial y}\frac{\partial N_{4}}{\partial y}\\ \frac{\partial N_{2}}{\partial x}\frac{\partial N_{1}}{\partial x}+\frac{\partial N_{2}}{\partial y}\frac{\partial N_{1}}{\partial y} & \frac{\partial N_{2}}{\partial x}\frac{\partial N_{2}}{\partial x}+\frac{\partial N_{2}}{\partial y}\frac{\partial N_{2}}{\partial y} & \frac{\partial N_{2}}{\partial x}\frac{\partial N_{3}}{\partial x}+\frac{\partial N_{2}}{\partial y}\frac{\partial N_{3}}{\partial y} & \frac{\partial N_{2}}{\partial x}\frac{\partial N_{4}}{\partial x}+\frac{\partial N_{2}}{\partial y}\frac{\partial N_{4}}{\partial y}\\ \frac{\partial N_{3}}{\partial x}\frac{\partial N_{1}}{\partial x}+\frac{\partial N_{3}}{\partial y}\frac{\partial N_{1}}{\partial y} & \frac{\partial N_{3}}{\partial x}\frac{\partial N_{2}}{\partial x}+\frac{\partial N_{3}}{\partial y}\frac{\partial N_{2}}{\partial y} & \frac{\partial N_{3}}{\partial x}\frac{\partial N_{3}}{\partial x}+\frac{\partial N_{3}}{\partial y}\frac{\partial N_{3}}{\partial y} & \frac{\partial N_{3}}{\partial x}\frac{\partial N_{4}}{\partial x}+\frac{\partial N_{3}}{\partial y}\frac{\partial N_{4}}{\partial y}\\ \frac{\partial N_{4}}{\partial x}\frac{\partial N_{1}}{\partial x}+\frac{\partial N_{4}}{\partial y}\frac{\partial N_{1}}{\partial y} & \frac{\partial N_{4}}{\partial x}\frac{\partial N_{2}}{\partial x}+\frac{\partial N_{4}}{\partial y}\frac{\partial N_{2}}{\partial y} & \frac{\partial N_{4}}{\partial x}\frac{\partial N_{3}}{\partial x}+\frac{\partial N_{4}}{\partial y}\frac{\partial N_{3}}{\partial y} & \frac{\partial N_{4}}{\partial x}\frac{\partial N_{4}}{\partial x}+\frac{\partial N_{4}}{\partial y}\frac{\partial N_{4}}{\partial y}\end{pmatrix}\begin{Bmatrix} \tilde{u}_{1}\\ \tilde{u}_{2}\\ \tilde{u}_{3}\\ \tilde{u}_{4}\end{Bmatrix} \\ & -\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin{Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end{Bmatrix} f\left ( x,y\right ) \,dxdy \end{align*}

In the above, we have the from $$I_{i}=\int \limits _{\Omega }-k_{i}\left \{ \tilde{u}\right \} d\Omega -\int \limits _{\Omega }\left \{ \tilde{u}\right \} f_{i}d\Omega$$, hence the element stiﬀness matrix is the ﬁrst integral, and the force vector comes from the second integral.

The integration is now carried out to obtain the element stiﬀness matrix. This was done using Mathematica. The local stiﬀness matrix for element $$i$$ is

\begin{align*} k_{i} & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\frac{\partial w}{\partial x_{i}}\frac{\partial \tilde{u}}{\partial x_{i}}dxdy\\ & =\begin{pmatrix} \frac{a^{2}+b^{2}}{3ab} & \frac{a}{6b}-\frac{b}{3a} & -\frac{a^{2}+b^{2}}{6ab} & -\frac{a}{3b}+\frac{b}{6a}\\ \frac{a}{6b}-\frac{b}{3a} & \frac{a^{2}+b^{2}}{3ab} & -\frac{a}{3b}+\frac{b}{6a} & -\frac{a^{2}+b^{2}}{6ab}\\ -\frac{a^{2}+b^{2}}{6ab} & -\frac{a}{3b}+\frac{b}{6a} & \frac{a^{2}+b^{2}}{3ab} & \frac{a}{6b}-\frac{b}{3a}\\ -\frac{a}{3b}+\frac{b}{6a} & -\frac{a^{2}+b^{2}}{6ab} & \frac{a}{6b}-\frac{b}{3a} & \frac{a^{2}+b^{2}}{3ab}\end{pmatrix} \end{align*}

For example, for element of width 1 and height 1, then $$a=\frac{1}{2},b=\frac{1}{2}$$ and the above becomes

$k_{i}= \begin{pmatrix} \frac{2}{3} & -\frac{1}{6} & -\frac{1}{3} & -\frac{1}{6}\\ -\frac{1}{6} & \frac{2}{3} & -\frac{1}{6} & -\frac{1}{3}\\ -\frac{1}{3} & -\frac{1}{6} & \frac{2}{3} & -\frac{1}{6}\\ -\frac{1}{6} & -\frac{1}{3} & -\frac{1}{6} & \frac{2}{3}\end{pmatrix}$

Hence

$I_{i}=-k_{i}\begin{Bmatrix} \tilde{u}_{1}\\ \tilde{u}_{2}\\ \tilde{u}_{3}\\ \tilde{u}_{4}\end{Bmatrix} -\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin{Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end{Bmatrix} f\left ( x,y\right ) \,dxdy$

Now the integration of the force vector is carried out.

\begin{align*} I_{i}^{f} & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin{Bmatrix} N_{1}\left ( x,y\right ) \\ N_{2}\left ( x,y\right ) \\ N_{3}\left ( x,y\right ) \\ N_{4}\left ( x,y\right ) \end{Bmatrix} f\left ( x,y\right ) \,dxdy\\ & =\int \limits _{x=-a}^{x=a}\int \limits _{y=-b}^{y=b}\begin{Bmatrix} \frac{1}{A}\left ( a-x\right ) \left ( b-y\right ) \\ \frac{1}{A}\left ( a+x\right ) \left ( b-y\right ) \\ \frac{1}{A}\left ( a+x\right ) \left ( b+y\right ) \\ \frac{1}{A}\left ( a-x\right ) \left ( b+y\right ) \end{Bmatrix} f\left ( x,y\right ) \,dxdy \end{align*}

In the above, the integration depends on the physical location of the element. We used a local coordinate system initially to obtain the shape functions. This has now to be transformed to the global coordinates system. Taking the center of each element as $$\left ( x_{0},y_{0}\right )$$ then we need to replace $$a\rightarrow x_{0}+a$$ and $$b\rightarrow y_{0}+b$$ everywhere in the above integral. We did not have to do this for ﬁnding the local stiﬀness matrix, since that did not depend on the physical location of the element (it was constant). But for the integration of the force function, we need to do this mapping. In the code, the center of each element is found, and the replacement is done.

#### 4.16.1 Integrating the force vector

4.16.1.1 $$f\left ( x,y\right ) =xy$$

This is normally done using Gaussian quadrature method. In this example, the integral is found oﬀ-line using Mathematica.

Evaluating the force vector requires integration over the element. This was done oﬀ-line using Mathematica. The result is a function of the coordinate of the center of the element. Then when running the code, this was evaluated for each element in the main loop.

##### 4.16.1.1 $$f\left ( x,y\right ) =xy$$

For $$f\left ( x,y\right ) =xy$$ the above gives a result as function of the center of the physical element.

$I_{i}^{f}=\int \limits _{x=\left ( x_{0}-a\right ) }^{x=x_{0}+a}\int \limits _{y=\left ( y_{0}-b\right ) }^{y=y_{0}+b}\begin{Bmatrix} \frac{1}{A}\left ( \left ( x_{0}-a\right ) -x\right ) \left ( y_{0}+b-y\right ) \\ \frac{1}{A}\left ( \left ( x_{0}-a\right ) +x\right ) \left ( y_{0}+b-y\right ) \\ \frac{1}{A}\left ( \left ( x_{0}-a\right ) +x\right ) \left ( y_{0}+b+y\right ) \\ \frac{1}{A}\left ( \left ( x_{0}-a\right ) -x\right ) \left ( y_{0}+b+y\right ) \end{Bmatrix} xy\,dxdy$

This was done using Mathematica

$I_{i}^{f}=\frac{1}{9\left ( x_{0}+a\right ) \left ( y_{0}+b\right ) }\begin{Bmatrix} a^{2}b^{2}\left ( a-3x_{0}\right ) \left ( b-3y_{0}\right ) \\ -ab^{2}\left ( a^{2}+3ax_{0}+6x_{0}^{2}\right ) \left ( b-3y_{0}\right ) \\ ab\left ( a^{2}+3ax_{0}+6x_{0}^{2}\right ) \left ( b^{2}+3by_{0}+6y_{0}^{2}\right ) \\ -ab^{2}\left ( a-3x_{0}\right ) \left ( b^{2}+3by_{0}+y_{0}^{2}\right ) \end{Bmatrix}$

#### 4.16.2 Case 1

4.16.2.1 case 2

Case $$f(x,y)=3\left ( \cos (4 \pi x)+\sin (3 \pi y)\right )$$

This was done using Mathematica

{(1/(96*Pi