从 Matlab 翻译：新问题（复杂的 ginzburg 朗道方程

Question

感谢您对我上一篇文章的所有建设性批评。我做了一些更改，但是我的代码仍然无法正常工作，我不知道为什么。当我运行这个版本时，我收到一个关于 matmul 中遇到的无效错误的运行时警告。 我的代码给出为

from __future__ import division
import numpy as np
from scipy.linalg import eig
from scipy.linalg import toeplitz


def poldif(*arg):
    """
    Calculate differentiation matrices on arbitrary nodes.
    Returns the differentiation matrices D1, D2, .. DM corresponding to the
    M-th derivative of the function f at arbitrarily specified nodes. The
    differentiation matrices can be computed with unit weights or
    with specified weights.
    Parameters
    ----------
    x       : ndarray
              vector of N distinct nodes
    M       : int
              maximum order of the derivative, 0 < M <= N - 1
    OR (when computing with specified weights)
    x       : ndarray
              vector of N distinct nodes
    alpha   : ndarray
              vector of weight values alpha(x), evaluated at x = x_j.
    B       : int
              matrix of size M x N, where M is the highest derivative required.
              It should contain the quantities B[l,j] = beta_{l,j} =
              l-th derivative of log(alpha(x)), evaluated at x = x_j.
    Returns
    -------
    DM : ndarray
         M x N x N  array of differentiation matrices
    Notes
    -----
    This function returns  M differentiation matrices corresponding to the
    1st, 2nd, ... M-th derivates on arbitrary nodes specified in the array
    x. The nodes must be distinct but are, otherwise, arbitrary. The
    matrices are constructed by differentiating N-th order Lagrange
    interpolating polynomial that passes through the speficied points.
    The M-th derivative of the grid function f is obtained by the matrix-
    vector multiplication
    .. math::
    f^{(m)}_i = D^{(m)}_{ij}f_j
    This function is based on code by Rex Fuzzle
    https://github.com/RexFuzzle/Python-Library
    References
    ----------
    ..[1] B. Fornberg, Generation of Finite Difference Formulas on Arbitrarily
    Spaced Grids, Mathematics of Computation 51, no. 184 (1988): 699-706.
    ..[2] J. A. C. Weidemann and S. C. Reddy, A MATLAB Differentiation Matrix
    Suite, ACM Transactions on Mathematical Software, 26, (2000) : 465-519
    """
    if len(arg) > 3:
        raise Exception('number of arguments is either two OR three')

    if len(arg) == 2:
        # unit weight function : arguments are nodes and derivative order
        x, M = arg[0], arg[1]
        N = np.size(x)
        # assert M<N, "Derivative order cannot be larger or equal to number of points"
        if M >= N:
            raise Exception("Derivative order cannot be larger or equal to number of points")
        alpha = np.ones(N)
        B = np.zeros((M, N))

    elif len(arg) == 3:
        # specified weight function : arguments are nodes, weights and B  matrix
        x, alpha, B = arg[0], arg[1], arg[2]
        N = np.size(x)
        M = B.shape[0]

    I = np.eye(N)  # identity matrix
    L = np.logical_or(I, np.zeros(N))  # logical identity matrix
    XX = np.transpose(np.array([x, ] * N))
    DX = XX - np.transpose(XX)  # DX contains entries x(k)-x(j)
    DX[L] = np.ones(N)  # put 1's one the main diagonal
    c = alpha * np.prod(DX, 1)  # quantities c(j)
    C = np.transpose(np.array([c, ] * N))
    C = C / np.transpose(C)  # matrix with entries c(k)/c(j).
    Z = 1 / DX  # Z contains entries 1/(x(k)-x(j)
    Z[L] = 0  # eye(N)*ZZ;                # with zeros on the diagonal.
    X = np.transpose(np.copy(Z))  # X is same as Z', but with ...
    Xnew = X

    for i in range(0, N):
        Xnew[i:N - 1, i] = X[i + 1:N, i]

    X = Xnew[0:N - 1, :]  # ... diagonal entries removed
    Y = np.ones([N - 1, N])  # initialize Y and D matrices.
    D = np.eye(N)  # Y is matrix of cumulative sums

    DM = np.empty((M, N, N))  # differentiation matrices

    for ell in range(1, M + 1):
        Y = np.cumsum(np.vstack((B[ell - 1, :], ell * (Y[0:N - 1, :]) * X)), 0)  # diags
        D = ell * Z * (C * np.transpose(np.tile(np.diag(D), (N, 1))) - D)  # off-diags
        D[L] = Y[N - 1, :]
        DM[ell - 1, :, :] = D

    return DM

def herdif(N, M, b=1):
    """
    Calculate differentiation matrices using Hermite collocation.
    Returns the differentiation matrices D1, D2, .. DM corresponding to the
    M-th derivative of the function f, at the N Chebyshev nodes in the
    interval [-1,1].
    Parameters
    ----------
    N   : int
          number of grid points
    M   : int
          maximum order of the derivative, 0 < M < N
    b   : float, optional
          scale parameter, real and positive
    Returns
    -------
    x  : ndarray
         N x 1 array of Hermite nodes which are zeros of the N-th degree
         Hermite polynomial, scaled by b
    DM : ndarray
         M x N x N  array of differentiation matrices
    Notes
    -----
    This function returns  M differentiation matrices corresponding to the
    1st, 2nd, ... M-th derivates on a Hermite grid of N points. The
    matrices are constructed by differentiating N-th order Hermite
    interpolants.
    The M-th derivative of the grid function f is obtained by the matrix-
    vector multiplication
    .. math::
    f^{(m)}_i = D^{(m)}_{ij}f_j
    References
    ----------
    ..[1] B. Fornberg, Generation of Finite Difference Formulas on Arbitrarily
    Spaced Grids, Mathematics of Computation 51, no. 184 (1988): 699-706.
    ..[2] J. A. C. Weidemann and S. C. Reddy, A MATLAB Differentiation Matrix
    Suite, ACM Transactions on Mathematical Software, 26, (2000) : 465-519
    ..[3] R. Baltensperger and M. R. Trummer, Spectral Differencing With A
    Twist, SIAM Journal on Scientific Computing 24, (2002) : 1465-1487
    """
    if M >= N - 1:
        raise Exception('number of nodes must be greater than M - 1')

    if M <= 0:
        raise Exception('derivative order must be at least 1')

    x = herroots(N)  # compute Hermite nodes
    alpha = np.exp(-x * x / 2)  # compute Hermite  weights.

    beta = np.zeros([M + 1, N])

    # construct beta(l,j) = d^l/dx^l (alpha(x)/alpha'(x))|x=x_j recursively
    beta[0, :] = np.ones(N)
    beta[1, :] = -x

    for ell in range(2, M + 1):
        beta[ell, :] = -x * beta[ell - 1, :] - (ell - 1) * beta[ell - 2, :]

    # remove initialising row from beta
    beta = np.delete(beta, 0, 0)

    # compute differentiation matrix (b=1)
    DM = poldif(x, alpha, beta)
    # scale nodes by the factor b
    x = x / b

    # scale the matrix by the factor b
    for ell in range(M):
        DM[ell, :, :] = (b ** (ell + 1)) * DM[ell, :, :]

    return x, DM

def herroots(N):
    """
    Compute roots of the Hermite polynomial of degree N
    Parameters
     ----------
    N   : int
          degree of the Hermite polynomial
    Returns
    -------
    x  : ndarray
         N x 1 array of Hermite roots
    """

    # Jacobi matrix
    d = np.sqrt(np.arange(1, N))
    J = np.diag(d, 1) + np.diag(d, -1)

    # compute eigenvalues
    mu = eig(J)[0]

    # return sorted, normalised eigenvalues
    # real part only since all roots must be real.
    return np.real(np.sort(mu) / np.sqrt(2))


a = 1-1j
b = 2+0.2j
c1 = 0.34
c2 = 0.005

alpha1 = (4*c2/a)**0.25
alpha2 = b/2*a

Nx = 220;

# hermite differentiation matrices
[x,D] = herdif(Nx, 2, np.real(alpha1))
D1 = D[0,:]
D2 = D[1,:]

# integration weights
diff = np.diff(x)
#print(len(diff))
p = np.concatenate([np.zeros(1), diff])
q = np.concatenate([diff, np.zeros(1)])
w = (p + q)/2
Q = np.diag(w)

#Discretised operator
const = c1*np.diag(np.ones(len(x)))-c2*(np.diag(x)*np.diag(x))
#print(const)
A = a*D2 - b*D1 + const

##### Timestepping

tmax = 200
tmin = 0
dt = 1
n = (tmax - tmin)/dt
tvec = np.linspace(0,tmax,n, endpoint = True)

#(len(tvec))

q = np.zeros((Nx, len(tvec)),dtype=complex)
f = np.zeros((Nx, len(tvec)),dtype=complex)

q0 = np.ones(Nx)*10**4
q[:,0] = q0
#print(q[:,0])
#print(q0)

# qnew - qold = dt*Aqold + dt*N(qold,qold,qold)
# qnew - qold = dt*Aqnew - dt*N(qold,qold,qold)
# therefore qnew - qold = 0.5*dtAqold + 0.5*dt*Aqnew + dtN(qold,qold,qold)
# rearranging to give qnew( 1- 0.5Adt) = (1 + 0.5Adt) + dt N(qold,qold,qold)

from numpy.linalg import inv

inverted = inv(np.eye(Nx)-0.5*A*dt)
forqold = (np.eye(Nx) + 0.5*A*dt)
firstterm = np.matmul(inverted,forqold)

for t in range(0, len(tvec)-1):
    nl = abs(np.square(q[:,t]))*q[:,t]
    q[:,t+1] = np.matmul(firstterm,q[:,t]) - dt*np.matmul(inverted,nl)

hermite 分化矩阵可以在网上找到并且在不同的文件中。此代码在五次交互后爆炸，我无法理解，因为我看不出它在此处找到的 matlab 中有何不同https://www.bagherigroup.com/research/open-source-codes/

我非常感谢任何帮助。

Answer 1

错误：

q[:,t+1] = inverted*forgold*np.array(q[:,t]) + inverted*dt*np.array(nl)

q[:, t+1] 索引一个二维数组（可能不是np.matrix，它更像 MATLAB）。此索引将维度数减少 1，因此错误消息中的 (220,) 形状。

错误消息显示 RHS 为 (220,220)。该形状可能来自inverted 和forgold。 np.array(q[:,t]) 是 1d。将 (220,220) 乘以 (220,) 是可以的，但您不能将该方形数组放入一维槽中。

错误行中np.array 的两种用法都是多余的。他们的论点已经是ndarray。

至于循环，可能是必要的。看起来q[:,t+1] 是 q[:,t], a serial, rather than parallel operation. Those are harder to render as 'vectorized' (unless you can usecumsum` 类操作的函数。

请注意，numpy * 是元素乘法，即 MATLAB 的 .*。 np.dot 和 @ 用于矩阵乘法。

q[:,t+1]= invert@q[:,t]

会工作