Matlab翻译：新问题（复ginzburg-landau方程）

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)

1条回答

网友

1楼 · 发布于 2024-04-26 04:45:12

错误在：

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

q[:, t+1]索引2d数组（可能不是np.matrix，后者更像MATLAB）。此索引将维度数减少1，因此错误消息中的（220，）形状。你知道吗

错误信息显示RHS为（220220）。这个形状可能来自inverted和forgold。np.array(q[:,t])是1d。将一个（220220）乘以一个（220，）是可以的，但是你不能将这个正方形数组放入1d槽中。你知道吗

在错误行中使用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`like运算的函数。你知道吗

注意，在numpy*中是元素乘法，即MATLAB的.*。np.dot和@用于矩阵乘法。你知道吗

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

会有用的

相关问题更多 >

编程相关推荐

热门问题

热门文章