Relaxation methods are a significant set of numerical analysis techniques utilized primarily in the fields of electrostatics and fluid dynamics.

Let ( Phi ) represent a function of ( x ) and ( y ), with ( A, B, ldots, G ) also being functions of these variables. Any second-order partial differential equation (PDE) with ( Phi ) as its solution can be expressed in a specific format. The subscripts denote partial derivatives; for instance, ( frac{partial^2 Phi}{partial x} ). Here, we see ( 2B ) instead of merely ( B ) because ( frac{partial^2 Phi}{partial x partial y} = frac{partial^2 Phi}{partial y partial x} ). Our primary focus will be on PDEs where ( B^2 - AC ) is strictly negative, classifying them as elliptic equations.

Among the simplest elliptic equations is Poisson's equation, ( nabla^2 Phi = f(x,y) ), with the special scenario ( nabla^2 Phi = 0 ) termed Laplace's equation. The operator ( nabla^2 ) is recognized as the Laplacian:

Given that ( B = 0 ) and ( A = C = 1 ), it’s evident that both Laplace's and Poisson's equations are elliptic. Poisson's equation is crucial in electrostatics, emerging from Gauss's law: since ( nabla cdot mathbf{E} = frac{rho}{epsilon_0} ) and ( mathbf{E} = -nabla Phi ), it follows that ( nabla cdot mathbf{E} = -nabla^2 Phi ), where ( rho ) denotes the charge density function. Consequently, Laplace's equation is applicable in regions devoid of charges, such as between capacitor plates.

Moreover, Laplace's equation is relevant in fluid dynamics. If the fluid's velocity field ( mathbf{V} ) is irrotational, then ( nabla times mathbf{V} = 0 ), implying ( mathbf{V} ) is conservative, and thus there exists a function ( Phi ) such that ( mathbf{V} = nabla Phi ). Taking the divergence on both sides yields ( nabla cdot mathbf{V} = nabla^2 Phi ). By employing the continuity equation, we derive ( -frac{partial rho}{partial t} = nabla^2 Phi ), with ( rho ) as the fluid's density. If the fluid is incompressible, ( rho ) remains constant, leading to ( nabla^2 Phi = 0 ). This potential is measured in ( m^2/s ).

In this discussion, we aim to solve Laplace's equation in spatial regions bounded by surfaces where ( Phi ) is specified. Such problems are known as Dirichlet problems concerning Laplace's equation.

Analytically solving Dirichlet problems can be labor-intensive, particularly with simple boundary conditions; for many practical scenarios, pursuing an analytical solution becomes impractical. Therefore, various numerical approximation schemes have been developed, among which relaxation techniques stand out as a crucial category. Broadly speaking, a relaxation scheme estimates the solution at each point by expressing it as a straightforward function of its neighboring points.

To construct a relaxation algorithm for Laplace's equation, we can impose a square grid onto the area of interest, with points spaced apart by a distance ( h ). We can then estimate ( Phi ) at ( (x,y) ) by averaging it over the neighboring points: ( (x+h,y) ), ( (x,y+h) ), ( (x-h,y) ), and ( (x,y-h) ). This approximation's accuracy can be validated using Taylor's theorem. We begin by writing the Taylor expansions of ( Phi ) for the closest neighbors on the grid:

The notation ( O(h^4) ) refers to terms where ( h ) appears raised to the fourth power or higher. By summing these four equations and rearranging ( Phi(x,y) ) to the left side, we achieve:

The term associated with ( h^2 ) corresponds to the Laplacian; since ( nabla^2 Phi = f ), this term simplifies to ( h^2 f ). Moreover, as ( h ) approaches 0, the higher powers of ( h ) diminish at a much quicker rate than the ( h^2 ) term. Thus, when ( h ) is exceedingly small compared to the problem's length scale, the higher power terms become negligible, and the following approximation holds:

Now, shifting our focus to the relaxation algorithm itself, we will only consider boundaries that are parallel to the ( x ) and ( y ) axes, as illustrated below, where the black dots indicate grid points spaced 1 centimeter apart.

We will simplify by assuming ( f = 0 ).

I will utilize MATLAB for my examples, as it is particularly suited for these tasks. If you lack MATLAB but wish to experiment with the code, Scilab, a free and open-source MATLAB alternative, can be used and is compatible with basic MATLAB scripts. It is important to note that in MATLAB, the coordinate (1,1) refers to the top left element of an array, and ( (i,j) ) corresponds to the element in row ( i ) and column ( j ).

The procedure is as follows:

• Initialize an ( (M+1) times (N+1) ) rectangular array named phi(i,j), where ( M ) and ( N ) are the width and height of the modeled area divided by the grid spacing. For instance, in the previous illustration, the width and height are 4 centimeters, and the grid spacing is 1 centimeter, yielding ( M = N = 4 ). Assign boundary values to the relevant points.
• For all points ( (i,j) ) that are not on the boundary, compute ( text{PHI}(i,j) = 0.25 times [text{PHI}(i+1,j) + text{PHI}(i-1,j) + text{PHI}(i,j+1) + text{PHI}(i,j-1)] ).
• Terminate after achieving the desired number of iterations or level of accuracy.
• If desired, both MATLAB and Scilab offer built-in functions to plot the results.

The subsequent images illustrate the results of applying relaxation techniques to Example 1, with grid spacing adjusted to a more practical 0.125 centimeters. To keep the article concise, I have placed the code on Pastebin, where the comments may be helpful for experimentation. The code is compatible with Scilab.

These visuals further exemplify the screening effect: a region surrounded by a grounded conductor is entirely shielded from external potentials. However, complete enclosure is not necessary to ensure effective protection from interference; this principle underpins the functionality of a Faraday cage.

Consider a one-meter long plate maintained at 10kV, with a Faraday cage constructed from centimeter-thick bars, spaced 10 centimeters apart, positioned 50 centimeters from the plate. The code I developed for this scenario can be found here.

First, let’s observe the potential without the Faraday cage.

As anticipated, the potential decreases approximately linearly as one moves further from the plate, approaching a perfect linearity if the plate were infinitely long in the ( y )-direction. Now, let’s introduce the Faraday cage.

The presence of the cage significantly impacts the potential, despite not perfectly isolating the two regions. Next, we will explore an example from fluid dynamics.

Consider a 1 meter wide L-shaped pipe, with sides measuring 3 meters, where the ends are set at 200 m²/s and -200 m²/s, while the walls maintain 0 m²/s. Our goal is to determine the flow within the pipe. Since ( mathbf{V} = nabla Phi ), the flow direction aligns with increasing ( Phi ), establishing that the bottom left boundary at -200 m²/s serves as the pipe's inlet, and the top left at 200 m²/s acts as the outlet.

The code for the relaxation process can be found here. Below is the resultant graph:

An intriguing aspect of this potential profile is its flatness near the corner, which persists even under high potential values at both inlet and outlet. This behavior is connected to the skew-symmetry of boundary conditions about the line ( y = 300 - x ). We can interpret this graph as indicating that a test particle introduced into the fluid at the midpoint of the bottom left corner (e.g., at ( x = 0, y = 50 )) initially moves swiftly, slows as it reaches the corner, and accelerates again toward the outlet at the top right. A streamline plot illustrates some of these trajectories, though it does not depict the particle's speed at each trajectory point:

Another noteworthy behavior is that in the horizontal section of the pipe, test particles are directed toward the wall, while in the vertical section, they are pushed away. Although the streamline plot does not illustrate this, particles do not cease movement upon contacting the wall; instead, they travel along the wall in a direction roughly parallel to it. Similarly, particles do not originate at the vertical walls, as the streamlines diverging from the walls represent particles previously flowing along the walls being pushed away.

This solution represents the least precise of the three I have presented. Specifically, the streamline shape is more sensitive to grid spacing than the mesh plot's potential shape. Unfortunately, I am currently using a 2011 MacBook Air until I can afford a new computer, so I couldn't reduce the spacing as much as desired. However, if you have access to a more powerful computer, I highly recommend experimenting with the MATLAB code settings, particularly trying a separation of 0.1 centimeters instead of the default 5 centimeters.