求解矩阵微分方程Ax'=Bx+b的Python实现求助

Got it, let's tackle this problem step by step. First, we need to rearrange your equation into a form that Python's numerical solvers can work with easily—since most ODE solvers expect the standard first-order linear form: (\mathbf{x}' = M\mathbf{x} + \mathbf{c}).

Step 1: Rewrite the Equation #

Assuming matrix (A) is invertible (we'll cover non-invertible cases later), multiply both sides of your original equation by (A^{-1}):

• (M = A^{-1}B) (the coefficient matrix for (\mathbf{x}))
• (\mathbf{c} = A^{-1}\mathbf{b}) (the constant vector term)

This transforms your equation into the standard form that tools like scipy.integrate.solve_ivp can handle.

Step 2: Set Up the Code #

First, make sure you have the necessary libraries installed:

pip install numpy scipy matplotlib

Here's a complete, runnable example with sample matrices/vectors (replace these with your actual values):

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Define your system parameters (replace these with your real matrices/vector)
N = 2  # Dimension of your system
A = np.array([[2, 1], [1, 2]])  # Example invertible matrix
B = np.array([[0, 1], [-1, 0]])
b = np.array([1, 0])

# Check if A is invertible first—critical step!
if np.isclose(np.linalg.det(A), 0):
    raise ValueError("Matrix A is singular (non-invertible). See note below for how to handle this.")
A_inv = np.linalg.inv(A)

# Compute M and c for the standard ODE form
M = A_inv @ B
c = A_inv @ b

# Define the differential equation function
def dxdt(t, x):
    return M @ x + c

# Set initial conditions and time range to solve over
x0 = np.array([0, 0])  # Your initial state vector
t_span = [0, 10]       # Start and end time
t_eval = np.linspace(t_span[0], t_span[1], 100)  # Points to evaluate the solution at

# Solve the ODE
solution = solve_ivp(dxdt, t_span, x0, t_eval=t_eval)

# Visualize the results
plt.figure(figsize=(10, 6))
for i in range(N):
    plt.plot(solution.t, solution.y[i], label=f'x_{i+1}(t)')
plt.xlabel('Time t')
plt.ylabel('State Variables')
plt.title('Solution to $A\\mathbf{x}\' = B\\mathbf{x} + \\mathbf{b}$')
plt.legend()
plt.grid(True)
plt.show()

Edge Cases & Additional Notes #

• If A is singular: If (A) doesn't have an inverse, you can use the Moore-Penrose pseudoinverse (np.linalg.pinv(A)) instead of the regular inverse. However, this doesn't guarantee a unique solution—you'll need to check if the original system is consistent first.
• Analytical Solution (optional): For linear systems, you can also compute the analytical solution using matrix exponentials. Here's a quick snippet for that:

  from scipy.linalg import expm
  
  def analytical_sol(t, x0, M, c):
      exp_Mt = expm(M * t)
      if np.isclose(np.linalg.det(M), 0):
          # Handle singular M (steady-state solution if applicable)
          return exp_Mt @ x0 + t * c
      else:
          M_inv = np.linalg.inv(M)
          return exp_Mt @ x0 + M_inv @ (exp_Mt - np.eye(N)) @ c
  
  # Compute and plot alongside numerical solution
  x_analytical = np.array([analytical_sol(t, x0, M, c) for t in t_eval]).T
  plt.figure(figsize=(10,6))
  for i in range(N):
      plt.plot(solution.t, solution.y[i], label=f'Numerical x_{i+1}')
      plt.plot(t_eval, x_analytical[i], '--', label=f'Analytical x_{i+1}')
  plt.legend()
  plt.grid(True)
  plt.show()
  

内容的提问来源于stack exchange，提问作者Ricku