计算以列主顺序存储的函数

Question

我试图在离散点计算一个函数，并以列的大顺序存储，如下所示：

import numpy as np; 

N = 3 ##
n = N * N
h = 1 / (N + 1) # step size
h2 = h**2  # 


deltaX = np.zeros(N)
deltaY = np.zeros(N);

def Function(x, y):
    output = -20. * np.pi * np.sin(2 * np.pi * x) * sin(4 * np.pi * y)
    return output
    
## Equally spaced delta: 
for i in range(1, N + 1): 
    deltaX[i - 1] = i * h; 
    deltaY[i - 1] = i * h; 

### Lexicographic Row order  ###

### Evaluation of function at deltaX and deltaY 

feval = np.zeros((n, 1))

如何评估此函数的离散化

Answer 1

好消息：您的函数正确地使用了numpy操作，因此完全矢量化了。这意味着您可以在输入数组的每个元素上计算它

输入的形状不必完全匹配。他们只需要在一起。这意味着只有非单例维度需要匹配

因此，从创建适当的输入数组开始。Numpy提供了一些工具，可以在不循环的情况下优雅地执行此操作：

N = 3
h = 1 / (N + 1)

delta_x = np.arange(1., N + 1.) * h
delta_y = np.linspace(h, N * h, N)[:, None]

作为示例，我特意使用了两种不同的方法来创建坐标数组。实际上，您可能希望使用这两种方法中的一种

索引[:, None]将delta_y转换为列向量None引入了一个新的单态轴。有很多其他方法可以做同样的事情，比如'delta_y=..重塑（-1，1）。 阅读我链接到的文档，以及我使用的所有功能

现在您在y方向有一列，在x方向有一行，您可以调用Function作为

val = Function(delta_x, delta_y)

将2D矩阵val排列成1D数组的操作称为raveling。默认情况下，它使用numpy在内存中使用的默认行主顺序。这种顺序也称为“C”顺序。另一种安排是按照列的主要顺序解释数组，就像Matlab那样。这就是所谓的Fortran顺序。它将需要一份数据的副本，因为这不是元素在内存中的布局方式

按Fortran顺序进行ravel的一种方法：

feval = val.ravel(order='F')

另一种方法是转置并使用C顺序：

feval = val.T.ravel()

最后两行可以合并，因此您将得到3行：

delta_x = h * np.arange(1., N + 1.)
delta_y = h * np.arange(1., N + 1.)[:, None]
feval = Function(delta_x, delta_y).ravel(order='F')

你可以把它做成一条直线，但那是在推动它