# Basic approach for numerical solution of PDE

I'm looking for some guidance on how to write a program to numerically solve a PDE. As an example for comparison in 1D:

$$\frac{d^2u}{dx^2} = f\;\;\;\;u(0) = 0\;\;\;\;u(1) = 0$$

We could try 6 equally spaced nodes between 0 and 1, interpolate within each element with a quadratic, and ensure continuity of $u$ and $du/dx$ at the boundaries, and to solve within each element we can integrate the differential equation over each element. I believe this constant-weight integration makes this FVM. Is that correct?

In any case, there were 3 DOFs for each of the 5 elements, 1 integral equation per element, 2 compatibility conditions for each of the 4 internal boundaries, and 2 boundary conditions given as part of the problem. So that's 15 equations in 15 variables which can then be solved. When trying to make the equivalent 2D method, I cannot work out how to generate the right number of compatibility equations. For example:

$$\nabla^2{u} = f\;\;\;\;u(0,y) = u(1,y) = u(x,0) = u(x,1) = 0$$

Imagine we mesh with a grid of 5 by 5 squares with each square divided into two triangles. Each triangular element can be represented with a quadratic in $x$ and $y$ (6 DOFs). The DE is integrated over for each element. Continuity of $u$ and $\nabla{u}$ could be ensured at the boundaries, though already this might be over/under-constrained.

There are 50 elements, 65 internal boundaries and 20 external boundaries. There are 300 DOFs. Compatibility of $u$ and $\nabla{u}$ on a single internal boundary yields 5 equations, but either the sets of 5 equations per internal boundary are not independent of each other or the system is over-constrained because already there are too many equations.

As far as I understand, representing the DOFs with the value of $u$ at the corners and edge midpoints of the triangles is a way of intrinsically satisfying the continuity of $u$ condition and reducing the number of DOFs, but it does not help intrinsically satisfy the other compatibility conditions. But I mention this in case proper representation of the DOFs is key to helping understand the problem.

If I am thinking about this in roughly the right way, how do I determine the compatibility conditions and generate the right number of equations to solve for the given number of DOFs? If I'm not thinking about the problem in the right way but you can help me funnel this into a more specific, generally-relevant question(s), please do.

• You have essentially discovered that continuity of $\nabla u_h$, the gradient of the approximation, cannot be enforced with $P_2$ Lagrange elements. If you need $\nabla u_h$ to be continuous, you will have to use a more advanced finite element. Please note that you do not have to enforce continuity of $\nabla u_h$ for a second-order method. In fact, you could even construct second-order methods where $u_h$ itself is discontinuous. – cthl Jul 26 '18 at 5:18
• Thanks, that's useful. I now understand the following. Piecewise polynomials means a certain number of DOFs per element. Continuity of the value or some derivative means a certain number of equations per boundary. In 1D, vertexes - edges = 1, which means there the number of boundary equations can be made to match the number of element DOFs. It works. But in 3D for example, the relationship is vertexes - edges + faces - cells = 1, so the number of boundary equations cannot be made to match the number of element DOFs in general. The naive generalisation does not work. – user2357 Aug 13 '18 at 14:46
• You haven't told us why you want continuity of $\nabla u_h$; at least for the Poisson equation, that is not required by standard theory. – Christian Clason Aug 14 '18 at 7:00
• I'm just trying to generalise the 1D method described, for which an $n^{th}$ order equation can approximated by a piecewise $n^{th}$ order polynomial with $C_{n-1}$ continuity which satisfies the equation when integrated over each element, leaving $n$ boundary conditions to fully determine the solution. I was not specifically trying to relate this to FEM or another standard theory (but perhaps this method does generalise on structured grids to a kind of FDM where there is a strict relationship between elements/cells and boundaries/facets). – user2357 Aug 14 '18 at 11:03