# Demystifying `meshgrid()`

We use the Numpy [meshgrid()](https://docs.scipy.org/doc/numpy/reference/generated/numpy.meshgrid.html) function to plot 2D arrays but it can be a bit mysterious as to what the "meshgrid" actually does or look like.

The docs say: 

_Return coordinate matrices from coordinate vectors._

_Make N-D coordinate arrays for vectorized evaluations of N-D scalar/vector fields over N-D grids, given one-dimensional coordinate arrays x1, x2,..., xn._


Things to keep in mind:
- NumPy functions evaluate input _element-wise_ : `f(v)` with 1D array `v` gives `[f(v[0]), f(v[1]), ...]` (a 1D array) or for a 2D array `a`, `f(a)` produces a 2D array `[[f(a[0,0]), f(a[0, 1]),...], [f(a[1,0]), f(a[1, 1]),...]]`
- matplotlib's 2D plotting functions take as input 2D arrays (together with optional arrays to specify coordinate values (as opposed to plotting "by index"))

## Visualization with a Gaussian
Testing with the initial condition, an isotropic 2D [Gaussian](http://mathworld.wolfram.com/GaussianFunction.html)

$$
g_2(x, y) = \frac{1}{2\pi\sigma^2} \exp\left[-\frac{(x-x_0)^2 + (y-y_0)^2}{2\sigma^2}\right]
$$

Note that the 2D Gaussian can be written as a product of two 1D Gaussians

\begin{align}
g_2(x, y) &= \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left[-\frac{(x-x_0)^2}{2\sigma^2}\right] \cdot 
          \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left[-\frac{(y-y_0)^2}{2\sigma^2}\right]\\
        &= g_1(x) \cdot g_1(y)
\end{align}

Plot $g_2(x, y)$ by (1) constructing it from the product of $g_1$s and (2) from the explicit formula. Use [`np.meshgrid`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.meshgrid.html) to evaluate the function on a grid.

Use an *asymmetric* grid ($L_x=0.5, N_x=50$ and $L_y=1, N_y=100$) to clearly see the shape of the underlying arrays.


```python
import numpy as np
import matplotlib.pyplot as plt
```


```python
%matplotlib inline
```


```python
L = 0.5
Lx, Ly = 0.5, 1
Nx, Ny = 50, 100

x = np.linspace(0, Lx, Nx+1)  # need N+1!
y = np.linspace(0, Ly, Ny+1)
X, Y = np.meshgrid(x, y)
    
def gaussian2D(x, y, u0=0.05, x0=None, y0=None, sigma=0.1*L):
    x0 = np.mean(x) if x0 is None else x0
    y0 = np.mean(y) if y0 is None else y0
    return u0/(2*np.pi*sigma) * np.exp(-((x-x0)**2 + (y-y0)**2)/(2*sigma**2))

def gaussian(x, u0=0.05, x0=None, sigma=0.1*L):
    x0 = np.mean(x) if x0 is None else x0
    return u0/np.sqrt(2*np.pi*sigma) * np.exp(-(x-x0)**2 / (2*sigma**2))
```

#### Meshgrid 


```python
print(x.shape, y.shape)
print(X.shape, Y.shape)
print(X)
print(Y)
```

    (51,) (101,)
    (101, 51) (101, 51)
    [[0.   0.01 0.02 ... 0.48 0.49 0.5 ]
     [0.   0.01 0.02 ... 0.48 0.49 0.5 ]
     [0.   0.01 0.02 ... 0.48 0.49 0.5 ]
     ...
     [0.   0.01 0.02 ... 0.48 0.49 0.5 ]
     [0.   0.01 0.02 ... 0.48 0.49 0.5 ]
     [0.   0.01 0.02 ... 0.48 0.49 0.5 ]]
    [[0.   0.   0.   ... 0.   0.   0.  ]
     [0.01 0.01 0.01 ... 0.01 0.01 0.01]
     [0.02 0.02 0.02 ... 0.02 0.02 0.02]
     ...
     [0.98 0.98 0.98 ... 0.98 0.98 0.98]
     [0.99 0.99 0.99 ... 0.99 0.99 0.99]
     [1.   1.   1.   ... 1.   1.   1.  ]]


Note that the two meshgrids have the same shape: `Nx x Ny`.

They also contain either constant rows or constant columns. We will see how this is useful — but you can already imagine what happens if you "overlay" `X` and `Y` and read of the entries in each "overlayed" cell: the coordinates of the cell:


```python
X[0, 0], Y[0, 0]
```




    (0.0, 0.0)




```python
X[1, 0], Y[1, 0]
```




    (0.0, 0.01)



#### Plot 1D Gaussians
Plot the 1D Gaussians using the meshgrid components.


```python
gaussian(X).shape
```




    (101, 51)




```python
fig = plt.figure()
ax = fig.add_subplot(111)
ax.contourf(X, Y, gaussian(X, sigma=0.1), 30)
ax.set_aspect(1)
plt.show()
```


```python
gaussian(Y).shape
```




    (101, 51)




```python
ax = plt.subplot(131)
ax.contourf(X, Y, gaussian(X, sigma=0.1), 30)
ax.set_aspect(1)
ax = plt.subplot(132)
ax.contourf(X, Y, gaussian(Y, sigma=0.1), 30)
ax.set_aspect(1)
ax = plt.subplot(133)
ax.contourf(X, Y, gaussian(X, sigma=0.1) * gaussian(Y, sigma=0.1), 30)
ax.set_aspect(1)
```

#### 2D gaussian 


```python
gaussian2D(X, Y).shape
```




    (101, 51)




```python
ax = plt.subplot(111)
ax.contourf(X, Y, gaussian2D(X, Y, sigma=0.1), 30)
ax.set_aspect(1)
```


```python

```

80c43c4e941d3a7cc516eb5d494bd3b4289edfe3
ipynb
Jupyter Notebook
15_PDEs/meshgrid.ipynb
ASU-CompMethodsPhysics-PHY494/PHY494-resources-2020
20e08c20995eab567063b1845487e84c0e690e96
CC-BY-4.0
15_PDEs/meshgrid.ipynb
ASU-CompMethodsPhysics-PHY494/PHY494-resources-2020
20e08c20995eab567063b1845487e84c0e690e96
CC-BY-4.0
15_PDEs/meshgrid.ipynb
ASU-CompMethodsPhysics-PHY494/PHY494-resources-2020
20e08c20995eab567063b1845487e84c0e690e96
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES