<jupyter_start><jupyter_text>Boundary value problems--- Poisson boundary value problemConsider the Poisson equation$$-\nabla^2\varphi = -\frac{\partial^2 \varphi}{\partial x^2} - \frac{\partial^2 \varphi}{\partial y^2} = f(\varphi,x,y)$$With Dirichlet boundary conditions\begin{align*}\varphi(0, y) &= 0, \quad 0\le y\le H \\\varphi(L, y) &= 0, \quad 0\le y\le H \\\varphi(x, 0) &= 0, \quad 0\le x\le L \\\varphi(x, H) &= 0, \quad 0\le x\le L\end{align*}**Applications:**- Model for laminar flow of Newtonian fluid through a rectangular duct- Steady state diffusion or heat conduction  Solution domain- 2D grid with boundary points in green and interior points in blue- The grid has $N_x$ and $N_y$ points in the $x$ and $y$ direction, respectively- The 2D grid is internally stored in one long array- The point $\phi(x_i, y_j) =\phi_{i,j}$ has the position $n=iN_y + j$  Method of finite differencesThe value on the grid boundaries are given by the Dirichlet boundary conditions $\phi_{i,j} = 0$ for\begin{align*}i=0 & \quad & 0\le j\le N_y-1 \\i=N_x-1v& \quad & 0\le j\le N_y-1 \\0\le i\le N_x-1 & \quad & j=0 \\0\le i\le N_x-1 & \quad & j=N_y-1\end{align*}In the interior we are using the approximation of the 2nd partial derivatives\begin{align*}\left.\frac{d^2\phi}{dx^2}\right|_{x_i, y_j} &= \frac{\phi(x_i+\Delta x, y_j) - 2 \phi(x_i,y_j) + \phi(x_i-\Delta x,y_j)}{(\Delta x)^2} = \frac{\phi(x_{i+1}, y_j) - 2 \phi(x_i,y_j) + \phi(x_{i-1},y_j)}{(\Delta x)^2} \\\left.\frac{d^2\phi}{dy^2}\right|_{x_i, y_j} &= \frac{\phi(x_i, y_j+\Delta y) - 2 \phi(x_i,y_j) + \phi(x_i,y_j-\Delta y)}{(\Delta y)^2} = \frac{\phi(x_i, y_{j+1}) - 2 \phi(x_i,y_j) + \phi(x_i,y_{j-1})}{(\Delta y)^2} \end{align*}So the discretisation is given by\begin{align*}-\frac{\phi(x_{i+1}, y_j) + 2 \phi(x_i,y_j) - \phi(x_{i-1},y_j)}{(\Delta x)^2} - \frac{\phi(x_i, y_{j+1}) + 2 \phi(x_i,y_j) - \phi(x_i,y_{j-1})}{(\Delta y)^2} &= f(x_i, y_j)\end{align*}for $1\le i \le N_x-2$ and $1\le j \le N_y-2$With the mapping $n=iN_y + j$ we get\begin{align*}A_{n,n-N_y} \phi_{n-N_y} &+ A_{n,n-1} \phi_{n-1} + A_{n,n} \phi_{n} \\ & + A_{n,n+1} \phi_{n+1} + A_{n,n+N_y} \phi_{n+N_y} = b_n\end{align*}for $1\le i \le N_x-2$ and $1\le j \le N_y-2$The parameters $A$ are given by\begin{align*}A_{n,n-N_y} = A_{n,n+N_y} &= \left[\frac{-1}{(\Delta x)^2}\right] \\A_{n,n-1} = A_{n,n+1} &= \left[\frac{-1}{(\Delta y)^2}\right] \\A_{n,n} &= \left[\frac{2}{(\Delta x)^2} + \frac{2}{(\Delta y)^2}\right] \\b_n &= f(x_i, y_j)\end{align*}  Numerical solutionFirst, we define the grid spacing in the $x$ and $y$ directions. Then we use the [`numpy.meshgrid`](https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html) function to generate 2D arrays for the $x$ and $y$ coordinates of the entire 2D grid. By default, `meshgrid` uses Cartesian indexing so that inputs of length $N_x$ and $N_y$ give outputs of shape ($N_y$, $N_x$). Finally, we plot the grid points.<jupyter_code>import numpy as np

# x and y length and number of grid points
x_length = 2
Nx = 45
y_length = 1
Ny = 30

# Define the meshgrid
x = np.linspace(0, x_length, Nx)
dx = x[1] - x[0]
y = np.linspace(0, y_length, Ny)
dy = y[1] - y[0]

X, Y = np.meshgrid(x, y)
print("The X and Y arrays are of shape {}".format(np.shape(X)))

# Plot the grid points
import matplotlib.pyplot as plt
plt.scatter(X, Y)
plt.xlabel("x axis")
plt.ylabel("y axis")
plt.title("Plot of the grid points")
plt.show()<jupyter_output><empty_output><jupyter_text>Next we allocate storage for the large array $A$ and the right hand side $b$. We will also define the boundary conditions.<jupyter_code># Total size
Ntot = Nx * Ny

# Allocate space for the arrays
A = np.zeros([Ntot, Ntot])
b = np.zeros(Ntot)

# Define function for the position
def get_position(i, j, Nx, Ny):
    return i * Ny + j

# Define boundary conditions
# Left boundary: x = 0
i = 0
for j in range(Ny):
    n = get_position(i, j, Nx, Ny)
    A[n, n] = 1
    b[n] = 0

# Right boundary: x = x_length
i = Nx-1
for j in range(Ny):
    n = get_position(i, j, Nx, Ny)
    A[n, n] = 1
    b[n] = 0

# Bottom boundary: y = 0
j = 0
for i in range(1, Nx-1):
    n = get_position(i, j, Nx, Ny)
    A[n, n] = 1
    b[n] = 0

# Top boundary: y = y_length
j = Ny-1
for i in range(1, Nx-1):
    n = get_position(i, j, Nx, Ny)
    A[n, n] = 1
    b[n] = 0<jupyter_output><empty_output><jupyter_text>Before we can solve the resulting equation system, we need to define the equations for the interior points.<jupyter_code># Define the factors
factor_x = 1/(dx**2)
factor_y = 1/(dy**2)
factor_cent = 2 * (factor_x + factor_y)

# Loop over the interior points to define A and b
for i in range(1,Nx-1):
    for j in range(1, Ny-1):
        n = get_position(i,j,Nx,Ny)
        A[n,n-Ny] = -factor_x
        A[n,n+Ny] = -factor_x
        A[n,n-1] = -factor_y
        A[n,n+1] = -factor_y
        A[n,n] = factor_cent
        b[n] = 1<jupyter_output><empty_output><jupyter_text>We can finally solve the system with the `numpy.linalg.solve` function.<jupyter_code># Solve with Gaussian elimination
phi = np.linalg.solve(A, b)<jupyter_output><empty_output><jupyter_text>Plot the resulting solution.<jupyter_code># Contour plot of the solution
import matplotlib.pyplot as plt

Z = np.zeros(np.shape(X))
for i in range(Nx):
    for j in range(Ny):
        n = get_position(i, j, Nx, Ny)
        # We need so swap the x and y indices due to the default meshgrid ordering
        Z[j, i] = phi[n]
        
h = plt.contourf(X, Y, Z)
plt.colorbar(h)
plt.xlabel("x axis")
plt.ylabel("y axis")
plt.show()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup