Python实现线性微分方程组f(k)函数报错求助

需要实现函数f(k)，满足以下要求：

• 函数定义：运动开始后前10秒内，系统状态向量x(t)的欧几里得范数平方的时间积分的平方根
• 系统模型：线性微分方程组dx/dt = Ax + bu，控制信号u = kx（x为2维向量）
• 输入参数：矩阵A、向量b、初始状态x0
• 特殊情况：若系统不稳定，f(k)返回-1

用户尝试的代码如下：

import numpy as np
from scipy.integrate import solve_ivp
from scipy.integrate import quad
from numpy.linalg import norm


def f(A, b, k, x0):
    M = A + b * k
    eigenvalues = np.linalg.eig(M)[0]
    if eigenvalues[0] > 0 and eigenvalues[1] > 0:
        return -1
    else:
        def odefun(t, x):
            x = x.reshape([2, 1])     
            dx = M @ x      
            return dx.reshape(-1)

        # solve system
        x0 = x0.reshape(-1)
        sol = solve_ivp(odefun, [0, 10], x0)['y']
        print(sol.shape) # (2, 41)
        return quad(norm(sol), 0, 10)


def contra(A, b, x0):
    vals = []
    for k in range(10):
        vals.append(f(A, b, k, x0))
    return vals

A = np.array(([1, 2], [3, 4]))
b = np.array([5, 6])
x0 = np.array([7, 8])
print(contra(A, b, x0))

运行时触发错误：

ValueError: invalid callable given

错误原因分析 #

1. 稳定性判断逻辑错误：原代码仅当两个特征值都大于0时才判定系统不稳定，但线性系统不稳定的条件是存在至少一个特征值的实部大于0，而非所有特征值都大于0。
2. 积分参数错误：quad函数要求传入一个可调用的函数（输入为时间t，输出为对应时刻的被积函数值），但原代码直接传入了solve_ivp返回的数值数组sol，导致无法被quad识别。
3. 矩阵维度匹配错误：原代码中b * k是元素级乘法，无法正确构造闭环系统矩阵M，需将b转为列向量后再进行乘法操作。

修正后的代码 #

import numpy as np
from scipy.integrate import solve_ivp
from scipy.integrate import quad
from numpy.linalg import norm
from scipy.interpolate import interp1d


def f(A, b, k, x0):
    # 构造闭环系统矩阵：将b转为列向量后与k相乘，保证矩阵加法维度匹配
    M = A + b.reshape(-1, 1) * k
    # 计算特征值并检查稳定性：所有特征值实部≤0才稳定（加小阈值避免浮点误差）
    eigenvalues = np.linalg.eig(M)[0]
    if any(np.real(eigenvalues) > 1e-8):
        return -1
    
    # 定义微分方程
    def odefun(t, x):
        x = x.reshape(-1, 1)
        dx = M @ x
        return dx.flatten()
    
    # 求解微分方程，启用dense_output以支持任意时刻状态查询
    sol = solve_ivp(odefun, [0, 10], x0.flatten(), dense_output=True)
    # 构造状态插值函数，支持通过任意时刻t获取对应状态x(t)
    x_interp = interp1d(sol.t, sol.y, axis=1, kind='cubic')
    
    # 定义被积函数：t时刻状态的欧几里得范数平方
    def integrand(t):
        x_t = x_interp(t)
        return norm(x_t) ** 2
    
    # 计算积分并取平方根得到最终结果
    integral, _ = quad(integrand, 0, 10)
    return np.sqrt(integral)


def contra(A, b, x0):
    vals = []
    for k in range(10):
        vals.append(f(A, b, k, x0))
    return vals


# 测试示例
A = np.array([[1, 2], [3, 4]])
b = np.array([5, 6])
x0 = np.array([7, 8])
print(contra(A, b, x0))

关键修正点说明 #

• 闭环矩阵构造：将b转为列向量后与标量k相乘，确保矩阵加法的维度合法性。
• 稳定性判断：检查所有特征值的实部，只要存在实部大于阈值（1e-8）的特征值，立即判定系统不稳定并返回-1。
• 状态插值：使用interp1d构造状态的插值函数，实现任意时刻t的状态查询，满足quad对可调用被积函数的要求。
• 积分逻辑：定义integrand函数计算每个时刻的范数平方，传入quad完成积分后，取平方根得到符合定义的最终结果。

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