用MATLAB解系统方程：线性齐次状态方程的解； 线性非齐次状态方程的解； 连续系统状态方程的离散化。

1. 线性齐次状态方程的解

假设我们有一个线性齐次状态方程：

$\dot{x} = Ax$

其中，$x$ 是一个 $n$ 维向量，$A$ 是一个 $n \times n$ 的矩阵。我们可以通过求解矩阵 $A$ 的特征值和特征向量，得到方程的通解：

$x(t) = c_1 e^{\lambda_1 t} v_1 + c_2 e^{\lambda_2 t} v_2 + \cdots + c_n e^{\lambda_n t} v_n$

其中，$\lambda_1, \lambda_2, \cdots, \lambda_n$ 是 $A$ 的特征值，$v_1, v_2, \cdots, v_n$ 是相应的特征向量，$c_1, c_2, \cdots, c_n$ 是待定的常数。在 MATLAB 中，我们可以使用 eig 函数求解特征值和特征向量，然后使用上面的通解公式求出 $x(t)$。

例如，假设我们有以下线性齐次状态方程：

$\dot{x} = \begin{bmatrix} -2 & 1 \\ 1 & -2 \end{bmatrix} x$

我们可以使用 eig 函数求解特征值和特征向量：

[A, V] = eig([-2 1; 1 -2])

得到的结果为：

A =

-3 0
0 -1

V =

-0.7071 -0.7071
0.7071 -0.7071

因此，我们可以得到方程的通解：

$x(t) = c_1 e^{-3t} \begin{bmatrix}-0.7071 \\ 0.7071 \end{bmatrix} + c_2 e^{-t} \begin{bmatrix}-0.7071 \\ -0.7071 \end{bmatrix}$

2. 线性非齐次状态方程的解

假设我们有一个线性非齐次状态方程：

$\dot{x} = Ax + f(t)$

其中，$x$ 是一个 $n$ 维向量，$A$ 是一个 $n \times n$ 的矩阵，$f(t)$ 是一个已知的 $n$ 维向量函数。我们可以通过求解矩阵 $A$ 的特征值和特征向量，以及求解非齐次项的特解，得到方程的通解：

$x(t) = x_h(t) + x_p(t)$

其中，$x_h(t)$ 是对应于齐次方程的通解，$x_p(t)$ 是非齐次方程的一个特解。我们可以使用常数变易法或者伯努利法求解非齐次项的特解。在 MATLAB 中，我们可以使用 ode45 函数求解非齐次方程的数值解，也可以使用 dsolve 函数求解精确解。

例如，假设我们有以下线性非齐次状态方程：

$\dot{x} = \begin{bmatrix} -2 & 1 \\ 1 & -2 \end{bmatrix} x + \begin{bmatrix} 1 \\ 1 \end{bmatrix} \sin(t)$

我们可以使用 dsolve 函数求解精确解：

syms t
A = [-2 1; 1 -2];
f = [1; 1]*sin(t);
x = dsolve('Dx = A*x + f', 'x(0) = [0; 0]');

得到的结果为：

x(t) = [ cos(t)/2 + sin(t)/2, cos(t)/2 - sin(t)/2]*exp(-2*t)

3. 连续系统状态方程的离散化

假设我们有一个连续系统状态方程：

$\dot{x}(t) = Ax(t) + Bu(t)$

其中，$x(t)$ 是一个 $n$ 维向量，$u(t)$ 是一个 $m$ 维向量，$A$ 是一个 $n \times n$ 的矩阵，$B$ 是一个 $n \times m$ 的矩阵。我们可以使用欧拉法或者梯形法等数值方法将连续时间的状态方程离散化为离散时间的状态方程。

欧拉法的离散化公式为：

$x_{k+1} = x_k + h(Ax_k + Bu_k)$

其中，$x_k$ 和 $u_k$ 分别表示时间 $t = kh$ 时刻的状态和输入，$h$ 表示采样周期。在 MATLAB 中，我们可以使用 for 循环和矩阵运算来实现欧拉法的离散化。

例如，假设我们有以下连续系统状态方程：

$\dot{x}(t) = \begin{bmatrix} -2 & 1 \\ 1 & -2 \end{bmatrix} x(t) + \begin{bmatrix} 1 \\ 0 \end{bmatrix} u(t)$

我们可以使用欧拉法将其离散化：

A = [-2 1; 1 -2];
B = [1; 0];
h = 0.01; % 采样周期
x0 = [1; 0]; % 初始状态
u = @(t) sin(t); % 输入函数
N = 1000; % 采样点数
x = zeros(2, N+1); % 存储状态历史记录
x(:,1) = x0;
for k = 1:N
x(:,k+1) = x(:,k) + h*(A*x(:,k) + B*u(k*h));
end

得到的结果为 $x$ 是一个 $2 \times 1001$ 的矩阵，每一列表示一个时刻的状态。我们可以将其可视化：

plot(0:h:N*h, x(1,:), 'r', 0:h:N*h, x(2,:), 'b')
legend('x_1', 'x_2')
xlabel('Time (s)')
ylabel('State')