How to solve an inhomogeneous two dimensional ODE system with initial conditions

Suppose we have \begin{align} \frac{\text{d}x}{\text{d}t} &= % -Ax+By+I\\ \frac{\text{d}y}{\text{d}t} &= % Cx-Dy\\ \end{align}

With $x(0)=0$ and $y(0)=y_0$. Let \begin{align} \vec{u}= \left( \begin{array}{c} x \\ y \end{array} \right), %%% \quad %%% M = \left( \begin{array}{cc} -A & B\\ C & -D \end{array} \right), %%% \quad %%% \vec{I}= \left( \begin{array}{c} I \\ 0 \end{array} \right), %%% \quad %%% \vec{u}(0)= \left( \begin{array}{c} 0 \\ y_0 \end{array} \right), \end{align}

Where all constants are $>0$. Therefore we can recast the system as \begin{align} \frac{\text{d}\vec{u}}{\text{d}t} = M\vec{u} + \vec{I} \end{align}

Let $\mu_+$ and $\mu_-$ be the eigenvalues of $M$. Let the homogeneous solution be given by $\vec{u}_h$ and the partial solution be given by $\vec{u}_p$, such that \begin{align} \frac{\text{d}\vec{u}_h}{\text{d}t} &= M\vec{u} \\ \end{align}

Therefore \begin{align} \vec{u_h}(t)= % \frac{y_0}{\mu_2-\mu_1} % \left[ %% \left(\begin{array}{c} B\\ A+\mu_1 \end{array}\right) e^{-\mu_1 t} %% - %% \left(\begin{array}{c} B\\ A+\mu_2 \end{array}\right) % e^{-\mu_2 t} % \right] \end{align}

where \begin{align} \mu_1 &= -\frac12 \left( \text{Tr}(M)+\sqrt{\text{Tr}(M)^2-4\det{M}} \right), % \\ % \mu_2 &= -\frac12\left(\text{Tr}(M)-\sqrt{\text{Tr}(M)^2-4\det{M}}\right), % \end{align}

Let \begin{align} U(t)= \frac{y_0}{\mu_2-\mu_1} % \left(\begin{array}{cc} Be^{-\mu_1 t} & Be^{-\mu_2 t} \\ \left(A+\mu_1\right)e^{-\mu_1 t} & \left(A+\mu_2\right)e^{-\mu_2 t} \end{array}\right) \end{align}

Notice that \begin{align} \frac{\text{d}U}{\text{d}t} = M U \end{align}

Therefore consider \begin{align} \vec{u}(t)=U(t)\vec{v}(t) \end{align}

Where $\vec{v}(0)=(1,1)^T$

giving \begin{align} \frac{\text{d}\vec{u}}{\text{d}t} % &= % \frac{\text{d}U}{\text{d}t}\vec{v}(t) + U\frac{\text{d}\vec{v}}{\text{d}t}\\ &= % MU\vec{v}(t) + U\frac{\text{d}\vec{v}}{\text{d}t} % = % M\vec{u} + \vec{I} % = % MU\vec{v}+\vec{I}\\ %% \Rightarrow \quad \frac{\text{d}\vec{v}}{\text{d}t} &= U^{-1}\vec{I}\\ %% \Rightarrow \quad \vec{v}(t) &= \int_0^t U^{-1}(s)\vec{I}\ \text{d}s\\ %% \Rightarrow \quad \vec{u}(t) &= U(t) \int_0^t U^{-1}(s)\vec{I}\ \text{d}s \end{align}

Is this the correct method? Have I got the correct initial conditions for $\vec{v}(t)$

• For a LTI system with a constant inhomogeneous input you will also get an additional constant in the output, which can also be described by a change in coordinates (changing coordinates can give you a homogeneous ODE). – Kwin van der Veen Aug 24 '16 at 13:17
• @fibonatic: changing coordinates sounds like a better way to solve this. A linear change in x and y? – Freeman Aug 24 '16 at 13:21
• Yes, a linear change in $x$ and $y$, such that $\hat{\vec{u}}=\vec{u}+\vec{w}$ is the linear transformation of $\vec{u}$ and satisfies $M\vec{u}+\vec{I}=M\hat{\vec{u}}$. – Kwin van der Veen Aug 24 '16 at 14:55
• You can actually do this for any polynomial inhomogeneity. – Kwin van der Veen Aug 29 '16 at 15:07

To solve:

$$\begin{cases} x'(t)=\text{I}+\text{B}\cdot y(t)-\text{A}\cdot x(t)\\ y'(t)=\text{C}\cdot x(t)-\text{D}\cdot y(t) \end{cases}$$

Use Laplace transform:

$$\begin{cases} s\text{X}(s)-x(0)=\frac{\text{I}}{s}+\text{B}\cdot\text{Y}(s)-\text{A}\cdot \text{X}(s)\\ s\text{Y}(s)-y(0)=\text{C}\cdot \text{X}(s)-\text{D}\cdot\text{Y}(s) \end{cases}$$

So, we can say:

$$\begin{cases} \text{X}(s)=\frac{\frac{\text{I}}{s}+\text{B}\cdot \text{Y}(s)+x(0)}{\text{A}+s}\\ \text{Y}(s)=\frac{\text{C}\cdot\text{X}(s)+y(0)}{\text{D}+s} \end{cases}$$

Now we can substiute them into each other to get:

• $$\text{X}(s)=\frac{s\text{B}y(0)+(\text{D}+s)(\text{I}+sx(0))}{s((\text{A}+s)(\text{D}+s)-\text{B}\text{C})}$$
• $$\text{Y}(s)=\frac{\text{C}(\text{I}+sx(0))+sy(0)(\text{A}+s)}{s((\text{A}+s)(\text{D}+s)-\text{B}\text{C})}$$
• Thank you for this! Give me a sec to work through it, it's been a couple of years since I used laplace transforms embarrassingly. – Freeman Aug 24 '16 at 13:24
• @Freeman You're welcome, I think this the most easy way of solving this LTI system. – Jan Aug 24 '16 at 13:24
• I agree, it is a nice method. I'm just attempting to transform it back, noting x(0)=0. – Freeman Aug 24 '16 at 13:28
• "The most easy way of solving this LTI system"... not sure. Because there is as much work for going back into the variable $t$ world ! In particular, due to the not simple discussion between simple and double poles. See my solution. I don't pretend it is shorter than yours ; I just wanted to say that Laplace Transform - that I have used a lot of times - may be more lengthy than it appears at first. – Jean Marie Aug 25 '16 at 15:52

First point: as it is an affine system, I transform it into a linear system by looking for constants $u$ and $v$ such that it is equivalent to \tag{1}\begin{align} \frac{\text{d}x}{\text{d}t} &= % -A(x+u)&+&B(y+v)&&\\ \frac{\text{d}y}{\text{d}t} &= % \ \ C(x+u)&-&D(y+v)&\\ \end{align}

A rapid computation shows that is possible using

$$\tag{2}u=\dfrac{CI}{BC-AD} \ \ \text{and} \ \ v=\dfrac{DI}{BC-AD}$$

(I have seen afterwards that you have done the same transformation in your solution).

Under the condition that the determinant $AD-BC$ of matrix $M$ is non zero, we are back to a linear system, that we will rename in the following way: $$Y'=MY$$ with $Y(0):=(u,y_0+v)^T.$

Second point: I directly use the matrix exponential form for the solution of this differential system under the form:

$$\tag{3}X=exp(t M).Y(0) \ \ \text{with} \ \ Y(0):=(u,y_0+v)^T$$

Here, we can rely on formulas that can be found in "Some Explicit Formulas for the Matrix Exponential" IEEE Trans. on Automatic Control, Vol 38, no.8, August 1993 (their corollary 3.3) making of course a distinction between the cases where the two eigenvalues $\lambda$ and $\mu$ of $M$ are equal or different, which corresponds to the fact that the discriminant of the characteristic equation of $M$, i.e.,

$$\tag{4}\delta:=(A-D)^2+4BC$$

is zero or non zero resp.:

• equal eigenvalues : $\tag{5}exp(tM)=e^{-(A+D)t/2}\begin{pmatrix}1-\frac{A-D}{2}t& Bt\\Ct&1+\frac{A-D}{2}t\end{pmatrix}.$

• different eigenvalues :

$\tag{5'}exp(tM)=e^{-(A+D)t/2}\begin{pmatrix}\gamma -\frac{A-D}{2}\sigma t& B \sigma t\\C \sigma t&\gamma +\frac{A-D}{2}\sigma t\end{pmatrix}$ using

$$\tag{6}\Delta:=\dfrac{1}{2}\sqrt{\delta}=\dfrac{1}{2}\sqrt{(A-D)^2+4BC}, \ \ \gamma:=\cosh{\Delta} \ \ \text{and} \ \ \sigma:=\dfrac{\sinh{\Delta}}{\Delta}.$$

The last straightforward step, that we will not do, amounts to plug either (5) or (5') into (3), using the values of $u$ and $v$ that are in relationship (2).

Remark: if $\delta<0$, the hyperbolic cosines/sines become naturally circular cosines/sines.

• Nice! Thanks for your input, it's nice to see different methods. – Freeman Aug 25 '16 at 12:58
• In fact, a good night after i wrote the first solution, I realized that there was a flaw in it (I couldn't write it into the $3 \times 3$ form I gave :-(...). Happily, I found another way to transform it at once into a linear system. For the rest, I think that you will agree that it is still a valuable solution where the two cases (equal and different eigenvalues) are neatly considered. – Jean Marie Aug 25 '16 at 15:26
• Thanks for the update!! Very cool :) I really appreciate the input. – Freeman Aug 26 '16 at 22:42

Suppose we have \begin{align} \frac{\text{d}x}{\text{d}t} &= % -Ax+By+I\\ \frac{\text{d}y}{\text{d}t} &= % Cx-Dy\\ \end{align}

With $x(0)=0$ and $y(0)=y_0$. Let us define $M$ as \begin{align} M = \left( \begin{array}{cc} -A & B\\ C & -D \end{array} \right) \end{align}

Now, by the transform on $x$ and $y$ such that \begin{align} x(t) &= X(t) + \frac{DI}{\det{M}}\\ y(t) &= Y(t) + \frac{CI}{\det{M}} \end{align}

the system becomes \begin{align} \frac{\text{d}X}{\text{d}t} &= % -AX+BY\\ \frac{\text{d}Y}{\text{d}t} &= % CX-DY\\ \end{align}

With \begin{align} X(0) = X_0 = -\frac{DI}{\det{M}},\quad Y(0) = Y_0 = y_0-\frac{CI}{\det{M}} \end{align}

Therefore, letting $\vec{U}(t) = (X,Y)^T$ \begin{align} \frac{\text{d}\vec{U}}{\text{d}t} &= M\vec{U} \\ \end{align}

Defining the absolute eigenvalues of $M$ as below, \begin{align} \mu_1 &= -\frac12 \left( \text{Tr}(M)+\sqrt{\text{Tr}(M)^2-4\det{M}} \right), % \\ % \mu_2 &= -\frac12\left(\text{Tr}(M)-\sqrt{\text{Tr}(M)^2-4\det{M}}\right), % \end{align}

Therefore \begin{align} \vec{U} &= % C_1 \left( \begin{array}{c} B \\ \mu_1+A \end{array} \right) e^{-\mu_1 t} % %% + %% C_2 \left( \begin{array}{c} B \\ \mu_2+A \end{array} \right) e^{-\mu_2 t} \end{align}

Therefore \begin{align} \left( \begin{array}{cc} B & B \\ \mu_1+A & \mu_2+A \end{array} \right) % \left( \begin{array}{c} C_1 \\ C_2 \end{array} \right) % &= % \left( \begin{array}{c} X_0 \\ Y_0 \end{array} \right)\\ %%%%%%%% \Rightarrow\quad \left( \begin{array}{c} C_1 \\ C_2 \end{array} \right) &= % \frac{1}{B(\mu_2-\mu_1)} % \left( \begin{array}{cc} \mu_2+A & -B \\ -\left(\mu_1+A\right) & B \end{array} \right) % \left( \begin{array}{c} X_0 \\ Y_0 \end{array} \right)\\ %%%%%%%%%%%%%%%% &= % \frac{1}{B(\mu_2-\mu_1)} % \left( \begin{array}{c} \left(\mu_2+A\right)X_0 - BY_0 \\ -\left(\mu_1+A\right)X_0 + BY_0 \end{array} \right) \end{align}

Therefore \begin{align} \vec{U} &= % \frac{1}{B(\mu_2-\mu_1)} \left[ % \left\{\left(\mu_2+A\right)X_0 - BY_0\right\} \left( \begin{array}{c} B \\ \mu_1+A \end{array} \right) e^{-\mu_1 t} % %% - %% \left\{\left(\mu_1+A\right)X_0 - BY_0\right\} \left( \begin{array}{c} B \\ \mu_2+A \end{array} \right) e^{-\mu_2 t} \right] \end{align}

Let us assume that $\mu_1\ll\mu_2$, giving, \begin{align} \vec{U} &\simeq % \frac{\left(\mu_2+A\right)X_0 - BY_0}{B(\mu_2-\mu_1)} \left( \begin{array}{c} B \\ \mu_1+A \end{array} \right) e^{-\mu_1 t} \end{align}

Therefore by the transformation performed earlier we find that \begin{align} \left( \begin{array}{c} x \\ y \end{array} \right) % &\simeq % \frac{\left(\mu_2+A\right)X_0 - BY_0}{B(\mu_2-\mu_1)} \left( \begin{array}{c} B \\ \mu_1+A \end{array} \right) e^{-\mu_1 t} %% +\frac{I}{\det{M}} % \left( \begin{array}{c} D \\ C \end{array} \right) \end{align}