This function solves the laplace equation over a square domain with Dirichlet boundary conditions that can be arbitrarily specified (homogeneous or not). In particular, we are interested in solving
\[ \nabla^2 u(x, y) = 0 | (x, y) \in (a, b) \times (a, b) := D \]
with Dirichlet boundary conditions
\[ g(x, y) | (x, y) \in \partial D \]
Input:
Domain begin and end (a and b), size of mesh, function for computing stencil, and a function for the Dirichlet boundary conditions.
Output:
A 2D matrix that represents the solution along the interior of the domain.
Usage/Example:
First of all, I am using a few type conventions that will prove useful in making our definitions more terse and readable. They are as follows:
The following code is inside of a rather large main function, so I’ll break it up part by part to make the explanation more clear.
First, we generate the basic parameters that define the type of analysis we wish to perform. Here we define the domain, the function evaluated at the boundary, and the size of the mesh we wish to use.
Then, we need to define what type of stencil we will be using for analysis. Here I define two stencils: the five and nine point equivalent discrete forms of the 2D laplace operator. Notice that the stencil points also include multipliers. To define a custom stencil, one must only create a new lambda function here and provide it to the solver.
Then, to solve the equation, we just need to provide these parameters to the solver function.
Output:
Since I required the upper boundary (where \( y = 0 \) since I am using matrix row as \(x\)) to have a value of \( 5.0 \), then we should expect to see a smooth distribution of values that are resemble a Gaussian distribution centered about the vertical that decreases in overall magnitude as we reach the bottom of the domain. Even with a small mesh for more terse output, it is apparent that the gaussian symmetry exists for both solutions:
Implementation/Code:
To implement this solver in a general way, I use lambda functions as parameters for both the Dirichlet boundary conditions and the stencil desired to be used. Both the right-hand-side and discrete Laplace operator are generated on the fly using these parameters. I make use of an additional method I added to my matrix class wherein I may flatten a matrix into a column vector according to lexigraphical ordering, or reshape it into an \(m\) by \(n\) matrix if it represents a column vector. This way I can mutate a matrix in a way that makes sense, then flatten it into the ordering required for our linear system. The comments in my code should help explain further: