Relaxation methods represent a significant category of numerical analysis strategies used in fields such as electrostatics and fluid dynamics. This article aims to provide an overview of these techniques.

Let ( phi ) be a function of ( x ) and ( y ), with ( A, B, ldots, G ) also defined as functions of ( x ) and ( y ). Any second-order partial differential equation (PDE) that has ( phi ) as a solution can be expressed in the following manner:

The subscripts indicate partial derivatives; for instance, ( frac{partial phi}{partial x} ). The presence of ( 2B ) instead of ( B ) is due to the equality ( frac{partial^2 phi}{partial x partial y} = frac{partial^2 phi}{partial y partial x} ). We will focus on PDEs where ( B^2 - AC ) is negative, categorizing these as elliptic equations.

A fundamental elliptic equation is Poisson's equation, ( nabla^2 phi = f(x, y) ), with the specific case of ( nabla^2 phi = 0 ) being referred to as Laplace's equation. The operator ( nabla^2 ) is known as the Laplacian:

Notably, with ( B = 0 ) and ( A = C = 1 ), both Laplace's and Poisson's equations are elliptic. The significance of Poisson's equation in electrostatics arises from Gauss's law: since ( nabla cdot mathbf{E} = rho/epsilon_0 ) and ( mathbf{E} = -nabla phi ), we derive ( nabla cdot mathbf{E} = -nabla^2 phi ), where ( rho ) is the charge density. Hence, Laplace's equation is applicable in regions void of charge, such as between capacitor plates.

In fluid dynamics, Laplace's equation also plays a role. If a fluid's velocity field ( mathbf{V} ) is irrotational, then ( nabla times mathbf{V} = 0 ), indicating that ( mathbf{V} ) is conservative, leading to the existence of a potential function ( phi ) such that ( mathbf{V} = nabla phi ). By taking the divergence of both sides, we obtain ( nabla cdot mathbf{V} = nabla^2 phi ). According to the continuity equation, ( -frac{partial rho}{partial t} = nabla^2 phi ), where ( rho ) represents fluid density. If the fluid is incompressible, then ( phi ) remains constant, resulting in ( nabla^2 phi = 0 ).

This article will focus on solving Laplace’s equation within regions bounded by surfaces where ( phi ) is predefined. Such problems are termed Dirichlet problems for Laplace’s equation.

Analytically tackling Dirichlet problems can often be laborious, particularly for simple boundary conditions, and for many realistic scenarios, pursuing analytical solutions may not be practical. Thus, various numerical approximation techniques are available for solving Dirichlet problems, with relaxation methods forming a vital subset of these techniques. A relaxation approach generally approximates the solution at a point by expressing it as a function of its surrounding points.

To create a relaxation algorithm for Laplace's equation, one can impose a square grid over the area of interest, where points are spaced a distance ( h ) apart. The value of ( phi ) at the point ( (x,y) ) can then be approximated as the average of the neighboring points: ( (x+h,y) ), ( (x,y+h) ), ( (x-h,y) ), and ( (x,y-h) ). This approximation can be validated using Taylor's theorem, where one first derives the Taylor expansions of ( phi ) for the nearest grid points:

The notation ( O(h^4) ) denotes terms in which ( h ) is raised to the fourth power or higher. By combining these four equations and isolating ( phi(x,y) ), we arrive at:

The term associated with ( h^2 ) corresponds to the Laplacian, and since ( nabla^2 phi = f ), this yields ( h^2 f ). Furthermore, as ( h ) approaches zero, the higher-order terms vanish significantly faster than the ( h^2 ) term. Therefore, for small ( h ) relative to the problem's length scale, the approximation:

is valid.

Now, shifting focus to the relaxation algorithm, we will restrict our discussion to boundaries aligned with the ( x ) and ( y ) axes. For instance, in the illustration below, the black dots signify grid points spaced 1 centimeter apart.

More intricate geometries involving curves or non-right angles fall outside the scope of this introductory discussion. For simplicity, we will also assume ( f = 0 ).

MATLAB will be utilized for the examples as it is exceptionally suited for such tasks. If MATLAB is unavailable, Scilab, a free and open-source alternative compatible with basic MATLAB code, can be used. It is important to note that in MATLAB, the coordinate (1,1) refers to the top left element of an array, while ( (i,j) ) denotes the element located in row ( i ) and column ( j ).

The procedure is outlined as follows:

1. Initialize a rectangular array ( text{phi}(i,j) ) of dimensions ( (M+1) times (N+1) ) where ( M ) and ( N ) are the width and height of the modeled area divided by the grid spacing. For example, in the earlier illustration, with a width and height of 4 centimeters and a grid spacing of 1 centimeter, ( M = N = 4 ). Assign the boundary values to the appropriate points.
2. For all points ( (i,j) ) not on the boundary, update ( 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)] ).
3. Cease iterations once the desired accuracy or iteration count is reached.
4. If desired, MATLAB and Scilab offer built-in functions for visualizing the results.

The following images showcase the results from applying relaxation techniques to Example 1, with the grid spacing adjusted to a practical 0.125 centimeters. To avoid clutter, the code is available on Pastebin, where comments can assist in experimentation. The code is also compatible with Scilab.

These images illustrate the screening effect: a region enclosed by a grounded conductor remains isolated from external potentials. Notably, complete enclosure is not required to ensure effective protection from interference, which is the foundational principle of a Faraday cage.

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

First, let's visualize the scenario without the Faraday cage.

As anticipated, the potential diminishes approximately linearly with increasing distance from the plate, converging to a perfect linearity for infinitely long plates in the y-direction. Now, let’s introduce the cage.

It is evident that the cage significantly alters the potential distribution, even without fully separating the regions. Next, we will explore an example from fluid dynamics.

In this scenario, the ends of an L-shaped pipe, measuring 1 meter wide and 3 meters long, are maintained at velocities of 200 m²/s and -200 m²/s, while the walls are at 0 m²/s. Our goal is to determine the flow in the pipe. Given that ( mathbf{V} = nabla phi ), the velocity aligns with increasing ( phi ), indicating that the bottom left boundary at -200 m²/s serves as the inlet, and the top left at 200 m²/s functions as the outlet.

The relaxation code for this example is available here, along with the resulting graph:

An intriguing characteristic of this potential is its flatness near the corner, even with significant potential values at both the inlet and outlet. This phenomenon is linked to the skew-symmetry of the boundary conditions around the line ( y = 300 - x ). Interpreting this graph suggests that a test particle entering the fluid at the center of the bottom left corner (e.g., at ( x=0, y=50 )) initially moves rapidly, experiences a delay upon reaching the corner, and then accelerates toward the outlet in the top right. A streamline plot can illustrate some of the trajectories, although it does not depict particle speed at each point along the path:

Another notable 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. While the streamline plot does not capture this, particles do not halt upon contact with the wall; instead, they travel along the wall in a path approximately parallel to it. Similarly, particles do not emerge from the vertical walls; the streamlines pointing away represent particles that were initially flowing along the walls before being redirected.

This solution is the least accurate among the three presented. Notably, the shape of the streamlines is more sensitive to grid spacing than the mesh plot of the potential. Currently, I am limited to my 2011 MacBook Air until I can acquire a new computer, which restricts my ability to reduce spacing as desired. However, if you possess a more powerful machine, I highly recommend experimenting with the MATLAB code settings, particularly trying a grid spacing of 0.1 centimeters instead of the default 5 centimeters.