# 16 Partial Differential Equations — 1 (Students)

## Solving Laplace's or Poisson's equation

**Poisson's equation** for the electric potential $\Phi(\mathbf{r})$ and the charge density $\rho(\mathbf{r})$:

$$
\nabla^2 \Phi(x, y, z) = -4\pi\rho(x, y, z)\\
$$

For a region of space without charges ($\rho = 0$) this reduces to **Laplace's equation**

$$
\nabla^2 \Phi(x, y, z) = 0
$$


Solutions depend on the **boundary conditions**: the value of the potential on the *boundary* (or the electric field normal to the surface, which directly follows from the charge distribution).

### Example: 2D Laplace equation
$$
\frac{\partial^2 \Phi(x,y)}{\partial x^2} + \frac{\partial^2 \Phi(x,y)}{\partial y^2} = 0
$$
("elliptic PDE")

Boundary conditions:
* square area surrounded by wires
* three wires at ground (0 V), one wire at 100 V

## Finite difference algorithm for Poisson's equation
Discretize space on a lattice (2D) and solve for $\Phi$ on each lattice site.

Taylor-expansion of the four neighbors of $\Phi(x, y)$:

\begin{align}
\Phi(x \pm \Delta x, y) &= \Phi(x, y) \pm \Phi_x \Delta x + \frac{1}{2} \Phi_{xx} \Delta x^2 + \dots\\
\Phi(x, y \pm \Delta y) &= \Phi(x, y) \pm \Phi_y \Delta x + \frac{1}{2} \Phi_{yy} \Delta x^2 + \dots\\
\end{align}

Add equations in pairs: odd terms cancel, and **central difference approximation** for 2nd order partial derivatives (to $\mathcal{O}(\Delta^4)$):

\begin{align}
\Phi_{xx}(x,y) = \frac{\partial^2 \Phi}{\partial x^2} & \approx 
  \frac{\Phi(x+\Delta x,y) + \Phi(x-\Delta x,y) - 2\Phi(x,y)}{\Delta x^2} \\
\Phi_{yy}(x,y) = \frac{\partial^2 \Phi}{\partial y^2} &\approx 
  \frac{\Phi(x,y+\Delta y) + \Phi(x,y-\Delta y) - 2\Phi(x,y)}{\Delta y^2}
\end{align}

Take $x$ and $y$ grids of equal spacing $\Delta$: Discretized Poisson equation

$$
\Phi(x+\Delta x,y) + \Phi(x-\Delta x,y) + \Phi(x,y+\Delta y) + \Phi(x,y-\Delta y) - 4\Phi(x,y) = -4\pi\rho(x,y)\,\Delta^2
$$

Defines a system of $N_x \times N_y$ simultaneous algebraic equations for $\Phi_{ij}$ to be solved.

Can be solved directly via matrix approaches (and then is the best solution) but can be unwieldy for large grids.

Alternatively: **iterative solution**:

$$
4\Phi(x,y) = \Phi(x+\Delta x,y) + \Phi(x-\Delta x,y) + \Phi(x,y+\Delta y) + \Phi(x,y-\Delta y) + 4\pi\rho(x,y)\,\Delta^2
$$

Or written for lattice sites $(i, j)$ where 

$$
x = x_0 + i\Delta\quad\text{and}\quad y = y_0 + j\Delta, \quad 0 \leq i,j < N_\text{max}
$$

$$
\Phi_{i,j} = \frac{1}{4}\Big(\Phi_{i+1,j} + \Phi_{i-1,j} + \Phi_{i,j+1} + \Phi_{i,j-1}\Big)
     + \pi\rho(i\Delta, j\Delta) \Delta^2
$$

* Converged solution at $(i, j)$ will be the average potential from the four neighbor sites + charge density contribution.
* *Not a direct solution*: iterate and hope for convergence.

#### Jacobi method
Do not change $\Phi_{i,j}$ until a complete sweep has been completed.

#### Gauss-Seidel method
Immediately use updated new values for $\Phi_{i-1, j}$ and $\Phi_{i, j-1}$ (if starting from $\Phi_{1, 1}$).

Leads to *accelerated convergence* and therefore *less round-off error* (but distorts symmetry of boundary conditions... hopefully irrelevant when converged but check!)

## Example: wire in a box: Solution via relaxation (Gauss-Seidel) 

Solve the box-wire problem on a lattice: The wire at $x=0$ (the $y$-axis) is at 100 V, the other three sides of the box are grounded (0 V).

Note: $\rho=0$ inside the box.


```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#%matplotlib inline
%matplotlib notebook
```

* set boundary conditions
* Implement Gauss-Seidel algorithm
* visualize solution
* does it make sense?
* try higher `Max_iter`


```python
Nmax = 100
Max_iter = 70
Phi = np.zeros((Nmax, Nmax), dtype=np.float64)

# initialize boundaries
# everything starts out zero so nothing special for the grounded wires

# wire at x=0 at 100 V
raise NotImplementedError     

for n_iter in range(Max_iter):
    for xi in range(1, Phi.shape[0]-1):
        for yj in range(1, Phi.shape[1]-1):
            raise NotImplementedError
```


```python
# plot Phi(x,y)
x = np.arange(Phi.shape[0])
y = np.arange(Phi.shape[1])
X, Y = np.meshgrid(x, y)

Z = Phi[X, Y]
```

### Simple wireframe plot 


```python
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(X, Y, Z)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel(r'potential $\Phi$ (V)')
```

(Note: The coordinate system has Y running in the opposite direction compared to the normal coordinate system because numpy arrays "have 0,0 in the upper left hand corner" and we decided to associate rows with X and columns with Y.)

### Surfaces and 2D contours 

Nicer plot (use this code for other projects):



```python
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.plot_wireframe(X, Y, Z, rstride=2, cstride=2, linewidth=0.5, color="gray")
surf = ax.plot_surface(X, Y, Z, cmap=plt.cm.coolwarm, alpha=0.6)
cset = ax.contourf(X, Y, Z, zdir='z', offset=-50, cmap=plt.cm.coolwarm)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel(r'potential $\Phi$ (V)')
ax.set_zlim(-50, 100)
ax.view_init(elev=40, azim=-65)

cb = fig.colorbar(surf, shrink=0.5, aspect=5)
cb.set_label(r"potential $\Phi$ (V)")
```

## Example: Parallel plate capacitor in a box
* two thin sheets of conducting materials at +100 V and -100 V
* width 50, distance 20, box 100
* box is grounded

### Implement solution with Gauss-Seidel
* implement boundary conditions (plates!)
* Gauss-Seidel algorithm
* convergence criterion: change in Frobenius norm (`np.linalg.norm(A)` = $||A|| = \sqrt{\sum_{i,j} a_{ij}^2}$) should be less than tol:
  $$
  || \Phi^{\text{new}} - \Phi^{\text{old}} || < \epsilon
  $$


```python
Nmax = 100
Max_iter = 200
tol = 1e-2
Phi = np.zeros((Nmax, Nmax), dtype=np.float64)
Phi_old = np.zeros_like(Phi)

def set_boundaries(Phi, w=50, d=20):
    # box (edges) always at 0 V
    Phi[:, 0] = Phi[:, -1] = 0
    Phi[0, :] = Phi[-1, :] = 0
    # plate 1 always at +100 V, plate 2 always -100V
    Nx, Ny = Phi.shape
    
    raise NotImplementedError
    
    return Phi
    
set_boundaries(Phi)
Phi_old[:, :] = Phi
print("Starting...")
for n_iter in range(Max_iter):
    if n_iter % 10 == 0:
        print("Iteration {0:5d}".format(n_iter), end="\r")
    for xi in range(1, Nmax-1):
        for yj in range(1, Nmax-1):
            raise NotImplementeError
            
    # enforce boundary conditions
    set_boundaries(Phi)
    
    # check convergence (using the Frobenius matrix norm, delta_norm)
    
    raise NotImplementeError
    
    Phi_old[:, :] = Phi
else:
    print("\nDid not converge in {0} steps, Delta FrobeniusNorm(Phi) = {1} "
          " > {2}".format(n_iter, delta_norm, tol))
```

### Check convergence and Visualize
* check initial results
* increase `Max_iter` if necessary


```python
x = np.arange(Nmax)
y = np.arange(Nmax)
X, Y = np.meshgrid(x, y)

Z = Phi[X, Y]
```


```python
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Y, Z, cmap=plt.cm.coolwarm, alpha=0.3)

cset = ax.contour(X, Y, Z, 20, zdir='z', offset=-100, cmap=plt.cm.coolwarm)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel(r'potential $\Phi$ (V)')
ax.set_zlim(-110, 100)

fig.colorbar(surf, shrink=0.5, aspect=5)
```

## Successive Over Relaxation (SOR)
Accelerate convergence with the scheme

$$
r_{i, j} = \Phi_{i,j}^\text{new} - \Phi_{i, j}^\text{old}\\
\Phi_{i,j}^\text{new} = \Phi_{i,j}^\text{old} + \omega r_{i,j}
$$

Values of $1 \leq \omega \leq 2$ may work well, $\omega > 2$ can lead to numerical instabilities. Experiment!

### Implement SOR with the capacitor problem
* use code above and add SOR
* try $\omega = 1.4$ for a start


```python
# plate capacitor with SOR
raise NotImplementedError
```

#### Check convergence and Visualize
* Find $\omega$ so that it converges in about half the number of steps than with Gauss-Seidel.
* Visualize


```python
x = np.arange(Nmax)
y = np.arange(Nmax)
X, Y = np.meshgrid(x, y)

Z = Phi[X, Y]
```


```python
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Y, Z, cmap=plt.cm.coolwarm, alpha=0.3)

cset = ax.contour(X, Y, Z, 20, zdir='z', offset=-100, cmap=plt.cm.coolwarm)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel(r'potential $\Phi$ (V)')
ax.set_zlim(-110, 100)

fig.colorbar(surf, shrink=0.5, aspect=5)
```

b86bbb052a98977e012c153c8d1dfe43c82e451d
ipynb
Jupyter Notebook
16_PDEs/16_PDEs-1-Students.ipynb
ASU-CompMethodsPhysics-PHY494/PHY494-resources-2019
e6114b49d28df887abe37c8144df8f4ae8cf6419
CC-BY-4.0
16_PDEs/16_PDEs-1-Students.ipynb
ASU-CompMethodsPhysics-PHY494/PHY494-resources-2019
e6114b49d28df887abe37c8144df8f4ae8cf6419
CC-BY-4.0
16_PDEs/16_PDEs-1-Students.ipynb
ASU-CompMethodsPhysics-PHY494/PHY494-resources-2019
e6114b49d28df887abe37c8144df8f4ae8cf6419
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES