<jupyter_start><jupyter_text>Chapter 3 Notes* Python functions are written with very narrow scope: they only address the individual problem-at-hand and _probably_ won't generalize.<jupyter_code>import sys
from importlib import reload
sys.path.append("..")
from sympy import *
init_printing(use_latex='mathjax')
s = symbols('s')
import Ch3.utilities as ch3_utils
reload(ch3_utils)<jupyter_output><empty_output><jupyter_text>Problem 3.1 Notes* The `sympy` package does algebra symbolically (check the `environment.yml` file)* The state transition matrix, STM, is represented by $\Phi(t)$* The helper function `computeResolvent` (see Ch3/utilities.py) does the heavy-lifting<jupyter_code>def helper3p1(A):
    res = ch3_utils.computeResolvent(A)
    pprint('Matrix A')
    pprint(A)
    pprint('has resolvent')
    pprint(res)<jupyter_output><empty_output><jupyter_text>a)<jupyter_code>A = Matrix([[-1, 0, 0],[1, -2, 0], [1, 2, -3]])
helper3p1(A)<jupyter_output>Matrix A
⎡-1  0   0 ⎤
⎢          ⎥
⎢1   -2  0 ⎥
⎢          ⎥
⎣1   2   -3⎦
has resolvent
⎡              1                                      ⎤
⎢            ─────                     0           0  ⎥
⎢            s + 1                                    ⎥
⎢                                                     ⎥
⎢           1       1                  1              ⎥
⎢       - ───── + ─────              ─────         0  ⎥
⎢         s + 2   s + 1              s + 2            ⎥
⎢                                                     ⎥
⎢    1         2         3          2       2      1  ⎥
⎢───────── - ───── + ─────────  - ───── + ─────  ─────⎥
⎣2⋅(s + 3)   s + 2   2⋅(s + 1)    s + 3   s + 2  s + 3⎦<jupyter_text>The inverse Laplace transform of the resolvent of A is:\begin{eqnarray}\Phi(t) = \begin{pmatrix}e^{-t} & 0 & 0 \\e^{-t} - e^{-2t} & e^{-2t} & 0 \\\frac{1}{2}e^{-3t} - 2 e^{-2t} + \frac{3}{2} e^{-t} & 2(e^{-2t} - e^{-3t}) & e^{-3t} \\\end{pmatrix}\end{eqnarray}  b)<jupyter_code>A = Matrix([[-1, 0, 0],[1, -1, 0], [0, 1, -1]])
helper3p1(A)<jupyter_output>Matrix A
⎡-1  0   0 ⎤
⎢          ⎥
⎢1   -1  0 ⎥
⎢          ⎥
⎣0   1   -1⎦
has resolvent
⎡   1                     ⎤
⎢ ─────       0        0  ⎥
⎢ s + 1                   ⎥
⎢                         ⎥
⎢   1         1           ⎥
⎢────────   ─────      0  ⎥
⎢       2   s + 1         ⎥
⎢(s + 1)                  ⎥
⎢                         ⎥
⎢   1         1        1  ⎥
⎢────────  ────────  ─────⎥
⎢       3         2  s + 1⎥
⎣(s + 1)   (s + 1)        ⎦<jupyter_text>The inverse Laplace transform of the resolvent of A is:\begin{eqnarray}\Phi(t) = \begin{pmatrix}e^{-t} & 0 & 0 \\t e^{-t} & e^{-t} & 0 \\\frac{t^2}{2} e^{-t} & t e^{-t} & e^{-t} \\\end{pmatrix}\end{eqnarray}The terms involving repeated roots require using the shift property of Laplace transforms, namely that\begin{eqnarray}\mathcal{L}^{-1} \bigl [ F(s-a) \bigr ] = e^{at} \mathcal{L}^{-1} \bigl [ F(s) \bigr ]\end{eqnarray}  c)<jupyter_code>A = Matrix([[-2, 1, 1],[1, -2, 1], [1, 1, -2]])
helper3p1(A)<jupyter_output>Matrix A
⎡-2  1   1 ⎤
⎢          ⎥
⎢1   -2  1 ⎥
⎢          ⎥
⎣1   1   -2⎦
has resolvent
⎡     2        1          1        1         1        1 ⎤
⎢ ───────── + ───   - ───────── + ───  - ───────── + ───⎥
⎢ 3⋅(s + 3)   3⋅s     3⋅(s + 3)   3⋅s    3⋅(s + 3)   3⋅s⎥
⎢                                                       ⎥
⎢      1        1        2        1          1        1 ⎥
⎢- ───────── + ───   ───────── + ───   - ───────── + ───⎥
⎢  3⋅(s + 3)   3⋅s   3⋅(s + 3)   3⋅s     3⋅(s + 3)   3⋅s⎥
⎢                                                       ⎥
⎢      1        1         1        1        2        1  ⎥
⎢- ───────── + ───  - ───────── + ───   ───────── + ─── ⎥
⎣  3⋅(s + 3)   3⋅s    3⋅(s + 3)   3⋅s   3⋅(s + 3)   3⋅s ⎦<jupyter_text>The inverse Laplace transform of the resolvent of A is:\begin{eqnarray}\Phi(t) = \frac{1}{3} \begin{pmatrix}2 e^{-3t} + 1 & 1 - e^{-3t} & 1 - e^{-3t} \\1 - e^{-3t} & 2 e^{-3t} + 1 & 1 - e^{-3t} \\1 - e^{-3t} & 1 - e^{-3t} & 2 e^{-3t} + 1 \\\end{pmatrix}\end{eqnarray}  Problem 3.2 Notes* The `sympy` package does algebra symbolically, as in previous problem* The helper function `firstCompanionForm` computes the first companion form described by Eqns. 3.88 and 3.94* The helper function `PFDforJordanForm` computes the partial fraction decomposition of an input transfer function.  This makes the Jordan form easier to calculate by Eqns. 3.108 and 3.116<jupyter_code>def helper3p2(num, den):
    pprint('Transfer function H(s):')
    pprint(factor(num)/factor(den))
    pprint('First Companion Form matrices A, B, C, D:')
    A, B, C, D = ch3_utils.firstCompanionForm(num, den)
    pprint(A)
    pprint(B)
    pprint(C)
    pprint(D)
    try:
        A, B, C, D = ch3_utils.jordanForm(num, den)
        pprint('Jordan (canonical) Form matrices A, B, C, D:')
        pprint(A)
        pprint(B)
        pprint(C)
        pprint(D)
    except NotImplementedError:
        pprint('Jordan form condition not implemented.')
        pprint('Here is the partial fraction decomposition to compute Jordan form manually.')
        num = factor(num)
        den = factor(den)
        pd = apart(num/den)
        pprint(pd)<jupyter_output><empty_output><jupyter_text>a)<jupyter_code>num = expand((s+2)*(s+4)).as_poly()
den = expand((s+1)*(s+3)*(s+5)).as_poly()
helper3p2(num, den)<jupyter_output>Transfer function H(s):
    (s + 2)⋅(s + 4)    
───────────────────────
(s + 1)⋅(s + 3)⋅(s + 5)
First Companion Form matrices A, B, C, D:
⎡ 0    1   0 ⎤
⎢            ⎥
⎢ 0    0   1 ⎥
⎢            ⎥
⎣-15  -23  -9⎦
⎡0⎤
⎢ ⎥
⎢0⎥
⎢ ⎥
⎣1⎦
[8  6  1]
[0]
3/(8*(s + 5)) + 1/(4*(s + 3)) + 3/(8*(s + 1))
Jordan (canonical) Form matrices A, B, C, D:
⎡-3  0   0 ⎤
⎢          ⎥
⎢0   -1  0 ⎥
⎢          ⎥
⎣0   0   -5⎦
⎡1⎤
⎢ ⎥
⎢1⎥
⎢ ⎥
⎣1⎦
[1/4  3/8  3/8]
[0]<jupyter_text>There are no repeated roots, so the Jordan form (Eqn. 3.108) is:\begin{eqnarray}A &=& \begin{pmatrix}-5 & 0 & 0 \\0 & -3 & 0 \\0 & 0 & -1 \\\end{pmatrix}, \B = \begin{pmatrix}1 & 1 & 1 \\\end{pmatrix}^T, \\C &=& \frac{1}{8} \begin{pmatrix}3 & 2 & 3\end{pmatrix}, \D = \begin{pmatrix}0\end{pmatrix}\end{eqnarray}_NOTE_: the form above is a bit different than the one returned by `jordanForm`.  A simple transformation to reorder the poles from smallest to largest was done offline because it was easiest.  b)<jupyter_code>num = expand((s+2)).as_poly()
den = expand(s*((s+1)**2 + 4)).as_poly()
helper3p2(num, den)<jupyter_output>Transfer function H(s):
     s + 2      
────────────────
  ⎛ 2          ⎞
s⋅⎝s  + 2⋅s + 5⎠
First Companion Form matrices A, B, C, D:
⎡0  1   0 ⎤
⎢         ⎥
⎢0  0   1 ⎥
⎢         ⎥
⎣0  -5  -2⎦
⎡0⎤
⎢ ⎥
⎢0⎥
⎢ ⎥
⎣1⎦
[2  1  0]
[0]
-(2*s - 1)/(5*(s**2 + 2*s + 5)) + 2/(5*s)
Jordan (canonical) Form matrices A, B, C, D:
⎡0   0     0  ⎤
⎢             ⎥
⎢0   0     1  ⎥
⎢             ⎥
⎣0  -1/5  -2/5⎦
⎡1⎤
⎢ ⎥
⎢0⎥
⎢ ⎥
⎣1⎦
[2/5  1/5  -2/5]
[0]<jupyter_text>c)<jupyter_code>num = expand((s+3)).as_poly()
den = expand((s+1)**2*(s+2)).as_poly()
helper3p2(num, den)<jupyter_output>Transfer function H(s):
     s + 3      
────────────────
       2        
(s + 1) ⋅(s + 2)
First Companion Form matrices A, B, C, D:
⎡0   1   0 ⎤
⎢          ⎥
⎢0   0   1 ⎥
⎢          ⎥
⎣-2  -5  -4⎦
⎡0⎤
⎢ ⎥
⎢0⎥
⎢ ⎥
⎣1⎦
[3  1  0]
[0]
Jordan form condition not implemented.
Here is the partial fraction decomposition to compute Jordan form manually.
  1       1        2    
───── - ───── + ────────
s + 2   s + 1          2
                (s + 1)<jupyter_text>We can use the partial fraction decomposition to do the Jordan form construction manually with $-1$ as a repeated root.The system can be decomposed into two subsystems. Subsystem 1: $H_1(s) = \frac{1}{s+2}$\begin{eqnarray}A_1 &=& \begin{pmatrix}-2\end{pmatrix}, B_1 = \begin{pmatrix}1\end{pmatrix}, \\C_1 &=& \begin{pmatrix}1\end{pmatrix}\end{eqnarray} Subsystem 2: $H_2(s) = \frac{1}{(s+1)^2} - \frac{1}{s+1}$Using Eqn 3.116, with $s_1 = -1$ and $r_{11}, r_{21} = -1, 1$:\begin{eqnarray}A_2 &=& \begin{pmatrix}-1 & 0 \\1 & -1 \\\end{pmatrix}, \B_2 = \begin{pmatrix}1 & 0 \\\end{pmatrix}^T, \\C_2 &=& \begin{pmatrix}-1 & 1\end{pmatrix}\end{eqnarray}Both systems are (again) combined block-diagonally:\begin{eqnarray}A &=& \begin{pmatrix}A_1 & 0_{1x2} \\0_{2x1} & A_2 \\\end{pmatrix}, \B = \begin{pmatrix}B_1^T & B_2^T \\\end{pmatrix}^T, \\C &=& \begin{pmatrix}C_1 & C_2\end{pmatrix}, \D = \begin{pmatrix}b_0\end{pmatrix}\end{eqnarray}which, for this problem, are:\begin{eqnarray}A &=& \begin{pmatrix}-2 & 0 & 0 \\0 & -1 & 0 \\0 & 1 & -1 \\\end{pmatrix}, \B = \begin{pmatrix}1 & 1 & 0 \\\end{pmatrix}^T, \\C &=& \begin{pmatrix}1 & -1 & 1\end{pmatrix}, \D = \begin{pmatrix}0\end{pmatrix}\end{eqnarray}  Problem 3.3 a)The matrices can be derived by just looking at the block-diagram as:\begin{eqnarray}A &=& \begin{pmatrix}-s_1 & 1 & 0 & \dots \\0 & -s_2 & 1 & 0 \dots \\0 & \dots & \ddots & \dots 1 \\0 & \dots & 0 & -s_k\end{pmatrix}, \B = \begin{pmatrix}0_{k-1} \\1\end{pmatrix}, \\C &=& \begin{pmatrix}C_1 & C_2 & \dots & C_k\end{pmatrix}, \D = \begin{pmatrix}0\end{pmatrix}\end{eqnarray} b)The dynamics of the system transformed to the normal form are _similar_ to the tandem form with $\mathbf{z}= T \mathbf{x}$ (where I have reversed the book's convention and am using $\mathbf{z}$ to refer to the normal form).  This amounts to finding a tranformation matrix, $T$ s.t. $T^{-1} \Lambda T = A$ from above, which is a Schur decomposition.  The upper-triangular Schur matrix, $U$, is simply $\Lambda$.  The Schur transformation matrix, $T$, is obtained by finding the eigenvectors of $A$ using $(A + s_iI)\mathbf{t}_i = \mathbf{0}$.  Since $A$ is upper-triangular, the eigenvalues are simply the diagonal elements.  To find the columns $\{\mathbf{t}_i \}_{1 \le i \le k}$ of $T$, it is necessary to solve linear equations of the form\begin{eqnarray}(s_i - s_1)t_i^1 + t_i^2 &=& 0 \\(s_i - s_2)t_i^2 + t_i^3 &=& 0 \\\vdots\end{eqnarray}for $0 \le i < k$.  **Note:  The superscripts above are not powers, they represent the numbered elements of the eigenvector $\mathbf{t}_i$**To find the first eigenvector, there is only a single equation, leaving $t_1^1$ free.  Choosing $t_1^1=1$ is fine (the application of $T^{-1}$ would rescale correctly regardless), so $\mathbf{t}_1 = (1 \ 0 \ \dots 0)^T$.  Setting the free parameter to 1 will be done subsequently for the same reason.To find the second eigenvector ($t_2^2 = 1$):\begin{eqnarray}(s_2 - s_1) t_2^1 + 1 &=& 0\end{eqnarray}which means $\mathbf{t}_2 = (-\frac{1}{s_2 - s_1} \ 1 \ 0 \ \dots 0)^T$To find the third eigenvector ($t_3^3 = 1$):\begin{eqnarray}(s_3 - s_1) t_3^1 + t_3^2 &=& 0 \\(s_3 - s_2) t_3^2 + 1 &=& 0\end{eqnarray}which means $\mathbf{t}_3 = (\frac{1}{(s_3 - s_2)(s_2 - s_1)} \ -\frac{1}{s_3 - s_2} \ 1 \ \dots 0)^T$To nail down the pattern, the fourth eigenvector is computed from the following set of linear equations:\begin{eqnarray}(s_4 - s_1) t_4^1 + t_4^2 &=& 0 \\(s_4 - s_2) t_4^2 + t_4^3 &=& 0 \\(s_4 - s_3) t_4^3 + 1 &=& 0\end{eqnarray}which means $\mathbf{t}_4 = (-\frac{1}{(s_4 - s_3)(s_4 - s_2)(s_4-s_1)} \ \frac{1}{(s_4 - s_3)(s_4-s_2)} \ -\frac{1}{s_4-s_3} \ 1 \dots 0)^T$Finally, putting it all together (and extending where necessary):\begin{eqnarray}T = \begin{pmatrix}1 & -\frac{1}{s_2 - s_1} & \frac{1}{(s_3 - s_2)(s_3-s_1)} & -\frac{1}{(s_4-s_3)(s_4 - s_2)(s_4-s_1)} & \dots & (-1)^{k-1} \prod_{i = 1}^k \frac{1}{(s_i - s_1)}\\0 & 1 & -\frac{1}{s_3 - s_2} & \frac{1}{(s_4 - s_3)(s_4-s_2)} & \dots & (-1)^{k-2} \prod_{i = 2}^k \frac{1}{(s_i - s_2)}\\0 & 0 & 1 & -\frac{1}{s_4 - s_3} & \dots & (-1)^{k-3} \prod_{i = 3}^k \frac{1}{(s_i - s_3)}\\0 & 0 & 0 & 1 & \dots & (-1)^{k-4} \prod_{i = 4}^k \frac{1}{(s_i - s_4)}\\ 0 & 0 & 0 & \ddots & \ddots & \vdots \\\vdots & \vdots & \vdots & \dots & \ddots & \vdots \\0 & 0 & 0 & \dots & \dots & 1\end{pmatrix}\end{eqnarray}<jupyter_code># perform a quick check of the above for a 3x3
S = numbered_symbols('s')
next(S)
s_array = [next(S) for i in range(3)]
A = Matrix([[-s_array[0], 1, 0], [0, -s_array[1], 1], [0, 0, -s_array[2]]])
T = Matrix([[1, -1/(s_array[1]-s_array[0]), 1/((s_array[2]-s_array[1])*(s_array[2]-s_array[0]))], 
            [0, 1, -1/(s_array[2]-s_array[1])], 
            [0, 0, 1]])
result = T**-1*A*T  # makes sure that T \Lambda T**-1 = A
for i in range(result.shape[0]):
    for j in range(result.shape[1]):
        result[i, j] = simplify(result[i, j])
pprint(result)<jupyter_output>⎡-s₁   0    0 ⎤
⎢             ⎥
⎢ 0   -s₂   0 ⎥
⎢             ⎥
⎣ 0    0   -s₃⎦<jupyter_text>The above check equals $\Lambda$ for the $3x3$ case, indicating the correct solution for $T$ was found.  Problem 3.4 From Eqn. 3.25:\begin{eqnarray}\frac{\partial \Phi(\tau, t)}{\partial \tau} = A(\tau) \Phi(\tau, t)\end{eqnarray}Using the fact that $\Phi(t, \tau) = \Phi^{-1}(\tau, t)$ and the hint given, it follows that:\begin{eqnarray}\frac{\partial \Phi(t, \tau)}{\partial \tau} &=& \frac{\partial \Phi^{-1}(\tau, t)}{\partial \tau} \\&=& -\Phi^{-1}(\tau, t) \frac{\partial \Phi(\tau, t)}{\partial \tau} \Phi^{-1}(\tau, t) \\&=& -\Phi^{-1}(\tau, t) A(\tau) \Phi(\tau, t) \Phi^{-1}(\tau, t) \\&=& -\Phi^{-1}(\tau, t) A(\tau) \\&=& -\Phi(t, \tau) A(\tau)\end{eqnarray}QED  Problem 3.5 The state space matrices can be quickly computed from the block-diagram:\begin{eqnarray}A &=& \begin{pmatrix}-a_1 & 1 & 0 & \dots \\-a_2 & 0 & 1 & 0 \dots \\\vdots & \dots & \ddots & \dots 1 \\-a_k & \dots & 0 & 0\end{pmatrix}, \B = \begin{pmatrix}0_{k-1} \\1\end{pmatrix}, \\C &=& \begin{pmatrix}q_0 - \sum_{i=1}^k q_i a_i & q_1 & \dots & q_{k-1}\end{pmatrix}, \D = \begin{pmatrix}q_k\end{pmatrix}\end{eqnarray}I used these matrices to find the transfer function via $H(s) = C (sI - A)^{-1} B + D$ (Eqn. 3.65) and the algorithm from Fig. 3.1 for computing the $E_i$ matrices in $\text{adj}(sI - A)$ (see Eqn. 3.67).  Some observations (I'm not going to work through all the details):* The form of $A$ means that $|sI - A|$, which is the denominator of $H(s)$, will be $s^k + a_1 s^{k-1} + a_2 s^{k-2} + \cdots + a_k$* The algorithm for computing $E_i$ recursively, combined with the sparseness of $B$, lead to incrementing of the index on the $q$ terms in lower powers of $s$ in the numerator of $H(s)$.The transfer function $H(s)$ is:\begin{eqnarray}H(s) = \frac{q_{k-1} s^{k-1} + [a_1 q_{k-1} + q_{k-2}]s^{k-2} + [(a_1 + a_2)q_{k-1} + q_{k-2} + q_{k-3}]s^{k-2} + \cdots}{s^k + a_1 s^{k-1} + a_2 s^{k-2} + \cdots + a_k}\end{eqnarray}The relationship between the $q_i$'s and the numerator coefficients of the standard transfer function representation (using $b_i$'s) is:\begin{eqnarray}b_0 &=& q_{k-1} \\b_1 &=& a_1 q_{k-1} + q_{k-2} \\b_2 &=& (a_1 + a_2)q_{k-1} + q_{k-2} + q_{k-3} \\\vdots \\b_j &=& q_{k-1} \sum_{i=1}^j a_i + \sum_{i=2}^{j+1} q_{k-i} \\\vdots \\b_k &=& q_{k-1} \sum_{i=1}^k a_i + \sum_{i=2}^k q_{k-i}\end{eqnarray}  Problem 3.6 a)The differential equations are:\begin{eqnarray}\ddot x + \frac{k^2}{Mr^2R} \dot x + \frac{mg}{M} \theta &=& \frac{k}{MRr} e \\\ddot \theta - \frac{M+m}{M \ell} g \dot \theta - \frac{k^2}{Mr^2R \ell} \dot x &=& -\frac{k}{MRr \ell} e\end{eqnarray}The 2 second-order ODEs above can be cast as 4 first-order ODEs with $x_1 \equiv x, x_2 \equiv \dot x, \theta_1 \equiv \theta, \theta_2 \equiv \dot \theta$:\begin{eqnarray}\dot x_1 &=& x_2 \\\dot x_2 &=& -\frac{k^2}{Mr^2R} x_2 - \frac{mg}{M} \theta_1 + \frac{k}{MRr} e \\\dot \theta_1 &=& \theta_2 \\\dot \theta_2 &=& \frac{M+m}{M \ell} g \dot \theta + \frac{k^2}{Mr^2R \ell} x_2 -\frac{k}{MRr \ell} e\end{eqnarray}With $\mathbf{x} \equiv [x_1 \ x_2 \ \theta_1 \ \theta_2]^T$, the state-space model (with output $\mathbf{y} = [x_1 \ \theta_1]^T$) can be written:\begin{eqnarray}\mathbf{\dot x} &=& A \mathbf{x} + B e \\\mathbf{y} &=& C \mathbf{x} + D e,\end{eqnarray}with\begin{eqnarray}A &=& \begin{pmatrix}0 & 1 & 0 & 0 \\0 & -\frac{k^2}{Mr^2R} & -\frac{mg}{M} & 0 \\0 & 0 & 0 & 1 \\0 & \frac{k^2}{Mr^2R \ell} & \frac{M+m}{M \ell} g & 0\end{pmatrix}, \B = \frac{k}{MRr}\begin{pmatrix}0 \\1 \\0 \\-\frac{1}{\ell}\end{pmatrix}, \\C &=& \begin{pmatrix}1 & 0 & 0 & 0 \\0 & 0 & 1 & 0\end{pmatrix}, \D = \begin{pmatrix}0 \\0\end{pmatrix}\end{eqnarray}  b)![title](problem3p6_smaller.jpg) c)<jupyter_code>a, b, c, d = symbols('a b c d')
# a = -k^2 / (Mr^2R)
# b = -mg / M
# c = k^2 / (Mr^2Rl)
# d = (M+m)g / (Ml)
A = Matrix([[0, 1, 0, 0],[0, a, b, 0], [0, 0, 0, 1], [0, c, d, 0]])
pprint(A)
res = simplify(ch3_utils.computeResolvent(A))
pprint(res)<jupyter_output>⎡0  1  0  0⎤
⎢          ⎥
⎢0  a  b  0⎥
⎢          ⎥
⎢0  0  0  1⎥
⎢          ⎥
⎣0  c  d  0⎦
⎡                     2                                                       
⎢1              -d + s                            b                           
⎢─  ───────────────────────────────  ───────────────────────────   ───────────
⎢s    ⎛         2                3⎞           2                3     ⎛        
⎢   s⋅⎝a⋅d - a⋅s  - b⋅c - d⋅s + s ⎠  a⋅d - a⋅s  - b⋅c - d⋅s + s    s⋅⎝a⋅d - a⋅
⎢                                                                             
⎢                    2                                                        
⎢               d - s                            b⋅s                          
⎢0   ────────────────────────────    ───────────────────────────     ─────────
⎢              2                3             2                3              
⎢    -a⋅d + a⋅s  + b⋅c + d⋅s - s     a⋅d - a⋅s  - b⋅c - d⋅s + s      a⋅d - a⋅s
⎢                                       [...]<jupyter_text>To get the STM, do a table-lookup for the inverse Laplace transform element-by-element (or evaluate numerically).  This would be tedious, so I am not going to do it. d)Split the matrices $C, D$ into two components: $C_x, D_x$ correspond to the first rows and $C_\theta$ the second rows.  The transfer function from $e$ to output $x$ is $H_x(s) = C_x(sI-A)^{-1}B + D_x$.  Similarly, the transfer function from $e$ to $\theta$ is $H_\theta(s) = C_\theta(sI-A)^{-1}+D_\theta$.  The resolvent $(sI-A)^{-1}$ was computed in part c) above.  Problem 3.7 a)From Problem 2.3's solution:\begin{eqnarray}A &=& \frac{1}{RC} \begin{pmatrix}-3 & 1 & 0 \\1 & -2 & 1 \\0 & 1 & -3\end{pmatrix}, \ B &=& \frac{1}{RC} \begin{pmatrix}2 \\0 \\0\end{pmatrix}, \E &=& \frac{1}{RC} \begin{pmatrix}0 \\0 \\2\end{pmatrix}\end{eqnarray} b)The solution to b) follows exactly the same as part c) above.<jupyter_code>R, C = symbols('R C')
A = Matrix([[-3/(R*C), 1/(R*C), 0],[1/(R*C), -2/(R*C), 1/(R*C)], [0, 1/(R*C), -3/(R*C)]])
pprint(A)
res = simplify(ch3_utils.computeResolvent(A))
pprint(res)<jupyter_output>⎡-3    1      ⎤
⎢───  ───   0 ⎥
⎢C⋅R  C⋅R     ⎥
⎢             ⎥
⎢ 1   -2    1 ⎥
⎢───  ───  ───⎥
⎢C⋅R  C⋅R  C⋅R⎥
⎢             ⎥
⎢      1   -3 ⎥
⎢ 0   ───  ───⎥
⎣     C⋅R  C⋅R⎦
⎡         ⎛ 2  2  2              ⎞                                            
⎢     C⋅R⋅⎝C ⋅R ⋅s  + 5⋅C⋅R⋅s + 5⎠                C⋅R                         
⎢─────────────────────────────────────  ───────────────────────  ─────────────
⎢ 3  3  3      2  2  2                  (C⋅R⋅s + 1)⋅(C⋅R⋅s + 4)   3  3  3     
⎢C ⋅R ⋅s  + 8⋅C ⋅R ⋅s  + 19⋅C⋅R⋅s + 12                           C ⋅R ⋅s  + 8⋅
⎢                                                                             
⎢                 C⋅R                       C⋅R⋅(C⋅R⋅s + 3)                   
⎢       ───────────────────────         ───────────────────────         ──────
⎢       (C⋅R⋅s + 1)⋅(C⋅R⋅s + 4)         (C⋅R⋅s + 1)⋅(C⋅R⋅s + 4)         (C⋅R⋅s
⎢                                                                             
⎢                                 [...]<jupyter_text>c)The transfer function from $e_0$ to $v_3$ is $H(s) = C_{v_3}(sI - A)^{-1}B + \frac{2}{RC}$, with $C_{v_3} = [0 \ 0 \ 1]^T$. Problem 3.8![3p8](kcl.jpg)With only resistors, capacitors, and voltage/current sources, any node in the network will look like the above.  By Kirchoff's loop law, the current through node $i$ must be the sum of the two currents due to the current source and $\frac{v}{R_i}$ ($v$ will be due to adjacent loops) minus the current drop across the capacitor.  This results in a differential equation $\dot v_i = \frac{1}{R_i C_i}(iR_i + v - v_i)$.  So, in the construction of the block submatrix in the full system matrix A, the entry at $(i,i)$ will be $-\frac{1}{R_i C_i}$, which is _negative_ clearly.  A consequence of Kirchoff's loop law is symmetry of the system matrix $A$: positive voltage difference across resistor $R_i$ into node $v_i$ is seen as a negative voltage difference into an adjacent node.From the problem above, the structure of the block due node $v_i$ becomes clear.  The magnitude of the element on the diagonal will be larger than the magnitude of the off-diagonal elements (terms due to current sources will be an additive bias not captured in the $A$ matrix).  This means that reduction to row-echelon form will result in an upper-triangular matrix with negative diagonal.  This proves that all eigenvalues lie on the _negative_ real axis, as required.  QED Problem 3.9 a)From Problem 2.5's solution:\begin{eqnarray}A &=& \begin{pmatrix}0 & 0 & 1 & 0 \\0 & 0 & 0 & 1 \\-\frac{K}{M} & \frac{K}{M} & -\frac{k^2}{r^2R} & 0 \\\frac{K}{M} & -\frac{K}{M} & 0 & -\frac{k^2}{r^2R}\end{pmatrix}, \B = \frac{k}{MRr}\begin{pmatrix}0 & 0 \\0 & 0 \\\frac{k}{Rr} & 0 \\0 & \frac{k}{Rr}\end{pmatrix}, \\C &=& \begin{pmatrix}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0\end{pmatrix}, \D = \begin{pmatrix}0 \\0\end{pmatrix}\end{eqnarray} b)<jupyter_code># with K, M, k, R, r = 40, 1, 2, 100, 0.02
A = Matrix([[0, 0, 1, 0], [0, 0, 0, 1], [-40, 40, -100, 0], [40, -40, 0, -100]])
res = computeResolvent(A)
C_x1 = Matrix([1, 0, 0, 0])
C_x2 = Matrix([0, 1, 0, 0])
B = Matrix([[0, 0], [0, 0], [1, 0], [0, 1]])
H_x1 = C_x1.T * res * B
H_x2 = C_x2.T * res * B
# H_x1
pprint(H_x1[0])
# H_x2
pprint(H_x2[0])
from numpy import sqrt
r1 = -100
r2 = 0
a, b, c = 1., 100., 80.
r3 = (-b + sqrt(b**2 - 4.*a*c)) / (2.*a)
r4 = (-b - sqrt(b**2 - 4.*a*c)) / (2.*a)
print('Roots of both transfer functions are nearly the same: {}, {}, {}, {}'.format(r1, r2, r3, r4))<jupyter_output>Roots of both transfer functions are nearly the same: -100, 0, -0.8065044950046243, -99.19349550499538
{'fraction_non_alphanumeric': 0.213540018981335, 'fraction_numerical': 0.05039092511411398, 'mean_word_length': 2.252204585537919, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 60, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '474670', 'license': 'MIT', 'n_tokens_mistral': 9825, 'n_tokens_neox': 9289, 'n_words': 2611, 'path': 'Ch3/FriedlandCh3.ipynb', 'repo_name': 'jwdinius/friedland-csd-solutions', 'size': '22127'}