In the 18.06 ODE lecture, we found that the position ($x$) and velocity ($v$) of a mass bouncing on a spring without friction are described by the equations:
$$ \frac{d}{dt} \underbrace{\begin{pmatrix} x \\ v \end{pmatrix}}_\vec{x} = \underbrace{\begin{pmatrix} 0 & 1 \\ -k/m & 0 \end{pmatrix}}_A \vec{x} $$for (real) spring constant $k > 0$ and mass $m > 0$.
(a) Show that by simply rescaling the variables $\vec{x} = (x,v)$ to new variables $\vec{y} = (a,b) = D \vec{x}$ for some diagonal matrix $D$, then you can get equations: $$ \frac{d\vec{y}}{dt} = B \vec{y} $$ where $B$ is (real) anti-Hermitian (also called skew-Hermitian): $B^H = -B$. (What is the formula relating $B$ and $A$ via $D$? $B$ and $A$ are ______ to one another.)
(b) For any anti-Hermitian $B = -B^H$, the matrix $iB$ is a ________ matrix because $(iB)^H = \_\_\_$. It follows that the eigenvalues of $B$ are purely ________ and hence the solutions $\vec{y}(t)$ of $d\vec{y}/dt = B \vec{y}$ are ________ (oscillating/growing/decaying).
(c) For any anti-Hermitian $B = -B^H$, the matrix $Q = e^{Bt}$ is unitary (for any real $t$) because ________. (Note: from the series definition of $e^A$, you can immediately see that $(e^A)^H = e^{A^H}$.)
(d) Recall from class that $\Vert Q \vec{y} \Vert = \Vert \vec{y} \Vert$ for any unitary $Q$ (or more generally any $Q$ with orthonormal columns). From this, derive that the solution $$\vec{y}(t) = e^{Bt} \vec{y}(0)$$ for your $B$ from part (a) conserves energy: $\frac{1}{2} mv^2 + \frac{1}{2} kx^2$ (kinetic + potential energy) is a constant as a function of time $t$.
(a) By rescaling $\vec{x} = (x,v)$ to new variables $\vec{y} = (a,b) = D \vec{x}$ for some diagonal matrix $D = \begin{pmatrix} \alpha & 0 \\ 0 & \beta \end{pmatrix}$, we have $$ \frac{d\vec{y}}{dt} = \frac{d (D\vec{x})}{dt}= D\frac{d \vec{x}}{dt} = DA \vec{x} = DA D^{-1} \vec{y} = \begin{pmatrix} \alpha & 0 \\ 0 & \beta \end{pmatrix} \begin{pmatrix} 0 & 1 \\ -k/m & 0 \end{pmatrix} \begin{pmatrix} 1/\alpha & 0 \\ 0 & 1/\beta \end{pmatrix} \vec{y} = \begin{pmatrix} 0 & \frac{\alpha}{\beta} \\ -\frac{\beta k}{\alpha m} & 0 \end{pmatrix} \vec{y} = B \vec{y}. $$ Now, $$ B^H = \begin{pmatrix} 0 & \frac{\alpha}{\beta} \\ -\frac{\beta k}{\alpha m} & 0 \end{pmatrix}^H = \begin{pmatrix} 0 & -\frac{\bar{\beta} k}{\bar{\alpha} m} \\ \frac{\bar{\alpha}}{\bar{\beta}} & 0 \end{pmatrix}. $$ For $B$ to be anti-Hermitian, i.e. $B = -B^H$, we have $$ \frac{\alpha}{\beta} = \frac{\bar{\beta} k}{\bar{\alpha} m} \Longleftrightarrow \alpha^2 m = \beta^2 k. $$ Therefore, by taking $\boxed{D = \begin{pmatrix} \sqrt{k} & 0 \\ 0 & \sqrt{m} \end{pmatrix}}$, we obtain an anti-Hermitian matrix $\boxed{B = \begin{pmatrix} 0 & \sqrt{\frac{k}{m}} \\ -\sqrt{\frac{k}{m}} & 0 \end{pmatrix}}$.
Note that the formula relating $B$ and $A$ is $\boxed{B = D A D^{-1}}$, and hence $B$ and $A$ are similar to one another.
(b) For any anti-Hermitian $B = -B^H$, the matrix $iB$ is a Hermitian matrix because $(iB)^H = \boxed{-i B^H = -i(-B) = iB}$.
By the property of Hermitian matrices, the eigenvalues of $iB$ are real. It follows that the eigenvalues of $B$ are purely imaginary and hence the solutions $\vec{y}(t)$ of $d\vec{y}/dt = B \vec{y}$ are oscillating.
(c) For any anti-Hermitian $B = -B^H$, we have $$\boxed{QQ^H = e^{Bt} (e^{Bt})^{H} = e^{Bt} e^{(Bt)^H} = e^{Bt} e^{-Bt} = e^{(B-B)t} = I}.$$ Note that the second last equality above follows from the fact that $Bt$ and $-Bt$ commute. Therefore, $Q$ is unitary.
(d) Since $\Vert Q \vec{y} \Vert = \Vert \vec{y} \Vert$ for any unitary $Q$ and $e^{Bt}$ is unitary, for the solution $\vec{y}(t) = e^{Bt} \vec{y}(0)$ we have $$\Vert \vec{y}(t) \Vert^2 = \Vert e^{Bt} \vec{y}(0) \Vert^2 = \Vert \vec{y}(0) \Vert^2.$$ Since $\vec{y} = D\vec{x} = \begin{pmatrix} \sqrt{k}x \\ \sqrt{m}v \end{pmatrix}$, from the above equation we have $$\boxed{\frac{1}{2} mv(t)^2 + \frac{1}{2} kx(t)^2 = \frac{1}{2} \Vert \vec{y}(t) \Vert^2 = \frac{1}{2} \Vert \vec{y}(0) \Vert^2 = \frac{1}{2} mv(0)^2 + \frac{1}{2} kx(0)^2}.$$ This shows that $\frac{1}{2} mv^2 + \frac{1}{2} kx^2$ is a constant as a function of $t$, i.e. the solution conserves energy.
(Based on Strang, section 6.3, problem 18.)
(a) Write five terms of the infinite series for $e^{At}$. Take the derivative $d/dt$, and show that you have four terms of $Ae^{At}$.
Conclusion: $\frac{d}{dt} e^{At} = Ae^{At}$ as claimed in class, and $x(t) = e^{At}x(0)$ solves $dx/dt = Ax$.
(b) Using the same five terms of the infinite series for $f(A) = e^A$ (no $t$), write down four terms of $f(A+dA) - f(A)$ as a linear operator $f'(A)[dA]$ acting on $dA$ as in Lecture 10 (i.e. drop all terms higher than first order in the "infinitesimal" change $dA$).
$f'(A)[dA]$ is not simply $e^A dA$ (as it would be for scalars in 18.01) except for very special perturbations $dA$ that ______________.
(a) Note that $e^{At} = I + At + \frac{(At)^2}{2} + \frac{(At)^3}{3!} + \frac{(At)^4}{4!} + \dots$. Taking the derivative $d/dt$ of the first five terms, we have $$\boxed{0 + A + \frac{A^2 2t }{2} + \frac{A^3 3t^2 }{3!} + \frac{ A^4 4t^3}{4!} = A(I + At +\frac{(At)^2}{2} + \frac{(At)^3}{3!})}.$$ This shows that we have four terms of $Ae^{At}$.
(b) Using the first five terms of the infinite series for $f(A) = e^A = I + A + \frac{A^2}{2} + \frac{A^3}{3!} + \frac{A^4}{4!} + \dots$, four terms of $f(A+dA)-f(A)$ are $$(I + (A+dA) + \frac{(A+dA)^2}{2} + \frac{(A+dA)^3}{3!} + \frac{(A+dA)^4}{4!}) - (I + A + \frac{A^2}{2} + \frac{A^3}{3!} + \frac{A^4}{4!}) \\ \approx dA+\frac{A(dA)+(dA) A}{2} + \frac{A^2 (dA) + A (dA) A+ (dA) A^2}{6} + \frac{A^3(dA) + A^2 (dA) A + A (dA) A^2 + (dA) A^3}{24},$$ where all terms higher than first order in $dA$ are dropped. Therefore, we have $$\boxed{f'(A)[dA] = dA+\frac{A(dA)+(dA) A}{2} + \frac{A^2 (dA) + A (dA) A+ (dA) A^2}{6} + \frac{A^3(dA) + A^2 (dA) A + A (dA) A^2 + (dA) A^3}{24}}.$$
$f'(A)[dA]$ is not simply $e^A dA$ except for very special perturbations $dA$ that satisfy $\boxed{A (dA) = (dA) A}$. To see this, note that if the condition is satisifed, we have $$dA+\frac{A(dA)+(dA) A}{2} + \frac{A^2 (dA) + A (dA) A+ (dA) A^2}{6} + \frac{A^3(dA) + A^2 (dA) A + A (dA) A^2 + (dA) A^3}{24} \\ = dA+\frac{2A(dA)}{2} + \frac{3A^2 (dA)}{6} + \frac{4A^3(dA)}{24} = (I + A + \frac{A^2}{2} + \frac{A^3}{3!}) dA,$$ which suggests that $f'(A)[dA] = e^A dA$.
(a) Write down two familiar functions that solve the equation $d^2 y /dt^2 = -9y$. What is a solution $y(t)$ that starts with $y(0)=3$ and $y'(0)=0$?
(b) Find (by hand) the eigenvalues $\lambda_1$ and $\lambda_2$ and corresponding eigenvectors $x_1$ and $x_2$ of $$A = \begin{pmatrix} 1 & -5\\ 10 & -14\\ \end{pmatrix} .$$
(c) Consider the system of ODEs $\frac{d^2x}{dt^2} = Ax$ (note 2nd derivative!) with the $A$ from part (b). If $x(t) = c_1(t) x_1 + c_2(t) x_2$, what equations do $c_1(t)$ and $c_2(t)$ satisfy, and hence what solutions $c_1(t)$ and $c_2(t)$ are possible? (See part (a).)
(d) For the ODE of part (c), solve for $x(t)$ given the initial conditions $x(0) = (2,3)$ and $x'(0) = (1,4)$.
(e) Alternatively, if we let $y$ be the 4-component vector $y = (x, dx/dt)$, then $dy/dt = By$ for $B = \_\_\_\_$. Since this must ultimately give the same solutions $x(t)$ as part (c), what must be the eigenvalues of $B$?
(a) It can be observed that two functions that solve the equation are $\boxed{\cos(3t)}$ and $\boxed{\sin(3t)}$, or any linear combination thereof. (You could also use $e^{\pm i 3t} = \cos(3t)\pm i\sin(3t)$, but this is redundant, just a linear combination of sines and cosines.) Therefore, the general solution is $$ y(t) = a \cos(3t) + b \sin(3t).$$ Since $y(0) = 3$, we have $a = 3$. Since $y'(0) = 0$, we have $3b = 0 \Longleftrightarrow b = 0$. Therefore, we have $$ \boxed{y(t) = 3\cos(3t)}.$$
(b) To find the eigenvalues, we consider $$0 = \det(A-\lambda I) = \det\begin{pmatrix} 1-\lambda & -5\\ 10 & -14-\lambda\\ \end{pmatrix} = (1-\lambda)(-14-\lambda)-(-5)(10) = \lambda^2 + 13\lambda + 36 = (\lambda + 9)(\lambda + 4).$$ Therefore, the eigenvalues are $\boxed{\lambda_1 = -9}$ and $\boxed{\lambda_2 = -4}$.
For $\lambda_1 = -9$, we have $$ (A-\lambda_1 I) x_1 = 0 \Longleftrightarrow \begin{pmatrix} 10 & -5\\ 10 & -5\\ \end{pmatrix}x_1 = 0,$$ from which we see that an eigenvector is $\boxed{x_1 = \begin{pmatrix} 1 \\ 2 \end{pmatrix}}$.
For $\lambda_1 = -4$, we have $$ (A-\lambda_2 I) x_2 = 0 \Longleftrightarrow \begin{pmatrix} 5 & -5\\ 10 & -10\\ \end{pmatrix}x_2 = 0,$$ from which we see that an eigenvector is $\boxed{x_2 = \begin{pmatrix} 1 \\ 1 \end{pmatrix}}$.
(c) Note that $A x_1 = -9 x_1$ and $A x_2 = -4 x_2$. Therefore, if $x(t) = c_1(t) x_1 + c_2(t) x_2$, we have $$Ax = A(c_1(t) x_1 + c_2(t) x_2) = -9c_1(t)x_1 -4 c_2(t) x_2.$$ Also, we have $$\frac{d^2x}{dt^2} = \frac{d^2(c_1(t) x_1 + c_2(t) x_2)}{dt^2} = \frac{d^2(c_1(t))}{dt^2} x_1 + \frac{d^2(c_2(t))}{dt^2} x_2.$$ Therefore, for the system of ODEs $\frac{d^2x}{dt^2} = Ax$, $c_1(t)$ and $c_2(t)$ satisfy $\boxed{\frac{d^2(c_1(t))}{dt^2} = -9c_1(t)}$ and $\boxed{\frac{d^2(c_2(t))}{dt^2} = -4c_2(t)}$.
From (a), the possible solutions are $$\boxed{c_1(t) = a_1 \cos(3t) + b_1 \sin(3t)}$$ and $$\boxed{c_2(t) = a_2 \cos(2t) + b_2 \sin(2t)},$$ where $a_1, b_1, a_2, b_2$ are constants.
(d) From the initial condition $x(0) = (2,3)$, we have $c_1(0)x_1 + c_2(0)x_2 = (2,3)$ and hence $$ a_1 \begin{pmatrix} 1 \\ 2 \end{pmatrix} + a_2 \begin{pmatrix} 1 \\ 1 \end{pmatrix} = \begin{pmatrix} 2 \\ 3 \end{pmatrix} \Longleftrightarrow a_1 = a_2 = 1. $$
From the initial condition $x'(0) = (1,4)$, we have $c_1'(0)x_1 + c_2'(0)x_2 = (1,4)$ and hence $$ 3b_1 \begin{pmatrix} 1 \\ 2 \end{pmatrix} + 2b_2 \begin{pmatrix} 1 \\ 1 \end{pmatrix} = \begin{pmatrix} 1 \\ 4 \end{pmatrix} \Longleftrightarrow b_1 = 1, b_2 = -1. $$
Therefore, the solution is $$\boxed{x(t) = (\cos(3t) + \sin(3t)) \begin{pmatrix} 1 \\ 2 \end{pmatrix} + (\cos(2t) - \sin(2t)) \begin{pmatrix} 1 \\ 1 \end{pmatrix} = \begin{pmatrix} \cos(3t) + \sin(3t) +\cos(2t) - \sin(2t) \\ 2\cos(3t) + 2\sin(3t) + \cos(2t) - \sin(2t) \end{pmatrix}} .$$
(e) If $y = (x, dx/dt)$, then $$ \frac{dy}{dt} = \begin{pmatrix} \frac{dx}{dt} \\ \frac{d^2x}{dt^2} \end{pmatrix} = \begin{pmatrix} \frac{dx}{dt} \\ Ax \end{pmatrix} = \begin{pmatrix} 0 & I \\ A & 0 \end{pmatrix} \begin{pmatrix} x \\ \frac{dx}{dt} \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & -5 & 0 & 0 \\ 10 & -14 & 0 & 0 \end{pmatrix}\begin{pmatrix} x \\ \frac{dx}{dt} \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & -5 & 0 & 0 \\ 10 & -14 & 0 & 0 \end{pmatrix}y.$$ Therefore, we have $\boxed{B = \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & -5 & 0 & 0 \\ 10 & -14 & 0 & 0 \end{pmatrix}}$.
Since this must ultimately give the same solutions $x(t)$ as part (c), but has solutions $y(t)=e^{Bt} y(0)$, we see that the eigenvalues of $B$ must be purely imaginary and consist of the frequencies of the problem, i.e. ± the square roots of the eigenvalues of $A$, i.e. $$\boxed{\pm 3i \mbox{ and } \pm 2i}.$$
In chemistry, the stoichiometry matrix is often used to describe a set of $m$ reactions among $n$ different chemical "species" (e.g. H₂O, C₈H₁₀N₄O₂, and so on).
For example, consider the following 3 (fictitious) chemical reactions among 4 species, labeled $x_1, x_2, x_3, x_4$:
$$ x_1 + 2x_2 \longleftrightarrow 3 x_2 + 2x_4 \\ 2x_2 + 4x_3 \longleftrightarrow x_1 + x_4 \\ x_1 + 4x_3 \longleftrightarrow 5x_4 $$which would be represented by the stoichiometry matrix
$$ S = \begin{pmatrix} -1 & 3-2 & 0 & 2 \\ 1 & -2 & -4 & 1 \\ -1 & 0 & -4 & 5 \end{pmatrix} $$whose rows are the reactions and whose columns are the species. (Some authors use the transpose of this matrix instead.)
If we use a vector $\vec{x} \in \mathbb{R}^4$ to represent the concentrations of each of these four species, and a vector $\vec{r} \in \mathbb{R}^3$ to represent the rates of each reaction, then the rate of change of the concentrations is given by the system of ordinary differential equations (ODEs):
$$ \frac{d\vec{x}}{dt} = S^T \vec{r} $$(where the rates $\vec{r}$ are not generally constant: they may depend on the concentrations $\vec{x}$ in a complicated way … so you can't solve this just by multiplying the right-hand side by $t$).
(a) Describe (find a basis for) all possible reaction rates $\vec{r}$ for which $\frac{d\vec{x}}{dt} = 0$ (the system is in steady state).
(b) Certain linear combinations of the concentrations are conserved: there are some (time-independent) vectors $\vec{v} \in \mathbb{R}^4$ for which $\frac{d(\vec{v}^T \vec{x})}{dt} = 0$ for all possible rate vectors $\vec{r}$. If $\vec{v}$ doesn't depend on $t$, then $\frac{d(\vec{v}^T \vec{x})}{dt}$ is ________ times $\frac{d\vec{x}}{dt}$. These vectors $\vec{v}$ all lie within the ___-dimensional ________ space of $S$. Why?
(a) The system is in steady state if $\frac{d\vec{x}}{dt} = S^T \vec{r} = 0$. This occurs if and only if $\vec{r}$ is in the left nullspace of $S$: $\boxed{\vec{r}\in N(S^T)}$. To find a basis for $N(S^T)$, we apply Gaussian elimination: \begin{align} S^T = \begin{pmatrix} -1 & 1 & -1\\ 1 & -2 & 0 \\0 & -4 & -4 \\ 2 & 1& 5 \end{pmatrix} \rightarrow \begin{pmatrix} -1 & 1 & -1\\ 0 & -1 & -1 \\0 & -4 & -4 \\ 0 & 3 & 3 \end{pmatrix} \rightarrow \begin{pmatrix} 1 & -1 & 1\\ 0 & -1 & -1 \\0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \end{align}
We see that $S^T$ has two pivot columns, and one free column. We can deduce that $\mathrm{rank}(S)=\mathrm{rank}(S^T)=2$. Since $S^T$ has three columns, it will have a one dimensional nullspace. We can then find the special solution of $S^Tx = 0$, which is $x= \begin{pmatrix} -2 \\ -1 \\ 1 \end{pmatrix}$.
Therefore the reaction rates for which $\vec{r}$ for which $\frac{d\vec{x}}{dt} = 0$ are of the form \begin{align} \boxed{\vec{r} = r_0\begin{pmatrix} -2 \\ -1 \\ 1 \end{pmatrix}} \end{align} for all scalars $r_0\in\mathbb{R}$.
(b) We want to find $\vec{v}$ independent of time, for which $\frac{d(\vec{v}^T \vec{x})}{dt} = 0$ for all possible rate vectors $\vec{r}$. If $\vec{v}$ doesn't depend on $t$, then \begin{align} \frac{d(\vec{v}^T \vec{x})}{dt} = \boxed{\vec{v}^T}\frac{d\vec{x}}{dt} = \vec{v}^T S^T \vec{r} = 0. \end{align} We seek $\vec{v}$ for which this is true for all $\vec{r}$, and so $\vec{v}^T S^T = 0 \implies S\vec{v} = 0$. Therefore, the vectors $\vec{v}$ all lie within the $4-2 = \boxed{2}$-dimensional nullspace of $S$.
Although you were not required to do so, we can find a basis for such $\vec{v}$ by finding the nullspace of $S$: \begin{align} S = \begin{pmatrix} -1 & 1 & 0 & 2 \\ 1 & -2 & -4 & 1 \\ -1 & 0 & -4 & 5 \end{pmatrix} \rightarrow \begin{pmatrix} -1 & 1 & 0 & 2 \\ 0 & -1 & -4 & 3 \\ 0 & -1 & -4 & 3 \end{pmatrix} \rightarrow \begin{pmatrix} 1 & -1 & 0 & -2 \\ 0 & 1 & 4 & -3 \\ 0 & 0 & 0 & 0 \end{pmatrix} \rightarrow \begin{pmatrix} 1 & 0 & 4 & -5 \\ 0 & 1 & 4 & -3 \\ 0 & 0 & 0 & 0 \end{pmatrix} \end{align}
We can then calculate the special solutions of $Sx = 0$, which are \begin{align} \begin{pmatrix} -4 \\ -4 \\ 1 \\ 0 \end{pmatrix} \;\; \text{and} \;\; \begin{pmatrix} 5 \\ 3 \\ 0 \\ 1 \end{pmatrix} \end{align}
We conclude that $\frac{d(\vec{v}^T \vec{x})}{dt} = 0$ for $\vec{v}$ of the form: \begin{align} \boxed{\vec{v} = c_1 \begin{pmatrix} -4 \\ -4 \\ 1 \\ 0 \end{pmatrix} + c_2\begin{pmatrix} 5 \\ 3 \\ 0 \\ 1 \end{pmatrix}} \end{align}