如何广播我的矩阵求幂

计算科学 Python scipy 麻木的
2021-12-15 03:44:48

我正在尝试-无济于事-在以下代码中广播以下代码Python：

import numpy as np
from scipy.linalg import expm

def ExactSol(c, omega0, x0, v0, T, p):
    # Docstring
    """
    Returns the solution of
    x'' + c x' + w0^2 x = 0
    embodied with initial condition x(0) = x0 and x'(0) = v0
    i.e. the function X(t) = exp(tA) * X0
    with X = (x, v), X0 = (x0, v0), A = (0       1)
                                        (-w0^2  -c)
    over [0,T] cut in p sub-intervals
    Positional arguments:
    c      -- damping parameter
    omega0 -- pulsation
    x0     -- initial position
    v0     -- initial speed
    T      -- final time
    p      -- the number of samples in [0,T]
    Keyword arguments:
    """
    # Body
    time = np.linspace(0, T, p + 1)
    time = time[:, np.newaxis, np.newaxis]
    X0   = np.array([x0, v0])
    A    = np.array([[0, 1], [- omega0**2, - c]])
    B    = time * A
    E    = np.zeros(B.shape)
    for i in range(0, B.shape[0]):
        E[i,:,:] = expm(B[i,:,:])
    return time[:,0,0], E@X0

我想避免使用for循环。

1个回答

自从A只是一个 2x2 矩阵，您可以对矩阵求幂进行硬编码。请注意，如果我们可以写A=UDU−1， 在哪里D是特征值的对角矩阵A和列U是相关的特征向量，那么eA=UeDU−1和etA=UetDU−1， 在哪里t是一个标量。

剩下的就是计算特征值和特征向量。它有点难看，但它可能比你的代码更快。基于此。请注意，我没有解决有缺陷的矩阵的情况。

import numpy as np
from numpy.linalg import eig,inv

def ExactSol(c, omega0, x0, v0, T, p):
    time = np.linspace(0, T, p + 1)
    A = np.array([[0, 1], [- omega0**2, - c]])
    e,U = eig(A)
    tol = 1.e-6
    E = np.zeros( (p+1,2,2),dtype=float)
    if abs(e[0]-e[1]) < tol:
        # this is a defective matrix, which must be handled differently
        raise UserWarning("Defective matrix not implemented")
    else:
        # the eigenvalues are distinct -- possibly complex, but 
        # E will always be real
        Uinv=inv(U)
        for k in range(2):
            E += np.real(np.exp(time*e[k])[:,np.newaxis,np.newaxis] *
                         np.dot(U[:,k,np.newaxis], Uinv[np.newaxis,k,:])[np.newaxis,:,:] )

    X0 = np.array([x0, v0])
    X = np.dot(E,X0[:,np.newaxis])
    return time, X

编辑：一些快速的时间信息。我设置omega0=1.0, c=0.1, x0=0.0, v0=1.0, T=100.0, 和p=10000. 使用你的代码，这个函数需要 2.95 秒，而使用我的代码，在我的机器上需要 0.0025 秒。因此，我将效率提高了 1000 多倍。

其它你可能感兴趣的问题
上一篇另一种评估均匀立方体引力的方法？ 下一篇“具有对角线优势”的稠密矩阵的预处理器