Nine-point stencil

Last updated

In numerical analysis, given a square grid in two dimensions, the nine-point stencil of a point in the grid is a stencil made up of the point itself together with its eight "neighbors". It is used to write finite difference approximations to derivatives at grid points. It is an example for numerical differentiation. This stencil is often used to approximate the Laplacian of a function of two variables.

Contents

An illustration of the nine-point stencil in two dimensions. 9-point stencil.png
An illustration of the nine-point stencil in two dimensions.

Motivation

If we discretize the 2D Laplacian by using central-difference methods, we obtain the commonly used five-point stencil, represented by the following convolution kernel:

Starting from homogeneous initial conditions with only one point-like perturbation, the correct growth process would yield a circle (left). If growth is faster along the coordinate axes, which is caused by the anisotropic effect of central difference discretization, the circle would turn into star-shaped structure, as the errors propagate (right): . Anisotropic cd.png
Starting from homogeneous initial conditions with only one point-like perturbation, the correct growth process would yield a circle (left). If growth is faster along the coordinate axes, which is caused by the anisotropic effect of central difference discretization, the circle would turn into star-shaped structure, as the errors propagate (right): .

Even though it is simple to obtain and computationally lighter, the central difference kernel possess an undesired intrinsic anisotropic property, since it doesn't take into account the diagonal neighbours. This intrinsic anisotropy poses a problem when applied on certain numerical simulations or when more accuracy is required, by propagating the Laplacian effect faster in the coordinate axes directions and slower in the other directions, thus distortiong the final result. [1]

This drawback calls for finding better methods for discretizing the Laplacian, reducing or eliminating the anisotropy.

Implementation

The two most commonly used isotropic nine-point stencils are displayed below, in their convolution kernel forms. They can be obtained by the following formula: [2] [3]

The first one is known by Oono-Puri, [4] [5] [6] [7] [8] and it is obtained when γ=1/2. [2]

The second one is known by Patra-Karttunen or Mehrstellen, [1] [7] [8] [9] [10] and it is obtained when γ=1/3. [2]

Both are isotropic forms of discrete Laplacian, [8] and in the limit of small Δx, they all become equivalent, [11] as Oono-Puri being described as the optimally isotropic form of discretization, [8] displaying reduced overall error, [2] and Patra-Karttunen having been systematically derived by imposing conditions of rotational invariance, [9] displaying smallest error around the origin. [2]

Desired anisotropy

On the other hand, if controlled anisotropic effects are a desired feature, when solving anisotropic diffusion problems for example, it is also possible to use the 9-point stencil combined with tensors to generate them.

Consider the laplacian in the following form:

Where c is just a constant coefficient. Now if we replace c by the 2nd rank tensor C:

Where c1 is the constant coefficient for the principal direction in x axis, and c2 is the constant coefficient for the secondary direction in y axis. In order to generate anisotropic effects, c1 and c2 must be different.

By multiplying it by the rotation matrix Q, we obtain C', allowing anisotropic propagations in arbitrary directions other than the coordinate axes. [12] [13]

The angle can be obtained by generating a vector field in order to orientate the pattern as desired. [13] Then:

Or, for different anisotropic effects using the same vector field [14]

It is important to note that, regardless of the values of , the anisotropic propagation will occur parallel to the secondary direcion c2 and perpendicular to the principal direction c1: [15] . The resulting convolution kernel is as follows [13]

If, for example, c1=c2=1, the cxy component will vanish, resulting in the simple five-point stencil, rendering no controlled anisotropy.

If c2>c1 and =0, the anisotropic effects will be more pronounced in the vertical axis.

if c2>c1 and =45 degrees, the anisotropic effects will be more pronounced in the upper-right / lower-left diagonal.

Related Research Articles

<span class="mw-page-title-main">Spherical coordinate system</span> 3-dimensional coordinate system

In mathematics, a spherical coordinate system is a coordinate system for three-dimensional space where the position of a point is specified by three numbers: the radial distance of that point from a fixed origin; its polar angle measured from a fixed polar axis or zenith direction; and the azimuthal angle of its orthogonal projection on a reference plane that passes through the origin and is orthogonal to the fixed axis, measured from another fixed reference direction on that plane. When radius is fixed, the two angular coordinates make a coordinate system on the sphere sometimes called spherical polar coordinates.

<span class="mw-page-title-main">Laplace operator</span> Differential operator

In mathematics, the Laplace operator or Laplacian is a differential operator given by the divergence of the gradient of a scalar function on Euclidean space. It is usually denoted by the symbols , (where is the nabla operator), or . In a Cartesian coordinate system, the Laplacian is given by the sum of second partial derivatives of the function with respect to each independent variable. In other coordinate systems, such as cylindrical and spherical coordinates, the Laplacian also has a useful form. Informally, the Laplacian Δf (p) of a function f at a point p measures by how much the average value of f over small spheres or balls centered at p deviates from f (p).

In mathematics, a Gaussian function, often simply referred to as a Gaussian, is a function of the base form

<span class="mw-page-title-main">Inverse trigonometric functions</span> Inverse functions of the trigonometric functions

In mathematics, the inverse trigonometric functions are the inverse functions of the trigonometric functions. Specifically, they are the inverses of the sine, cosine, tangent, cotangent, secant, and cosecant functions, and are used to obtain an angle from any of the angle's trigonometric ratios. Inverse trigonometric functions are widely used in engineering, navigation, physics, and geometry.

In mathematics, a function defined on an inner product space is said to have rotational invariance if its value does not change when arbitrary rotations are applied to its argument.

In numerical linear algebra, a Givens rotation is a rotation in the plane spanned by two coordinates axes. Givens rotations are named after Wallace Givens, who introduced them to numerical analysts in the 1950s while he was working at Argonne National Laboratory.

In linear algebra, linear transformations can be represented by matrices. If is a linear transformation mapping to and is a column vector with entries, then

<span class="mw-page-title-main">Vector fields in cylindrical and spherical coordinates</span> Vector field representation in 3D curvilinear coordinate systems

Note: This page uses common physics notation for spherical coordinates, in which is the angle between the z axis and the radius vector connecting the origin to the point in question, while is the angle between the projection of the radius vector onto the x-y plane and the x axis. Several other definitions are in use, and so care must be taken in comparing different sources.

<span class="mw-page-title-main">Pedal curve</span> Curve generated by the projections of a fixed point on the tangents of another curve

In mathematics, a pedal curve of a given curve results from the orthogonal projection of a fixed point on the tangent lines of this curve. More precisely, for a plane curve C and a given fixed pedal pointP, the pedal curve of C is the locus of points X so that the line PX is perpendicular to a tangent T to the curve passing through the point X. Conversely, at any point R on the curve C, let T be the tangent line at that point R; then there is a unique point X on the tangent T which forms with the pedal point P a line perpendicular to the tangent T – the pedal curve is the set of such points X, called the foot of the perpendicular to the tangent T from the fixed point P, as the variable point R ranges over the curve C.

In linear algebra, a rotation matrix is a transformation matrix that is used to perform a rotation in Euclidean space. For example, using the convention below, the matrix

<span class="mw-page-title-main">Cissoid</span> Plane curve constructed from two other curves and a fixed point

In geometry, a cissoid ( is a plane curve generated from two given curves C1, C2 and a point O. Let L be a variable line passing through O and intersecting C1 at P1 and C2 at P2. Let P be the point on L so that Then the locus of such points P is defined to be the cissoid of the curves C1, C2 relative to O.

<span class="mw-page-title-main">Transverse isotropy</span>

A transversely isotropic material is one with physical properties that are symmetric about an axis that is normal to a plane of isotropy. This transverse plane has infinite planes of symmetry and thus, within this plane, the material properties are the same in all directions. Hence, such materials are also known as "polar anisotropic" materials. In geophysics, vertically transverse isotropy (VTI) is also known as radial anisotropy.

<span class="mw-page-title-main">Osculating circle</span> Circle of immediate corresponding curvature of a curve at a point

In differential geometry of curves, the osculating circle of a sufficiently smooth plane curve at a given point p on the curve has been traditionally defined as the circle passing through p and a pair of additional points on the curve infinitesimally close to p. Its center lies on the inner normal line, and its curvature defines the curvature of the given curve at that point. This circle, which is the one among all tangent circles at the given point that approaches the curve most tightly, was named circulus osculans by Leibniz.

Spatial rotations in three dimensions can be parametrized using both Euler angles and unit quaternions. This article explains how to convert between the two representations. Actually this simple use of "quaternions" was first presented by Euler some seventy years earlier than Hamilton to solve the problem of magic squares. For this reason the dynamics community commonly refers to quaternions in this application as "Euler parameters".

<span class="mw-page-title-main">Plane stress</span> When the stress vector within a material is zero across a particular plane

In continuum mechanics, a material is said to be under plane stress if the stress vector is zero across a particular plane. When that situation occurs over an entire element of a structure, as is often the case for thin plates, the stress analysis is considerably simplified, as the stress state can be represented by a tensor of dimension 2. A related notion, plane strain, is often applicable to very thick members.

<span class="mw-page-title-main">Rotation of axes in two dimensions</span> Transformation of coordinates through an angle

In mathematics, a rotation of axes in two dimensions is a mapping from an xy-Cartesian coordinate system to an x′y′-Cartesian coordinate system in which the origin is kept fixed and the x′ and y′ axes are obtained by rotating the x and y axes counterclockwise through an angle . A point P has coordinates (x, y) with respect to the original system and coordinates (x′, y′) with respect to the new system. In the new coordinate system, the point P will appear to have been rotated in the opposite direction, that is, clockwise through the angle . A rotation of axes in more than two dimensions is defined similarly. A rotation of axes is a linear map and a rigid transformation.

In geometry, various formalisms exist to express a rotation in three dimensions as a mathematical transformation. In physics, this concept is applied to classical mechanics where rotational kinematics is the science of quantitative description of a purely rotational motion. The orientation of an object at a given instant is described with the same tools, as it is defined as an imaginary rotation from a reference placement in space, rather than an actually observed rotation from a previous placement in space.

The direct-quadrature-zerotransformation or zero-direct-quadraturetransformation is a tensor that rotates the reference frame of a three-element vector or a three-by-three element matrix in an effort to simplify analysis. The DQZ transform is the product of the Clarke transform and the Park transform, first proposed in 1929 by Robert H. Park.

<span class="mw-page-title-main">Composite laminate</span>

In materials science, a composite laminate is an assembly of layers of fibrous composite materials which can be joined to provide required engineering properties, including in-plane stiffness, bending stiffness, strength, and coefficient of thermal expansion.

References

  1. 1 2 3 Patra, Michael; Karttunen, Mikko (2006). "Stencils with isotropic discretization error for differential operators". Numerical Methods for Partial Differential Equations. 22 (4): 3–7. doi:10.1002/num.20129. S2CID   123145969.
  2. 1 2 3 4 5 Fomel, Sergey; Claerbout, Jon F. (9 October 1997). "Constructing an isotropic Laplacian operator". Exploring three-dimensional implicit wavefield extrapolation with the helix transform.{{cite book}}: |website= ignored (help)
  3. Lindeberg, Tony. "Scale-Space for Discrete Signals" (PDF). Stockholm: Royal Institute of Technology. pp. 22–23.
  4. Oono, Y.; Puri, S. (1987). "Computationally efficient modeling of ordering of quenched phases" (PDF). Physical Review Letters. 58 (8): 837. Bibcode:1987PhRvL..58..836O. doi:10.1103/PhysRevLett.58.836. PMID   10035049.  Lock-green.svg
  5. Oono, Y.; Puri, S. (1988). "Study of phase-separation dynamics by use of cell dynamical systems. I. Modeling" (PDF). Physical Review A. 38 (1): 436. Bibcode:1988PhRvA..38..434O. doi:10.1103/PhysRevA.38.434. PMID   9900182.  Lock-green.svg
  6. Sevink, G. J. A. (2015). "Rigorous embedding of cell dynamics simulations in the Cahn-Hilliard-Cook framework: Imposing stability and isotropy". Physical Review E. 91 (5): 4. Bibcode:2015PhRvE..91e3309S. doi:10.1103/PhysRevE.91.053309. hdl: 1887/3194035 . PMID   26066281.  Lock-green.svg
  7. 1 2 "Rotation-invariant Laplacian for 2D grids". 21 March 2021.
  8. 1 2 3 4 "Investigation of Isotropic Laplacian Operators by Computer Simulation for A-B Diblock Copolymers" (PDF). IJCSNS International Journal of Computer Science and Network Security. 18 (5): 102–103. May 2018. Lock-green.svg
  9. 1 2 Thampi, Sumesh P.; Ansumali, Santosh; Adhikari, R.; Succi, Sauro (2013), "Isotropic discrete Laplacian operators from lattice hydrodynamics", Journal of Computational Physics, 234: 3–7, arXiv: 1202.3299 , Bibcode:2013JCoPh.234....1T, doi:10.1016/j.jcp.2012.07.037, S2CID   14633171
  10. Lynch, Robert (7 January 1992), "Fundamental Solutions of 9-point Discrete Laplacians; Derivation and Tables", Department of Computer Science Technical Reports: 2  Lock-green.svg
  11. Provatas, Nikolas; Elder, Ken (2010), Phase-Field Methods in Materials Science and Engineering (PDF), p. 219, doi:10.1002/9783527631520, ISBN   9783527631520   Lock-green.svg
  12. McGinty, Bob (2012). "Coordinate Transforms". continuummechanics.org.
  13. 1 2 3 Witkin, Andrew; Kass, Michael (1991), "Reaction-diffusion textures" (PDF), ACM SIGGRAPH Computer Graphics, 25 (4): 3–4, doi:10.1145/127719.122750 Lock-green.svg
  14. Sanderson, Allen R.; Kirby, Robert M.; Johnson, Chris R.; Yang, Lingfa (2006), "Advanced Reaction-Diffusion Models for Texture Synthesis" (PDF), Journal of Graphics Tools, 11 (3): 4, doi:10.1080/2151237X.2006.10129222, S2CID   13132043
  15. Garnier, David-Henri; Schmidt, Martin-Pierre; Rohmer, Damien (2022), "Growth of oriented orthotropic structures with reaction/Diffusion", Structural and Multidisciplinary Optimization, 65 (11): 7–8, doi:10.1007/s00158-022-03395-7, S2CID   253304840