From earlier lessons in this chapter, we assume you can: (i) rewrite an n-th order scalar ODE as a first-order vector model, (ii) construct the phase-variable (companion) form, and (iii) include multiple inputs/outputs by selecting appropriate output equations.

By the end of this lesson you will be able to:

• Convert common second- and third-order ODEs into a state-space model \( \dot{\mathbf{x}}=A\mathbf{x}+B\mathbf{u} \), \( \mathbf{y}=C\mathbf{x}+D\mathbf{u} \).
• Choose outputs that are not necessarily the original dependent variable (e.g., acceleration output).
• Prove equivalence: the state model implies the original ODE and vice versa (under consistent initialization).
• Implement the conversion in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

flowchart TD
  A["Start: given ODE and chosen outputs"] --> B["Normalize ODE (leading coeff = 1)"]
  B --> C["Pick states: x1=y, x2=y', ..., xn=y^(n-1)"]
  C --> D["Write x1' = x2, x2' = x3, ... , x(n-1)' = xn"]
  D --> E["Use ODE to write xn' = -a0 x1 - ... - a(n-1) xn + b*u"]
  E --> F["Assemble A and B matrices"]
  F --> G["Write outputs y = C x + D u"]
  G --> H["Check equivalence + consistent initial conditions"]
  H --> I["Implement in code (Python/C++/Java/MATLAB/Mathematica)"]
        

Consider a standard second-order forced ODE (mass–spring–damper): \( m\ddot{y} + c\dot{y} + ky = u(t) \), where \( y(t) \) is displacement and \( u(t) \) is force. We will define two outputs: \( y_1(t)=y(t) \) (position) and \( y_2(t)=\ddot{y}(t) \) (acceleration).

2.1 Normalize and choose states

Divide by \( m \) (assume \( m \neq 0 \)):

\[ \ddot{y} + \frac{c}{m}\dot{y} + \frac{k}{m}y = \frac{1}{m}u. \]

Choose the phase-variable states: \( x_1 \) and \( x_2 \) via

\[ x_1 \leftarrow y, \qquad x_2 \leftarrow \dot{y}. \]

2.2 Derive state equations

By definition, \( \dot{x}_1 = x_2 \). Also,

\[ \dot{x}_2 = \ddot{y} = -\frac{k}{m}x_1 - \frac{c}{m}x_2 + \frac{1}{m}u. \]

Therefore,

\[ \dot{\mathbf{x}} = \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = \underbrace{\begin{bmatrix} 0 & 1 \\ -\frac{k}{m} & -\frac{c}{m} \end{bmatrix}}_{A} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} + \underbrace{\begin{bmatrix} 0 \\ \frac{1}{m} \end{bmatrix}}_{B} u. \]

2.3 Output equations (position and acceleration)

The first output is simply \( y_1 = y = x_1 \). The second output uses the normalized dynamics:

\[ y_2 = \ddot{y} = -\frac{k}{m}x_1 - \frac{c}{m}x_2 + \frac{1}{m}u. \]

Thus,

\[ \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} = \underbrace{\begin{bmatrix} 1 & 0 \\ -\frac{k}{m} & -\frac{c}{m} \end{bmatrix}}_{C} \mathbf{x} + \underbrace{\begin{bmatrix} 0 \\ \frac{1}{m} \end{bmatrix}}_{D} u. \]

flowchart LR
  U["u(t)"] --> SUM["sum"]
  X1["x1 = y"] -->|"-k/m * x1"| SUM
  X2["x2 = y'"] -->|"-c/m * x2"| SUM
  SUM --> ACC["y2 = y''"]
  X1 --> POS["y1 = y"]
  SUM --> INT1["Integrate: x2"]
  INT1 --> X2
  X2 --> INT2["Integrate: x1"]
  INT2 --> X1
        

2.4 Equivalence proof

We prove that the state-space model implies the original ODE and that the ODE implies the state model under consistent state definitions.

Claim A (State model → ODE): If \( \dot{x}_1 = x_2 \) and \( \dot{x}_2 = -\frac{k}{m}x_1 - \frac{c}{m}x_2 + \frac{1}{m}u \), and if \( y \leftarrow x_1 \), then \( y \) satisfies \( m\ddot{y} + c\dot{y} + ky = u \).

Proof: From \( y \leftarrow x_1 \) and \( \dot{x}_1 = x_2 \), we have \( \dot{y} = x_2 \). Differentiating again gives \( \ddot{y} = \dot{x}_2 \). Substitute the second state equation:

\[ \ddot{y} = -\frac{k}{m}y - \frac{c}{m}\dot{y} + \frac{1}{m}u \;\;\Longrightarrow\;\; m\ddot{y} + c\dot{y} + ky = u. \]

This is exactly the given ODE. \(\square\)

Claim B (ODE → State model): If \( y \) satisfies \( m\ddot{y} + c\dot{y} + ky = u \), define \( x_1 \leftarrow y \), \( x_2 \leftarrow \dot{y} \). Then \( \mathbf{x} \) satisfies the derived state equations.

Proof: By definition, \( \dot{x}_1 = \dot{y} = x_2 \). Also, from the ODE (normalized), \( \ddot{y} = -\frac{k}{m}y - \frac{c}{m}\dot{y} + \frac{1}{m}u \). Hence, \( \dot{x}_2 = \ddot{y} = -\frac{k}{m}x_1 - \frac{c}{m}x_2 + \frac{1}{m}u \). \(\square\)

The two claims establish that the models are equivalent for all inputs \( u(t) \) for which the derivatives exist, and for consistent initial conditions: \( x_1(0)=y(0) \), \( x_2(0)=\dot{y}(0) \).

Consider the third-order ODE with two inputs \( u_1(t),u_2(t) \):

\[ y^{(3)} + a_2 \ddot{y} + a_1 \dot{y} + a_0 y = b_1 u_1 + b_2 u_2, \quad a_0,a_1,a_2,b_1,b_2 \in \mathbb{R}. \]

Choose states in phase-variable form:

\[ x_1 \leftarrow y,\qquad x_2 \leftarrow \dot{y},\qquad x_3 \leftarrow \ddot{y}. \]

3.1 State equations

The first two are definitional: \( \dot{x}_1 = x_2 \), \( \dot{x}_2 = x_3 \). Use the ODE to write \( \dot{x}_3 = y^{(3)} \):

\[ \dot{x}_3 = -a_0 x_1 - a_1 x_2 - a_2 x_3 + b_1 u_1 + b_2 u_2. \]

Stack them:

\[ \dot{\mathbf{x}} = \underbrace{\begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ -a_0 & -a_1 & -a_2 \end{bmatrix}}_{A} \mathbf{x} + \underbrace{\begin{bmatrix} 0 & 0\\ 0 & 0\\ b_1 & b_2 \end{bmatrix}}_{B} \begin{bmatrix} u_1\\u_2 \end{bmatrix}. \]

3.2 Two outputs

Let the measured outputs be \( y_1 \leftarrow y \) and \( y_2 \leftarrow \dot{y} \). Then

\[ \mathbf{y} = \begin{bmatrix} y_1\\y_2 \end{bmatrix} = \underbrace{\begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \end{bmatrix}}_{C}\mathbf{x} + \underbrace{\begin{bmatrix} 0 & 0\\ 0 & 0 \end{bmatrix}}_{D} \begin{bmatrix} u_1\\u_2 \end{bmatrix}. \]

3.3 Equivalence (short proof)

The same structure as Example 1 applies. The chain \( x_1=y \), \( \dot{x}_1=x_2 \), \( \dot{x}_2=x_3 \), and \( \dot{x}_3 \) from the ODE implies:

\[ y^{(3)} = \dot{x}_3 = -a_0 y - a_1 \dot{y} - a_2 \ddot{y} + b_1 u_1 + b_2 u_2, \]

which rearranges to the original ODE. Conversely, defining the states from a solution \( y(t) \) of the ODE yields the state model by differentiation and substitution.

Consider two coupled second-order ODEs with two inputs (already consistent with earlier coverage of vector/multi-channel modeling):

\[ \begin{aligned} \ddot{y}_1 + \alpha \dot{y}_1 + \beta y_1 + \gamma y_2 &= u_1,\\ \ddot{y}_2 + \alpha \dot{y}_2 + \gamma y_1 + \beta y_2 &= u_2. \end{aligned} \]

4.1 State choice

Define a 4-dimensional state vector:

\[ x_1 \leftarrow y_1,\quad x_2 \leftarrow \dot{y}_1,\quad x_3 \leftarrow y_2,\quad x_4 \leftarrow \dot{y}_2. \]

4.2 State equations

The definitional equations are \( \dot{x}_1=x_2 \), \( \dot{x}_3=x_4 \). From the ODEs:

\[ \dot{x}_2 = -\beta x_1 - \alpha x_2 - \gamma x_3 + u_1,\qquad \dot{x}_4 = -\gamma x_1 - \beta x_3 - \alpha x_4 + u_2. \]

In matrix form with \( \mathbf{x}=[x_1,x_2,x_3,x_4]^T \) and \( \mathbf{u}=[u_1,u_2]^T \):

\[ \dot{\mathbf{x}} = \underbrace{\begin{bmatrix} 0 & 1 & 0 & 0\\ -\beta & -\alpha & -\gamma & 0\\ 0 & 0 & 0 & 1\\ -\gamma & 0 & -\beta & -\alpha \end{bmatrix}}_{A} \mathbf{x} + \underbrace{\begin{bmatrix} 0 & 0\\ 1 & 0\\ 0 & 0\\ 0 & 1 \end{bmatrix}}_{B}\mathbf{u}. \]

4.3 Example outputs

If we measure both positions, set \( \mathbf{y}=[y_1,y_2]^T \):

\[ \mathbf{y} = \underbrace{\begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 \end{bmatrix}}_{C}\mathbf{x} + \underbrace{\begin{bmatrix} 0 & 0\\ 0 & 0 \end{bmatrix}}_{D}\mathbf{u}. \]

4.4 Equivalence remark

Each coupled ODE is recovered by differentiating the appropriate position state once to get its velocity state, twice to get acceleration, and then substituting the state equations for \( \dot{x}_2 \) and \( \dot{x}_4 \). The only requirement for exact equivalence is consistent initialization: \( x_1(0)=y_1(0) \), \( x_2(0)=\dot{y}_1(0) \), \( x_3(0)=y_2(0) \), \( x_4(0)=\dot{y}_2(0) \).

For the normalized scalar n-th order ODE \( y^{(n)} + a_{n-1}y^{(n-1)} + \cdots + a_1\dot{y} + a_0 y = \sum_{j=1}^{m} b_j u_j \), the phase-variable state choice \( x_1 \leftarrow y, x_2 \leftarrow \dot{y}, \ldots, x_n \leftarrow y^{(n-1)} \) yields:

\[ \dot{\mathbf{x}} = \underbrace{\begin{bmatrix} 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1\\ -a_0 & -a_1 & -a_2 & \cdots & -a_{n-1} \end{bmatrix}}_{A} \mathbf{x} + \underbrace{\begin{bmatrix} 0 & 0 & \cdots & 0\\ 0 & 0 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & 0\\ b_1 & b_2 & \cdots & b_m \end{bmatrix}}_{B} \mathbf{u}. \]

Output selection is problem-dependent but must be expressible using the states and inputs (as in the acceleration output in Example 1). The next lesson discusses scaling/normalization choices that improve numerical conditioning; here we focus on correctness and structural construction.

The function below constructs the phase-variable realization for a normalized n-th order ODE with multiple inputs. It also shows how to package the result as a state-space model using the widely used control library.

import numpy as np

def phase_variable_realization(a, b):
    """
    Build (A,B) for y^{(n)} + a[n-1] y^{(n-1)} + ... + a[0] y = sum_j b[j] u_j.
    Inputs:
      a: array-like of length n (a[0]=a0, ..., a[n-1]=a_{n-1})
      b: array-like of shape (m,) or (m,1) for m inputs
    Returns:
      A: (n,n), B: (n,m)
    """
    a = np.asarray(a, dtype=float).reshape(-1)
    n = a.size
    b = np.asarray(b, dtype=float).reshape(-1)
    m = b.size

    A = np.zeros((n, n), dtype=float)
    # shift structure
    for i in range(n - 1):
        A[i, i + 1] = 1.0
    # last row
    A[n - 1, :] = -a  # [-a0, -a1, ..., -a_{n-1}]

    B = np.zeros((n, m), dtype=float)
    B[n - 1, :] = b
    return A, B

# Example 2 parameters: y''' + a2 y'' + a1 y' + a0 y = b1 u1 + b2 u2
a0, a1, a2 = 2.0, 3.0, 1.0
b1, b2 = 5.0, -1.0
A, B = phase_variable_realization([a0, a1, a2], [b1, b2])

# Outputs y1=y, y2=y'
C = np.array([[1.0, 0.0, 0.0],
              [0.0, 1.0, 0.0]])
D = np.zeros((2, 2))

print("A=\n", A)
print("B=\n", B)
print("C=\n", C)
print("D=\n", D)

# Optional packaging with python-control (pip install control)
try:
    import control
    sys = control.ss(A, B, C, D)
    print(sys)
except Exception as e:
    print("python-control not available or not installed:", e)
      

Notes: (i) This implements the companion (phase-variable) construction introduced earlier. (ii) Simulation of \( \dot{\mathbf{x}}=A\mathbf{x}+B\mathbf{u} \) is deferred conceptually to Chapter 7; however, you can numerically integrate with basic ODE solvers if needed.

In C++, a practical approach is to use Eigen for matrices and implement a small constructor for the phase-variable form. (Eigen is a standard linear-algebra dependency in control-oriented C++ workflows.)

#include <iostream>
#include <vector>
#include <Eigen/Dense>

struct StateSpace {
  Eigen::MatrixXd A, B, C, D;
};

StateSpace phaseVariableRealization(const std::vector<double>& a,
                                    const std::vector<double>& b,
                                    const Eigen::MatrixXd& C,
                                    const Eigen::MatrixXd& D) {
  const int n = static_cast<int>(a.size());
  const int m = static_cast<int>(b.size());

  Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n, n);
  for (int i = 0; i < n - 1; ++i) A(i, i + 1) = 1.0;

  // Last row: [-a0, -a1, ..., -a_{n-1}]
  for (int j = 0; j < n; ++j) A(n - 1, j) = -a[j];

  Eigen::MatrixXd B = Eigen::MatrixXd::Zero(n, m);
  for (int j = 0; j < m; ++j) B(n - 1, j) = b[j];

  StateSpace sys;
  sys.A = A; sys.B = B; sys.C = C; sys.D = D;
  return sys;
}

int main() {
  // Example 2: y''' + a2 y'' + a1 y' + a0 y = b1 u1 + b2 u2
  std::vector<double> a = {2.0, 3.0, 1.0};   // a0, a1, a2
  std::vector<double> b = {5.0, -1.0};       // b1, b2

  Eigen::MatrixXd C(2,3);
  C << 1.0, 0.0, 0.0,
       0.0, 1.0, 0.0;
  Eigen::MatrixXd D = Eigen::MatrixXd::Zero(2,2);

  StateSpace sys = phaseVariableRealization(a, b, C, D);

  std::cout << "A=\n" << sys.A << "\n\n";
  std::cout << "B=\n" << sys.B << "\n\n";
  std::cout << "C=\n" << sys.C << "\n\n";
  std::cout << "D=\n" << sys.D << "\n\n";
  return 0;
}
      

If you later want time simulation before Chapter 7, a minimal numerical integrator (e.g., explicit Euler) can be coded directly, or you can use Boost.Odeint; however, the conceptual solution structure is introduced in Chapter 7.

Java has no universally dominant “control” package, so a common approach is to: (i) build (A,B,C,D) explicitly, and (ii) use a linear-algebra library such as Apache Commons Math or EJML. Below is a compact implementation using Apache Commons Math matrices.

import java.util.Arrays;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RealMatrix;

public class PhaseVariableSS {

  public static class StateSpace {
    public final RealMatrix A, B, C, D;
    public StateSpace(RealMatrix A, RealMatrix B, RealMatrix C, RealMatrix D) {
      this.A = A; this.B = B; this.C = C; this.D = D;
    }
  }

  // a: [a0, a1, ..., a_{n-1}], b: [b1, ..., bm]
  public static StateSpace phaseVariable(double[] a, double[] b, RealMatrix C, RealMatrix D) {
    int n = a.length;
    int m = b.length;

    RealMatrix A = MatrixUtils.createRealMatrix(n, n);
    for (int i = 0; i < n - 1; i++) {
      A.setEntry(i, i + 1, 1.0);
    }
    for (int j = 0; j < n; j++) {
      A.setEntry(n - 1, j, -a[j]);
    }

    RealMatrix B = MatrixUtils.createRealMatrix(n, m);
    for (int j = 0; j < m; j++) {
      B.setEntry(n - 1, j, b[j]);
    }

    return new StateSpace(A, B, C, D);
  }

  public static void main(String[] args) {
    // Example 2
    double[] a = new double[] {2.0, 3.0, 1.0};  // a0, a1, a2
    double[] b = new double[] {5.0, -1.0};      // b1, b2

    RealMatrix C = MatrixUtils.createRealMatrix(new double[][] {
      {1.0, 0.0, 0.0},
      {0.0, 1.0, 0.0}
    });
    RealMatrix D = MatrixUtils.createRealMatrix(new double[][] {
      {0.0, 0.0},
      {0.0, 0.0}
    });

    StateSpace sys = phaseVariable(a, b, C, D);

    System.out.println("A=\n" + Arrays.deepToString(sys.A.getData()));
    System.out.println("B=\n" + Arrays.deepToString(sys.B.getData()));
    System.out.println("C=\n" + Arrays.deepToString(sys.C.getData()));
    System.out.println("D=\n" + Arrays.deepToString(sys.D.getData()));
  }
}
      

Later in the course, these matrices will be used for controllability/observability tests and state-feedback design; for now, the focus is correct modeling.

MATLAB provides native state-space objects via ss. Below we implement Example 1 and Example 2.

9.1 MATLAB (Example 1)

% Example 1: m*ydd + c*yd + k*y = u
m = 2.0; c = 3.0; k = 5.0;

A = [0 1;
    -k/m  -c/m];
B = [0;
     1/m];

% Outputs: y1 = y, y2 = ydd = -(k/m)*x1 -(c/m)*x2 + (1/m)*u
C = [ 1    0;
     -k/m -c/m ];
D = [0;
     1/m];

sys = ss(A,B,C,D);
disp(sys)
      

9.2 MATLAB (Example 2)

% Example 2: y''' + a2 y'' + a1 y' + a0 y = b1 u1 + b2 u2
a0 = 2.0; a1 = 3.0; a2 = 1.0;
b1 = 5.0; b2 = -1.0;

A = [0 1 0;
     0 0 1;
    -a0 -a1 -a2];
B = [0 0;
     0 0;
     b1 b2];

% Outputs: y1 = y = x1, y2 = y' = x2
C = [1 0 0;
     0 1 0];
D = zeros(2,2);

sys = ss(A,B,C,D);
disp(sys)
      

9.3 Simulink construction (textual block recipe)

For Example 1, build a model using standard blocks (no images):

• Two integrators in cascade to generate \( x_2=\dot{y} \) and \( x_1=y \).
• A summation block computing \( \dot{x}_2 = -(k/m)x_1 -(c/m)x_2 + (1/m)u \).
• Gain blocks for \( -k/m \), \( -c/m \), and \( 1/m \).
• Output blocks: (i) \( y_1=x_1 \) and (ii) \( y_2=\dot{x}_2 \) (the summation output).

This structure matches the state equations exactly; it is a direct signal-flow implementation of the derived state model.

Mathematica supports state-space models through StateSpaceModel. Below is Example 1 in symbolic/numeric form.

(* Example 1: m y'' + c y' + k y = u *)
m = 2.0; c = 3.0; k = 5.0;

A = { {0, 1}, {-k/m, -c/m} };
B = { {0}, {1/m} };

(* Outputs: y1 = y = x1, y2 = y'' = -(k/m)x1 -(c/m)x2 + (1/m)u *)
C = { {1, 0}, {-k/m, -c/m} };
D = { {0}, {1/m} };

sys = StateSpaceModel[{A, B, C, D}];
sys
      

For Example 2 (two inputs), use \( B \) with two columns and \( D \) with appropriate dimensions. The construction is identical.

Problem 1 (Second-order with acceleration output): Convert \( \ddot{y} + 4\dot{y} + 3y = 2u \) to state space using \( x_1 \leftarrow y \), \( x_2 \leftarrow \dot{y} \). Let outputs be \( y_1 \leftarrow y \) and \( y_2 \leftarrow \ddot{y} \).

Solution:

\[ \dot{x}_1 = x_2,\qquad \dot{x}_2 = -3x_1 - 4x_2 + 2u. \]

\[ A=\begin{bmatrix}0&1\\-3&-4\end{bmatrix},\quad B=\begin{bmatrix}0\\2\end{bmatrix}. \]

\[ y_1 = x_1,\qquad y_2 = \ddot{y} = -3x_1 - 4x_2 + 2u \;\;\Longrightarrow\;\; C=\begin{bmatrix}1&0\\-3&-4\end{bmatrix},\quad D=\begin{bmatrix}0\\2\end{bmatrix}. \]

Problem 2 (Third-order, two inputs): Convert \( y^{(3)} + 2\ddot{y} + 5\dot{y} + y = 3u_1 - 4u_2 \) into state space in phase-variable form. Use outputs \( y_1 \leftarrow y \), \( y_2 \leftarrow \dot{y} \).

Solution:

\[ x_1 \leftarrow y,\; x_2 \leftarrow \dot{y},\; x_3 \leftarrow \ddot{y} \;\;\Longrightarrow\;\; \dot{x}_1=x_2,\; \dot{x}_2=x_3,\; \dot{x}_3=-x_1-5x_2-2x_3+3u_1-4u_2. \]

\[ A=\begin{bmatrix} 0&1&0\\ 0&0&1\\ -1&-5&-2 \end{bmatrix},\quad B=\begin{bmatrix} 0&0\\ 0&0\\ 3&-4 \end{bmatrix},\quad C=\begin{bmatrix} 1&0&0\\ 0&1&0 \end{bmatrix},\quad D=\begin{bmatrix} 0&0\\ 0&0 \end{bmatrix}. \]

Problem 3 (Coupled second-order): For the coupled system \( \ddot{y}_1 + y_1 + 2y_2 = u_1 \), \( \ddot{y}_2 + 2y_1 + y_2 = u_2 \), build a 4-state model with \( x_1 \leftarrow y_1 \), \( x_2 \leftarrow \dot{y}_1 \), \( x_3 \leftarrow y_2 \), \( x_4 \leftarrow \dot{y}_2 \). Take outputs \( y_1 \) and \( y_2 \).

Solution:

\[ \dot{x}_1=x_2,\quad \dot{x}_2=-x_1-2x_3+u_1,\quad \dot{x}_3=x_4,\quad \dot{x}_4=-2x_1-x_3+u_2. \]

\[ A=\begin{bmatrix} 0&1&0&0\\ -1&0&-2&0\\ 0&0&0&1\\ -2&0&-1&0 \end{bmatrix},\quad B=\begin{bmatrix} 0&0\\ 1&0\\ 0&0\\ 0&1 \end{bmatrix},\quad C=\begin{bmatrix} 1&0&0&0\\ 0&0&1&0 \end{bmatrix},\quad D=\mathbf{0}_{2\times 2}. \]

Problem 4 (Equivalence check by elimination): For Example 2 in this lesson, start from the state equations and eliminate the states to recover the original third-order ODE.

Solution:

Using \( x_1=y \), \( \dot{x}_1=x_2 \Rightarrow \dot{y}=x_2 \), \( \dot{x}_2=x_3 \Rightarrow \ddot{y}=x_3 \), and \( \dot{x}_3=-a_0x_1-a_1x_2-a_2x_3+b_1u_1+b_2u_2 \Rightarrow y^{(3)}=\dot{x}_3 \), substitute:

\[ y^{(3)} = -a_0 y - a_1 \dot{y} - a_2 \ddot{y} + b_1 u_1 + b_2 u_2 \;\;\Longrightarrow\;\; y^{(3)} + a_2 \ddot{y} + a_1 \dot{y} + a_0 y = b_1 u_1 + b_2 u_2. \]

Problem 5 (Output design): For Example 1, suppose you instead want outputs \( y_1 \leftarrow y \) and \( y_2 \leftarrow \dot{y} \). Write the new \( C \) and \( D \).

Solution:

\[ \mathbf{y}=\begin{bmatrix}y\\\dot{y}\end{bmatrix} = \underbrace{\begin{bmatrix}1&0\\0&1\end{bmatrix}}_{C}\mathbf{x} + \underbrace{\begin{bmatrix}0\\0\end{bmatrix}}_{D}u. \]

We converted representative higher-order ODEs into state-space models using phase-variable (companion) state choices and verified equivalence formally. The examples demonstrated how outputs can be selected as state components (position/velocity) or as derived variables (acceleration) expressible through \( C \) and \( D \). Implementations in multiple programming environments reinforce that the conversion is algorithmic and reproducible.