计算按列的主要顺序存储的函数

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 = ....reshape(-1，1)。并读取我链接到的文档，以及我使用的所有函数。

现在有了y方向的列和x中的行，您可以将Function调用为

val = Function(delta_x, delta_y)

将二维矩阵val排列成一维阵列的操作称为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')

你可以把它做成一个单线，但这是在推动它。