# Laplace equation solution for electrical potential

I'm trying to simulate a lightning in a programming language and I started to read how it can be done, and I found that it could be done by using the Laplacian Growing Model. My experience with mathematics is limited.

I think I understood the algorithm described here in section 4.1. My problem is that I don't know how to solve the Laplacian equation for the lightning model seen at (b) Lightning configuration.

I tried to catch something up from here but that configuration is a rectangle and the lightning configuration seems like a triangle or maybe a circular arc.

I know my boundary conditions would be the potential at the origin and the potential at the ground but I really don't know how to translate this in mathematical language and how to resolve the laplace for it.

Can you help me with a solution (along with an explanation) or a tutorial "for dummies"?

-
Do you have a preferred language? Would matlab be okay? – Omnomnomnom Jul 28 '13 at 20:20
@Omnomnomnom I want to implement it in Java, but I need the theoretical solution for the problem, but I think I can understand some Matlab too. But my question is theoretical, language-adnostic. – DaJackal Jul 29 '13 at 10:51

Here's a description of how to solve the laplacian equation for an arbitrary $2$-D grid with boundary conditions using the so-called "iteration method":

## Initialization:

Create a matrix $M$ representing each point on the grid, with one entry per "pixel". So, for example, to work with a $7\times 7$ grid, we would set $$M = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}$$ We then set permanent boundary conditions that won't be changed as we improve our approximation. To get something kind of circular, we could set $$M = \begin{bmatrix} \mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1\\ \mathbf1 & \mathbf1 & 0 & 0 & 0 & \mathbf1 & \mathbf1\\ \mathbf1 & 0 & 0 & 0 & 0 & 0 & \mathbf1\\ \mathbf1 & 0 & 0 & \mathbf0 & 0 & 0 & \mathbf1\\ \mathbf1 & 0 & 0 & 0 & 0 & 0 & \mathbf1\\ \mathbf1 & \mathbf1 & 0 & 0 & 0 & \mathbf1 & \mathbf1\\ \mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 \end{bmatrix}$$ The bold entries remain constant, while the rest will update as we improve our approximation.

We then normally put in some non-zero initial guess for each point to be computed, just a very rough approximation. For instance, in this case we could put $$M = \begin{bmatrix} \mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1\\ \mathbf1 & \mathbf1 & 0.66 & 0.66 & 0.66 & \mathbf1 & \mathbf1\\ \mathbf1 & 0.66 & 0.33 & 0.33 & 0.33 & 0.66 & \mathbf1\\ \mathbf1 & 0.66 & 0.33 & \mathbf0 & 0.33 & 0.66 & \mathbf1\\ \mathbf1 & 0.66 & 0.33 & 0.33 & 0.33 & 0.66 & \mathbf1\\ \mathbf1 & \mathbf1 & 0.66 & 0.66 & 0.66 & \mathbf1 & \mathbf1\\ \mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 \end{bmatrix}$$

We are now ready to begin the iterative process

## Iteration:

What we use now (and will justify later, maybe) is the fact that on a $2$D plane where $\phi$ is a function satisfying $$\nabla^2\phi=\frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}=0$$ the following is approximately true: Suppose that $\phi(x_0,y_0)=\phi_0$ Let $\Delta x$ and $\Delta y$ be "sufficiently small" and define: $$\phi_{x^+}=\phi(x_0+\Delta x,y_0)\\ \phi_{x^-}=\phi(x_0-\Delta x,y_0)\\ \phi_{y^+}=\phi(x_0,y_0+\Delta y)\\ \phi_{y^-}=\phi(x_0,y_0-\Delta y)$$ Then we have $$\phi_0 \approx \frac14(\phi_{x^+}+\phi_{x^-}+\phi_{y^+}+\phi_{y^-})$$ In other words, we could calculate the potential at a point by averaging the potential at the four neighboring points. We do exactly that, using the improved values as we go.

For our matrix $M$, let's begin with the first changeable point in the upper-left corner.

\begin{align} \phi_{2,3} &=\frac14(\phi_{2,4}+\phi_{2,1}+\phi_{1,3}+\phi_{3,3}) \\&= 0.25\cdot(0.66 + 1 + 1 + 0.33) \\&= 0.7475 \end{align}

Making the new matrix

$$M = \begin{bmatrix} \mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1\\ \mathbf1 & \mathbf1 & 0.7475 & 0.66 & 0.66 & \mathbf1 & \mathbf1\\ \mathbf1 & 0.66 & 0.33 & 0.33 & 0.33 & 0.66 & \mathbf1\\ \mathbf1 & 0.66 & 0.33 & \mathbf0 & 0.33 & 0.66 & \mathbf1\\ \mathbf1 & 0.66 & 0.33 & 0.33 & 0.33 & 0.66 & \mathbf1\\ \mathbf1 & \mathbf1 & 0.66 & 0.66 & 0.66 & \mathbf1 & \mathbf1\\ \mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 \end{bmatrix}$$ In fact, if we know that we will be using a symmetrical configuration, we could make things more efficient by updating all points that will have the same outcome. That is, we could update $M$ to get

$$M = \begin{bmatrix} \mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1\\ \mathbf1 & \mathbf1 & 0.7475 & 0.66 & 0.7475 & \mathbf1 & \mathbf1\\ \mathbf1 & 0.7475 & 0.33 & 0.33 & 0.33 & 0.7475 & \mathbf1\\ \mathbf1 & 0.66 & 0.33 & \mathbf0 & 0.33 & 0.66 & \mathbf1\\ \mathbf1 & 0.7475 & 0.33 & 0.33 & 0.33 & 0.7475 & \mathbf1\\ \mathbf1 & \mathbf1 & 0.7475 & 0.66 & 0.7475 & \mathbf1 & \mathbf1\\ \mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 \end{bmatrix}$$

We would then proceed to the next point $$\phi_{2,4}=0.25\cdot(0.7475+0.7475+1+0.33)= 0.5194$$ And update the matrix accordingly

$$M = \begin{bmatrix} \mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1\\ \mathbf1 & \mathbf1 & 0.7475 & 0.5194 & 0.7475 & \mathbf1 & \mathbf1\\ \mathbf1 & 0.7475 & 0.33 & 0.33 & 0.33 & 0.7475 & \mathbf1\\ \mathbf1 & 0.5194 & 0.33 & \mathbf0 & 0.33 & 0.5194 & \mathbf1\\ \mathbf1 & 0.7475 & 0.33 & 0.33 & 0.33 & 0.7475 & \mathbf1\\ \mathbf1 & \mathbf1 & 0.7475 & 0.5194 & 0.7475 & \mathbf1 & \mathbf1\\ \mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 &\mathbf1 \end{bmatrix}$$

Proceed in such a fashion until each changeable value is updated once (in this case, that will take $4$ steps if you take advantage of symmetry, $20$ if you make no such simplifications), then repeat. Do this enough, and eventually the updates won't significantly change the value. When that happens, you've reached the final approximation and can stop the process. You now have an approximate solution to the Laplacian equation, given the boundary conditions of your matrix.

A similar method to this one is the relaxation method, which is faster for systems with relatively few pixels and slightly more accurate.

EDIT:

I found this derivation of the iteration method (apparently, also known as the finite difference method) if you want a mathematical justification for the process.

-
Let me know if the method is clear as explained, and let me know if you want the mathematical justification – Omnomnomnom Jul 29 '13 at 17:32
Thank you very much @Omnomnomnom. At the first glance, I think I udnerstood your solution, I can't wait to implement it. – DaJackal Jul 30 '13 at 7:00

Hmm.. maybe I am 2 years late to the party, but what the heck.

You can also pose this as a Global Linear Optimization. Setting up a matrix $L_2$ norm minimization problem: $$\min({\bf v})\left\{{\bf C}\|{\bf v}\|_2 + \|{\bf L(v+d)}\|_2\right\}$$ where $\bf C$ contains a binary diagonal matrix mask like the one nomnom called "boundary condition", the $\bf d$ is a vectorization of boundary condition values. So $\bf v+d$ is the result, and the $\bf L$ is a matrix representation of the laplacian operator, typically some suitable discrete filter. First term punishes changes to boundary conditions. Second term punishes laplacian of result from differing from 0.

So why bother doing all this instead of just taking the first numerical scheme we can find? It is because it is very convenient framework if we want to attack bigger and more involving problems. There is no limit as to how many terms we can have and we can have any constraints we like as long as we can express them as matrix operations on the data. Also as there exist so many numerical methods to solve linear optimizations we are not restricted to one specific numerical method but can use the whole arsenal of linear solving techniques.

-