The biharmonic problem comes up in important application areas such as fluid dynamics. Methods to solve it have become more complex as computers have improved in speed and storage capacity. The methods are based on a finite-difference approximation to the given partial differential equation (PDE), leading to a system of linear algebraic equations.

This paper presents a method from Stephenson [1] that combines a second-order approximation of the solution with a fourth-order approximation of the gradient, using a nine-point stencil. Ben-Artzi et al. modify the lower order part so that the entire approximation is fourth order. Having defined the basic discrete operators, they construct a matrix representation from simple difference operators such as the symmetric tridiagonal matrix usually associated with the standard second difference operator, namely 2 on the diagonal and -1 on the off-diagonal. The modified Stephenson operator is then written as the sum of a matrix constructed from this standard matrix and a lower-rank matrix. Kronecker products are used to construct two-dimensional operators from the one-dimensional difference matrices. The resulting system is solved with a combination of the fast Fourier transform and the Sherman-Morrison formula. The cost of the solution is O(N²log N), where N is the number of grid points in one direction. Numerical results from a Fortran90 implementation of the method are given, together with error estimates.