# Laplace's equation for electrostatics

- toc:false
- branch: master
- badges: true
- comments: false
- categories: [mathematics, numerical recipes]
- hide: true

-----
Questions:
- How do I use a finite difference method to calculate derivatives? 
- How do I use the relaxation method to solve Laplace's equation?

Objectives:
- Use the finite difference method to convert Laplace's equation into a set of linear simultaneous equations
- Use the relaxation method to solve Laplace's equation

-----

### The method of finite differences is often used to solve partial differential equations

Consider the two-dimensional Laplace equation for the electric potential $\phi$ subject to appropriate boundary conditions:

\begin{equation}
\frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial y^2} = 0
\end{equation}

Real physical problems are in three dimensions, but we can more easily visualise the method of finite differences - and the extension to three dimensions is straight forward.

The method of finite differences involves dividing the space into a grid of discrete points $[x,y]$ and calculating numerical derivatives or at each of these points. 


 
In this case we consider a 2-dimensional sheet with a fixed voltage $V$ at the top side, and all other sides fixed at $0V$.  But how do we calculate these numerical derivatives?

### A numerical derivative can be calculated using the forward, backward or central difference

The standard definition of a derivative is 

\begin{equation}
\frac{\mathrm{d} f}{\mathrm{d} x} = \lim_{h\to0}\frac{f(x+h)-f(x)}{h}.
\end{equation}

To calculate the derivative numerically we make $h$ very small and calculate

\begin{equation}
\frac{\mathrm{d} f}{\mathrm{d} x} \simeq \frac{f(x+h)-f(x)}{h}.
\end{equation}

This is the <mark>forward difference</mark> because it is measured in the forward direction from $x$.

The <mark>backward difference</mark> is measured in the backward direction from $x$:

\begin{equation}
\frac{\mathrm{d} f}{\mathrm{d} x} \simeq \frac{f(x)-f(x-h)}{h},
\end{equation}

and the <mark>central difference</mark> uses both the forwards and backwards directions around $x$:

\begin{equation}
\frac{\mathrm{d} f}{\mathrm{d} x} \simeq \frac{f(x+\frac{h}{h2})-f(x-\frac{h}{2})}{h},
\end{equation}

> Important: The central difference is in general more accurate than the forward or backward differences, but im all cases we need to carefully converge with respect to the step size $h$

### Second-order derivatives, partial derivatives and the Laplacian can be calculated using the same approach

The second derivative is a derivative of a derivative, and so we can calculate it be applying the first derivative formulas twice. The resulting expression (after application of central differences) is:

\begin{equation}
\frac{\mathrm{d} ^2f}{\mathrm{d} x^2} \simeq \frac{f(x+h)-2f(x)+f(x-h)}{h^2}.
\end{equation}

The extension to partial derivatives is straight-forward:

\begin{equation}
\frac{\partial f}{\partial x} \simeq \frac{f(x+\frac{h}{2},y)-f(x-\frac{h}{2},y)}{h},
\end{equation}

\begin{equation}
\frac{\partial f}{\partial y} \simeq \frac{f(x,y+\frac{h}{2})-f(x,y-\frac{h}{2})}{h},
\end{equation}

\begin{equation}
\frac{\partial ^2f}{\partial x^2} \simeq \frac{f(x+h,y)-2f(x,y)+f(x-h,y)}{h^2},
\end{equation}

\begin{equation}
\frac{\partial ^2f}{\partial y^2} \simeq \frac{f(x,y+h)-2f(x,y)+f(x,y-h)}{h^2}.
\end{equation}

By adding the two equations above, the Laplacian in two dimensions is:

\begin{equation}
\frac{\partial ^2f}{\partial x^2} + \frac{\partial ^2f}{\partial y^2} \simeq \frac{f(x+h,y)+f(x-h,y)+f(x,y+h)+f(x,y-h)-4f(x,y)}{h^2},
\end{equation}

The expression above is known as a <mark>five-point stencil</mark> as it uses five points to calculate the Laplacian.



### The finite difference method turns our partial differential equation into a set of linear simulatenous equation

Returning to our Laplace equation for for the electric potential $\phi$:

\begin{equation}
\frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial y^2} = 0
\end{equation}



The numerical Laplacian can be substituted into the equation above, giving us a set of $n$ simulatenous equations for the $n$ grid points.

\begin{equation}
 \frac{\phi(x+h,y)+\phi(x-h,y)+\phi(x,y+h)+\phi(x,y-h)-4\phi(x,y)}{h^2} = 0,
\end{equation}

where $h$ is the distance between each grid point. 

### To  solve this set of equations we use the relaxation method

To calculate $\phi(x,y)$ we use the relaxation method, also known as the Jacobi method in the context of the Laplace equation.
First we re-arrange the equation above:

\begin{equation}
 \phi(x,y)=\frac{1}{4}\left(\phi(x+h,y)+\phi(x-h,y)+\phi(x,y+h)+\phi(x,y-h)\right).
\end{equation}

<mark> Note that because we set the Laplacian equal to zero in Equation 13 (for this particular example), the $h^2$ term has dropped out of the expression - this might not be the case for other examples. </mark>

This tells us that $\phi(x,y)$ is the average of the surrounding grid points, which can be represented visually as:


 
 Second, we fix $\phi(x,y)$ at the boundaries using the boundary conditions. Third, we guess the initial values of the interior $\phi(x,y)$ points - our guesses do not need to be good, and can be zero.
 
Finally we use Equation 14 to calculate new values of $\phi'(x,y)$ at all points in space. We take these new $\phi'(x,y)$ values and feed them into Equation 14 again to calculate new values. We repeat this iterative process until the $\phi(x,y)$ values converge, and that is our solution.

Convergence can be tested by specifying what the maximum difference should be between iterations. For example, that $\phi'(x,y)-\phi(x,y)< 1e-5$ for all grid points.

### The relaxation and finite difference methods can be applied using the Python skills we have developed

We will now use our Python Skillz to solve Laplace's equation with the boundary conditions outlined above. Let's also imagine that the sheet is 1m along each side and that we want a grid spacing of 1cm. First let's import the libraries we will use:


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

Now let's specify our simulation parameters...


```python
grid_width = 101 # number of grid points for width of simulation
grid_height = 101 # number of grid points for height of simulation
target = 1e-6 # target accuracy to complete convergence
```

...and our boundary conditions:


```python
V_top = 1.0 # top wall voltage, in volts
V_left = 0.0 # left wall voltage
V_right = 0.0 # right wall voltage
V_bottom = 0.0 # bottom wall voltage
```

Next let's create a NumPy array to hold our $\phi(x,y)$ and $\phi'(x,y)$ values:


```python
phi = np.zeros([grid_height,grid_width], float)
phi_prime = np.empty([grid_height,grid_width], float)
```

Now we apply the boundary conditions to our array:


```python
phi[0,:] = V_top
```

Now we write a function which uses the finite difference method to calculate an updated $\phi'(x,y)$ .


```python
def finite_difference(phi):

    for i in range(grid_height): # for each grid point
        for j in range(grid_width):
            if i==0 or i==grid_height-1 or j==0 or j==grid_width-1:
                phi_prime[i,j] = phi[i,j] # if at boundary, keep fixed
            else: # otherwise apply finite difference
                phi_prime[i,j] = (phi[i+1,j]+phi[i-1,j]+phi[i,j+1]+phi[i,j-1]) / 4
    
    return phi_prime
```

Finally let's use the relaxation method. We repeatedly call the function `finite_difference` until all values of $\phi(x,y)$ are converged. 


```python
delta = 1.0 # create delta. It can take any value larger than the target accuracy

while delta > 1e-6:  # keep running the following code until delta < 1e-6
    phi_prime = finite_difference(phi)  # calculate phi_prime
    delta = np.max(np.abs(phi-phi_prime)) # calculate the maximum difference between phi and phi_prime
    phi,phi_prime = phi_prime,phi # Swap phi and phi-prime, ready for the next iteration
```

We can visualise our result using the function `matplotlib.pyplot.imshow` which displays our data as an image:



```python
plt.imshow(phi_prime)
```

This result makes sense: there is a region of high electric potential around the top side of the sheet, where the voltage is fixed at 1V, and regions of low potential around the other three walls. If we would like a colour bar to indicate the $\phi(x,y)$ values across the image then we can use the function `matplotlib.pyplot.contourf` to produce a filled contour plot. Note that this function flips our image (plotting the values held in the array from left to right, top to bottom) so we use `numpy.flip` to achieve the expected result.


```python
plt.contourf(np.flip(phi_prime),levels=np.linspace(0,1,101))
plt.colorbar().set_label("Electric potential")
```

---
Keypoints:
- The finite difference method for numerical derivatives is often used to solve partial differential equations
- A numerical derivative can be calculated using the forward, backward or central difference
- Second-order derivatives, partial derivatives and the Laplacian can be calculated using the same approach
- The finite difference method turns our partial differential equation into a set of linear simulatenous equation
- To solve this set of equations we use the relaxation method
- The relaxation and finite difference methods can be applied using the Python skills we have developed

---

Do [the quick-test](https://nu-cem.github.io/CompPhys/2021/08/02/Finite-Difference-Qs.html).

Back to [Modelling with Partial Differential Equations](https://nu-cem.github.io/CompPhys/2021/08/02/PDEs.html).

---

f04573bb8f10b062d6c1f73b37cb436bc0f24b49
ipynb
Jupyter Notebook
_notebooks/2021-08-02-Finite-Difference.ipynb
NU-CEM/CompPhys
7a7e8ab672797d6db0dbbd673fdc66c6e7bad971
CC-BY-4.0
_notebooks/2021-08-02-Finite-Difference.ipynb
NU-CEM/CompPhys
7a7e8ab672797d6db0dbbd673fdc66c6e7bad971
CC-BY-4.0
2021-10-06T08:11:43.000Z
2021-11-24T13:46:08.000Z
_notebooks/2021-08-02-Finite-Difference.ipynb
NU-CEM/CompPhys
7a7e8ab672797d6db0dbbd673fdc66c6e7bad971
CC-BY-4.0
2021-11-03T10:11:28.000Z
2021-11-03T10:11:28.000Z
Qwen/Qwen-72B
1. YES
2. YES