Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part VI of the course APMA0340
Introduction to Linear Algebra with Mathematica

Glossary

Applications

Example 1: A balance between protein and mRNA synthesis and degredation/inactivation can be modeled to help in laboratory testing of processes involving biological information. A certain study uses an E. coli cell-free expression system, where the informational elements (mRNA and proteins) are studied outside of the body and their respective cells. Testing processes in these conditions is key to developing new biological engineering applications.

In this example, we examine endogenous mRNA inactivation to better understand gene expression dynamics. The protein used is deGFP, which is a highly translatable green fluorescent protein. The following differential equations describe the rate of change over time of the concentration of deGFP mRNA (m), dark deGFP (deGFPd), and fluorescent deGFP (deGFPf).

\begin{align*} \frac{\text d}{{\text d} t}\, m(t) &= -B\,m(t) , \\ \frac{\text d}{{\text d} t}\,\mbox{deGRPd}(t) &= a\,m(t) -k\,\mbox{deGFPd} (t), \\ \frac{\text d}{{\text d} t} \,\mbox{deGFPf} (t)&= k\,\mbox{deGFPd}(t) . \end{align*} 
clear;
syms B m(t) deGFPd(t) a k deGFPf(t);
eqns=[diff(m)==-B*m;diff(deGFPd)==a*m-k*deGFPd;diff(deGFPf)==k*deGFPd]
The various constants are defined as the protein production rate a, the mRNA inactivation rate B, and the maturation time of the protein 1/k (including folding, etc). The following initial conditions are given in the study as the starting point for the respective assays:
\[ m(0) = \frac{1}{4}, \qquad \mbox{deGFPd} (0) = \frac{1}{4} , \qquad \mbox{deGFPf}(0) =0 . \] 
cond=[deGFPf(0)==0;deGFPd(0)==0.25;m(0)==0.25]
We solve the above differential equations using these initial conditions to find equations for the concentration of each item in terms of time. 
soln=dsolve([eqns;cond]);
deGFPd_(t)=soln.deGFPd
\[ \mbox{deGFPd} (t) = e^{-kt}\,\frac{B+a-k}{4\left( B-k \right)} - \frac{a}{4\left( B-k \right)} \, e^{-Bt} . \] 
m_(t)=soln.m
\[ m(t) = \frac{1}{4}\, e^{-Bt} . \] 
deGFPf_(t)=soln.deGFPf
\[ \mbox{deGFPf}(t) = \frac{B+a}{4\,B} - \frac{B+a-k}{4\left( B-k \right)}\, e^{-kt} + \frac{ak}{4B\left( B-k \right)}\, e^{-Bt} . \]
The study gives the final solved equation as
\[ [\mbox{deGFP}_f (t) = \frac{a\,m_0}{B\left( B-k \right)} \left( k\,e^{-Bt} - k +B - B\, e^{-kt} \right) . \]
where m0 is the concentration of active mRNA when transcription is ended. This equation is algebraically equivalent to the result given by MATLAB (once some higher order terms are removed). This model follows the collected data closely as shown on page 5 of Shin and Noireaux (2010).
Graphs of mRNA and proteins
   ■
End of Example 1

Example 2: The differential equation governing the motion of a particle with electric charge q in an electromagnetic field is

\[ m\,\ddot{\bf r} = q \left( {\bf E} + {\bf v} \times {\bf B} \right) , \]
where m r, v, E, and B are, respectively, the mass of the particle, its position, its velocity, the electric field, and the magnetic field. This Lorentz force law (derived by Oliver Heaviside in 1889 a few years earlier than Hendrik Lorentz) can be written in component form as
\[ \begin{split} \ddot{x} (t) &= q \left( E_x + \dot{y}\, B_z - \dot{z}\, B_y \right) , \\ \ddot{y} (t) &= q \left( E_y + \dot{z}\, B_x - \dot{x}\, B_z \right) , \\ \ddot{z} (t) &= q \left( E_z + \dot{x}\, B_y - \dot{y}\, B_x \right) . \end{split} \]
We rewrite the above equation in vector form:
\[ \frac{{\text d}^2}{{\text d}t^2} \begin{bmatrix} x (t) \\ y(t) \\ z(t) \end{bmatrix} = q \begin{bmatrix} 0&B_z &-B_y \\ -B_z &0&B_x \\ B_y & -B_x &0 \end{bmatrix} \frac{{\text d}}{{\text d}t} \begin{bmatrix} x (t) \\ y(t) \\ z(t) \end{bmatrix} + q \begin{bmatrix} E_x \\ E_y \\ E_z \end{bmatrix} . \]
   ■
End of Example 2

Example 3: Suppose that a ball of mass m is released with some initial velocity from a particular point over a surface. Then its position (x,y) is determined from the follwoing initial avlue problem

\[ \begin{split} \ddot{x} &= 0 , \\ \ddot{y} &= g , \end{split} \qquad \begin{bmatrix} x(0) \\ y(0) \end{bmatrix} = \begin{bmatrix} x0 \\ y0 \end{bmatrix} , \qquad \begin{bmatrix} \dot{x}(0) \\ \dot{y}(0) \end{bmatrix} = \begin{bmatrix} vx0 \\ vy0 \end{bmatrix} . \]
When the ball reaches the surface, it bounces back with velocity of k value of the original hit.

 

(*Calculates the position and velocity of the ball after one bounce \ on a given curve*)
OneBounce[k_, ramp_][{t0_, x0_, xp0_, y0_, yp0_}] :=
Module[{sol, t1, x1, xp1, y1, yp1, gramp, gp},
sol = First[ NDSolve[{x''[t] == 0, x'[t0] == xp0, x[t0] == x0, y''[t] == -9.8,
y'[t0] == yp0, y[t0] == y0}, {x, y}, {t, t0, Infinity},
Method -> {"EventLocator", "Event" :> y[t] - ramp[x[t]]},
MaxStepSize -> 0.01]];
t1 = InterpolatingFunctionDomain[x /. sol][[1, -1]];
{x1, xp1, y1, yp1} = Reflection[k, ramp][{x[t1], x'[t1], y[t1], y'[t1]} /. sol];
Sow[{x[t] /. sol, t0 <= t <= t1}, "X"];
Sow[{y[t] /. sol, t0 <= t <= t1}, "Y"];
Sow[{x1, y1}, "Bounces"];
{t1, x1, xp1, y1, yp1}]

 

(*Calculates the position and velocity of the ball after reflecting \ off of the curve*)
Reflection[k_, ramp_][{x_, xp_, y_, yp_}] :=
Module[{gramp, gp, xpnew, ypnew}, gramp = -ramp'[x];
If[Not[NumberQ[gramp]], Print["Could not compute derivative"];
Throw[$Failed]];
gramp = {-ramp'[x], 1};
If[gramp.{xp, yp} == 0, Print["No reflection"];
Throw[$Failed]];
gp = {1, -1} Reverse[gramp];
{xpnew, ypnew} = (k/(gramp.gramp)) (gp gp.{xp, yp} - gramp gramp.{xp, yp});
{x, xpnew, y, ypnew}]

 

(*Calculates the path of the ball over a specific range of x bounce \ by bounce and plots it*)
BouncingBall[k_, ramp_, {x0_, y0_}] :=
Module[{data, end, bounces, xmin, xmax, ymin, ymax}, If[y0 < ramp[x0], Print["Start above the ramp"];
Return[$Failed]];
data = Reap[Catch[Sow[{x0, y0}, "Bounces"];
NestWhile[OneBounce[k, ramp], {0, x0, 0, y0, 0}, Function[1 - #1[[1]] / #2[[1]] > 0.01], 2, 25]], _, Rule];
end = data[[1, 1]];
data = Last[data];
bounces = ("Bounces" /. data);
xmax = Max[bounces[[All, 1]]];
xmin = Min[bounces[[All, 1]]];
ymax = Max[bounces[[All, 2]]];
ymin = Min[bounces[[All, 2]]];
Show[{Plot[ramp[x], {x, xmin, xmax}, PlotRange -> {{xmin, xmax}, {ymin, ymax}}, AspectRatio -> (ymax - ymin) / (xmax - xmin)],
ParametricPlot[ Evaluate[{Piecewise["X" /. data], Piecewise["Y" /. data]}], {t, 0, end}, PlotStyle -> RGBColor[1, 0, 0]]}]]
Now we apply these subroutings to different surfaces.
ramp[x_] := If[x < 1, 1 - x, 0];
BouncingBall[.7, ramp, {0, 1.25}]
(* 0.7 means the coefficient of rebouncing *)
circle[x_] := If[x < 1, Sqrt[1 - x^2], 0];
BouncingBall[.7, circle, {.1, 1.25}]
wavyramp[x_] := If[x < 1, 1 - x + .05 Cos[11 Pi*x], 0];
BouncingBall[.75, wavyramp, {0, 1.25}]
(*Small slope linear with high friction*)

lowslopelr[x_] := If[x < 1, 1 - 0.1 x, 0];
BouncingBall[.5, lowslopelr, {0, 1.25}]
(*Sine wave ramp*)
sineramp[x_] := If[x < 5, 1 - .5 Sin[x], 0];
BouncingBall[.75, sineramp, {0, 1.25}]
arccosramp[x_] := If[-1 < x < 1, ArcCos[x], 0];
BouncingBall[.5, arccosramp, {-0.1, 2.3}]
(*Logarithmic ramp*)
logramp[x_] := If[x < 2, 2 - Log[x + 1], 0];
BouncingBall[.65, logramp, {0, 3}]
(*Parabolic ramp*)
squareramp[x_] := If[x < 4, x^2, 0];
BouncingBall[.7, squareramp, {.5, 1}]
(*Cubic ramp*)
cuberamp[x_] := If[-1 < x < 5, 4 (x^3) - x, 0];
BouncingBall[.9, cuberamp, {-0.25, .5}]
   ■
End of Example 3

Example 4: Now we consider a bouncing ball on stairs.

(*standard stairs*)
c = .75;
sol = NDSolve[{y''[t] == -9.8, y[0] == 13.5, y'[0] == 5, a[0] == 13, WhenEvent[y[t] - a[t] == 0, y'[t] -> -c y'[t]], WhenEvent[Mod[t, 1], a[t] -> a[t] - 1]}, {y, a}, {t, 0, 8}, DiscreteVariables -> {a}];
Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 8}, PlotStyle -> Thick]
(*tall steps*)
c= .6;
sol = NDSolve[{y''[t] == -9.81, y[0] == 13.5, y'[0] == 5, a[0] == 13, WhenEvent[y[t] - a[t] == 0, y'[t] -> -c y'[t]], WhenEvent[Mod[t, 2], a[t] -> a[t] - 3]}, {y, a}, {t, 0, 8}, DiscreteVariables -> {a}];
Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 8}, PlotStyle -> Thick]
(*positive to negative steps*)
c = .5;
sol = NDSolve[{y''[t] == -9.81, y[0] == 10, y'[0] == 5, a[0] == 4, WhenEvent[y[t] == a[t], y'[t] -> -c y'[t]], WhenEvent[Mod[t, 2], a[t] -> 1 - a[t]], WhenEvent[ y[t] < a[t] , y''[t] = -y'[t]]}, {y, a}, {t, 0, 4}, DiscreteVariables -> {a}];
Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 8}, PlotStyle -> {Thick, Thick}]
   standard stairs    tall steps    positive to negative steps
(*corner to corner bounces (c)*)
c= .7;
sol = NDSolve[{y''[t] == -9.81, y[0] == 10, y'[0] == 5, a[0] == 3, WhenEvent[y[t] == a[t], y'[t] -> -c y'[t]], WhenEvent[a[t] == y[t], a[t] -> 1 - (a[t])^2], WhenEvent[ y[t] < a[t] , y''[t] = -y'[t]]}, {y, a}, {t, 0, 10}, DiscreteVariables -> {a}];
Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 8}, PlotStyle -> {Thick, Thick}]

(*time based steps*)
c = .5;
sol = NDSolve[{y''[t] == -9.81, y[0] == 10, y'[0] == 5, a[0] == -1, WhenEvent[y[t] == a[t], y'[t] -> -c y'[t]], WhenEvent[t == -a[t], a[t] -> a[t] - 3], WhenEvent[ y[t] < a[t] , y''[t] = -y'[t]]}, {y, a}, {t, 0, 5}, DiscreteVariables -> {a}];
Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 8}, PlotStyle -> {Thick, Thick}]

c = 0.8;
sol = NDSolve[{y''[t] == -9.81, y[0] == 10, y'[0] == 6, a[0] == -1, WhenEvent[y[t] == a[t], y'[t] -> -c y'[t]], WhenEvent[t == -a[t], a[t] -> a[t] - 3], WhenEvent[ y[t] < a[t] , y''[t] = -y'[t]]}, {y, a}, {t, 0, 5}, DiscreteVariables -> {a}];
Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 8}, PlotStyle -> {Thick, Thick}]
   corner to corner bounces    time based steps    elastic bounce
   ■
End of Example 4

Controllability of Linear Systems

The concept of controllability analysis plays a crucial role in many regulation problems, such as the stabilization of unstable systems using feedback, tracking problems, obtaining optimal control strategies, or, simply prescribing an input that has a desired effect on the state. It considers whether a system can be steered from an arbitrary initial state to any desired final state using appropriate control inputs within a given time frame. Controllability in linear systems is governed by the qualities of the system’s state matrix and control matrix. A system is considered controllable if the reachable states can span its state space under the effect of control inputs. This feature is critical in engineering and control theory because it underpins the design and execution of successful control techniques for a wide range of applications, including robotics and aerospace, as well as economics and chemical processes.

The system . \( \displaystyle \quad \dot{\bf x}(t) = \mathbf{A}\,\mathbf{x}(t) + \mathbf{B}\,\mathbf{u}(t) \quad \) with initial condition x(t₀) = x₀) is said to be controllable in the interval [t₀ , t₁] if for every x₀, x₁ ∈ ℝn × 1 there exists a control input u. ∈ 𝔏²([t₀ , t₁]; ℝᵐ) such that the corresponding solution starting from x(t₀) = x₀ also satisfies x(t₁) = x₁.
Using the variation of parameter method the solution of the initial value problem
\begin{equation} \label{EqA.1} \dot{\bf x}(t) = \mathbf{A}\,\mathbf{x}(t) + \mathbf{B}\,\mathbf{u}(t) , \qquad \mathbf{x}(t_0 ) = \mathbf{x}_0 \end{equation}
can be written in the form;
\begin{equation} \label{EqA.2} \mathbf{x}(t) = e^{\mathbf{A}\left( t - t_0 \right)} \,\mathbf{x}_0 + \int_{t_0}^t \, e^{\mathbf{A}\left( t - \tau \right)} \,\mathbf{B}\,\mathbf{u}(\tau )\,{\treext d}\tau . \end{equation}
It follows that the system \eqref{EqA.1} is controllable if and only if there exists a control function u ∈ 𝔏²([t₀ , t₁]; ℝᵐ) such that
\[ \mathbf{x}\left( t_1\right) = \mathbf{x}_1 = e^{\mathbf{A}\left( t_1 - t_0 \right)} \,\mathbf{x}_0 + \int_{t_0}^{t_1} \, e^{\mathbf{A}\left( t_1 - \tau \right)} \,\mathbf{B}\,\mathbf{u}(\tau )\,{\text d}\tau . \]
That is,
\[ \mathbf{x}_1 - e^{\mathbf{A}\left( t_1 - t_0 \right)} \,\mathbf{x}_0 = \int_{t_0}^{t_1} \, e^{\mathbf{A}\left( t_1 - \tau \right)} \,\mathbf{B}\,\mathbf{u}(\tau )\,{\text d}\tau . \]
If we set \( \displaystyle \quad \mathbf{w} = e^{\mathbf{A}\left( t_1 - \tau \right)} \,\mathbf{B}\,\mathbf{u}(\tau )\,{\text d}\tau ,\quad \) system \eqref{EqA.1} is controllable if and only if ∀w ∈ ℝn × 1 there exists a control function u ∈ 𝔏²([t₀ , t₁]; ℝᵐ) such that
\[ \mathbf{w} = \int_{t_0}^{t_1} \, e^{\mathbf{A}\left( t_1 - \tau \right)} \,\mathbf{B}\,\mathbf{u}(\tau )\,{\text d}\tau . \]
Let us define the linear operator    C : 𝔏²([t₀ , t₁]; ℝᵐ) ⇾ ℝn × 1 by
\[ C\,\mathbf{u} = \int_{t_0}^{t_1} \, e^{\mathbf{A}\left( t_1 - \tau \right)} \,\mathbf{B}\,\mathbf{u}(\tau )\,{\text d}\tau . \]
Observe that . C is a bounded linear operator and system \eqref{EqA.1} is controllable if and only if Cu = w has a solution for every w ∈ ℝn × 1. That is, controllability of system \eqref{EqA.1} is equivalent to the surjectivity of the operator C. The operator C defines its adjoint C✶ : ℝn × 1 → 𝔏²([t₀ , t₁]; ℝᵐ) in the following way:
\begin{align*} \left\langle C\,\mathbf{u} , v \right\rangle &= \left\langle \int_{t_0}^{t_1} \, e^{\mathbf{A}\left( t_1 - \tau \right)} \,\mathbf{B}\,\mathbf{u}(\tau )\,{\text d}\tau , v \right\rangle \\ &= \int_{t_0}^{t_1} \, \left\langle e^{\mathbf{A}\left( t_1 - \tau \right)} \,\mathbf{B}\,\mathbf{u}(\tau )\,{\text d}\tau , v \right\rangle \\ &= \int_{t_0}^{t_1} \, \left\langle \mathbf{u}(\tau ) , \mathbf{B}^{\mathrm T}\,e^{\mathbf{A}^{\mathrm T}\left( t_1 - \tau \right)} \,\mathbf{v} \right\rangle\, {\text d}\tau \\ &= \left\langle \mathbf{u} , C^{\ast}\,\mathbf{v} \right\rangle . \end{align*}
That is,
\[ \left( C^{\ast} \,\mathbf{v} \right) (\tau ) = \mathbf{B}^{\mathrm T}\, e^{\mathbf{A}^{\mathrm T}\left( t_1 - \tau \right)} \,\mathbf{v} . \]
The following theorem explains the relation between controllability of the system \eqref{EqA.1} with the operators C and C✶.

Theorem 1: System \eqref{EqA.1} is controllable if and only if one of the following conditions holds.
  1. The operator C is onto (surjective).
  2. The operator C✶ is onto.
  3. The Gramian matrix \[ \mathbf{W}\left( t_0 , t_1 \right) = C\,C^{\ast} = \int_{t_0}^{t_1} \, e^{\mathbf{A}\left( t_1 - \tau \right)} \,\mathbf{B}\,\mathbf{B}^{\mathrm T} \,e^{\mathbf{A}^{\mathrm T} \left( t_1 - \tau \right)} \,{\text d}\tau \] is non-singular.
Thus, we have observed that system \eqref{EqA.1} is controllable if and only if there exists a control function u ∈ 𝔏²([t₀ , t₁]; ℝᵐ) such that
\[ \mathbf{w} = C\,\mathbf{u} = \int_{t_0}^{t_1} \, e^{\mathbf{A}\left( t_1 - \tau \right)} \,\mathbf{B}\,\mathbf{u}(\tau )\,{\text d}\tau . \]

Example 5:    ■

End of Example 5
By Cayley–Hamilton theorem, we can write the above equation as
\begin{align*} \mathbf{w} &= \int_{t_0}^{t_1} \, \left[ \mathcal{P}_0 (\tau )\,\mathbf{I}_n + \mathcal{P}_1 (\tau )\,\mathbf{A} + \cdots + \mathcal{P}_{n-1} (\tau )\,\mathbf{A}^{n-1} \right] \,\mathbf{B}\,\mathbf{u}(\tau )\,{\text d}\tau \\ &= \int_{t_0}^{t_1} \,\left[ \mathbf{B} \ \mathbf{A}\,\mathbf{B} \ \cdots \ \mathbf{A}^{n-1} \mathbf{B} \right] \begin{pmatrix} \mathcal{P}_0 (\tau ) \\ \mathcal{P}_1 (\tau ) \\ \vdots \\ \mathcal{P}_{n-1} (\tau ) \end{pmatrix} {\text d}\tau , \end{align*}
where 𝒫i(τ), i = 0, 1, … n−1, are polynomial functions that appear in the exponential expansion of \( \displaystyle \quad e^{\mathbf{A}\left( t_1 - \tau \right)} . \quad \) Thus, we can say that system \eqref{EqA.1} is controllable if and only if the range (image) of [B A B … An−1B] spans the who;e space ℝn × 1; that is, if and only if Rank([B A B … An−1B]) = n. This result is proposed by the Hungarian-American electrical engineer and math- ematician Rudolf Emil Kálmán (1930–2016) and is known as Kálmán’s Rank Condition for controllability.

Theorem 2 (Kálmán, 1960): A linear time-invariant system \( \displaystyle \quad \dot{\bf x}(t) = \mathbf{A}\,\mathbf{x}(t) + \mathbf{B}\,\mathbf{u}(t) \quad \) is completely state controllable if and only if its controllability matrix
\[ \mathbf{Q} = \left[ \mathbf{B}\ \mathbf{A}\,\mathbf{B}\ \cdots \ \mathbf{A}^{n-1}\mathbf{B} \right] \]
has full rank equal to n, the dimension of the state vector x.

Example 6:    ■

End of Example 6
Now, suppose that there exists a row vector v ∈ ℝⁿ such that v A = λv and v B = 0. Then observe that
\[ \mathbf{v} \left[ \mathbf{B} \ \mathbf{A}\,\mathbf{B} \ \mathbf{A}^2 \mathbf{B} \ \cdots \ \mathbf{A}^{n-1}\,\mathbf{B} \right] = 0 , \]
and hence Rank([B A B ⋯ An-1B]) < n, which implies that system \eqref{EqA.1} is not controllable. Thus, for the controllability of system \eqref{EqA.1}, no row vector . v ∈ ℝⁿ with v A = λv should be orthogonal to the columns of B. This method is known after the mathematicians https://en.wikipedia.org/wiki/Vasile_M._Popov">Vasile M. Popov (1928-), https://en.wikipedia.org/wiki/Vitold_Belevitch">Vitold Belevitch (1921–1999) and Jr. Malo L.J. Hautus (1940–).

Theorem 3: System \eqref{EqA.1} is controllable if and only if one of the following conditions holds.
  1. PBH Rank Condition: . Rank[s Iₙ − A, B] = n,     ∀ s ∈ ℂ.
  2. PBH Eigenvector Condition: the relationship v A = λv implies v B ≠ 0, where v is a left eigenvector of A associated with the eigenvalue . λ.

Example 7: the satellite system

page 280 at in A Course in Linear Algebra by George &    ■
End of Example 7

 

 

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations