输入一个常系数齐次线性微分方程组并在不使用第三方库的前提下用python计算

假设有一个常系数齐次线性微分方程组：

$$\frac{d}{dt}\mathbf{x}(t) = A\mathbf{x}(t)$$

其中 $\mathbf{x}(t)$ 是一个 $n$ 维向量函数，$A$ 是一个 $n\times n$ 的常数矩阵。

我们可以将其写成矩阵形式：

$$\frac{d}{dt}\begin{pmatrix}x_1(t) \\ x_2(t) \\ \vdots \\ x_n(t)\end{pmatrix} = \begin{pmatrix}a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn}\end{pmatrix}\begin{pmatrix}x_1(t) \\ x_2(t) \\ \vdots \\ x_n(t)\end{pmatrix}$$

我们可以使用 NumPy 来计算这个微分方程组的解。首先，我们需要定义矩阵 $A$ 和初始向量 $\mathbf{x}(0)$：

import numpy as np

A = np.array([[a11, a12, ..., a1n], [a21, a22, ..., a2n], ..., [an1, an2, ..., ann]])
x0 = np.array([x1_0, x2_0, ..., xn_0])

接下来，我们需要定义一个函数来计算微分方程组的导数：

def dxdt(x, t, A):
    return np.dot(A, x)

注意，这个函数的输入是一个向量 $x$ 和一个时间 $t$，输出是向量 $x$ 的导数。我们使用 `np.dot` 函数来计算矩阵乘法。

最后，我们可以使用 SciPy 的 `solve_ivp` 函数来求解微分方程组：

from scipy.integrate import solve_ivp

t_span = (t0, tf)
sol = solve_ivp(dxdt, t_span, x0, args=(A,))

其中 `t0` 和 `tf` 是时间范围，`args` 参数传递了矩阵 $A$。`sol` 是一个 `scipy.integrate.OdeResult` 对象，可以使用 `sol.t` 和 `sol.y` 属性来获得时间和解向量。