Evaluation a function stored in column major order

Question

I am trying to evaluate a function at discretized point and stored in column-major order like this:

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))

How could I approach to evaluate the discretization for this function?

user_2988730 · Answer

Good news: your function properly uses numpy operations, so is completely vectorized. That means that you can evaluate it at every element of the input arrays.

The shape of the inputs don't have to match exactly. They just have to broadcast together. That means that only non-singleton dimensions need to match.

So start by creating the appropriate input arrays. Numpy provides the tools to do this elegantly without looping:

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

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

I deliberately used two different ways to create the coordinate arrays, to serve as an example. In practice, you'd want to use one of the two methods.

The index [:, None] turns delta_y into a column vector. None introduces a new singleton axis. There are any number of Other ways to do the same thing, like `delta_y = ....reshape(-1, 1). And read the docs I linked to, and for all the functions I used.

Now that you have a column in the y direction and a row in the x, you can call Function as just

val = Function(delta_x, delta_y)

The operation of arranging the 2D matrix val into a 1D array is called raveling. By default, it uses the default row-major order that numpy uses in memory. This order is also called "C" order. The alternative arrangement is to interpret the array in column major order, like Matlab does. This is called Fortran order. It will require a copy of the data since that's not how the elements are laid out in memory.

One way to ravel in Fortran order:

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

An alternative is to transpose and use C order:

feval = val.T.ravel()

The last two lines can be combined, so you end up with 3 lines:

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')

You could make it into a one-liner, but that's pushing it.

Evaluation a function stored in column major order

1 Answers