Solving 1D Poisson

Finite Elements

Now we solve. We have the weak form, we have a mesh, and we have the hat-function basis. This chapter assembles those three things into a linear system and solves it. The pieces fit together exactly the way the math said they would.

The element stiffness matrix

Inside one element $[x_{e}, x_{e + 1}]$ of width $h$ , only two basis functions are non-zero: the left-edge hat $ϕ_{e}$ and the right-edge hat $ϕ_{e + 1}$ . Their derivatives on this element are constants: $ϕ_{e}^{'} = - 1/ h$ and $ϕ_{e + 1}^{'} = + 1/ h$ . So the four integrals $\int_{x_{e}}^{x_{e + 1}} ϕ_{i}^{'} ϕ_{j}^{'} d x$ are products of constants times $h$ . Doing the arithmetic:

K^{e} = \frac{1}{h} (1 - 1 - 1 1) .

Symmetric, singular (its row sums are zero, which means it kills constant functions — good, those have zero derivative), and the same for every element of the same width. This little 2×2 is the entire FEM contribution from one element.

def element_stiffness(h):
    """Local 2x2 stiffness matrix for one hat-element of width h."""
    return (1.0 / h) * np.array([[ 1.0, -1.0],
                                 [-1.0,  1.0]])

The element load vector

The right-hand side $b_{i} = \int f (x) ϕ_{i} (x) d x$ breaks the same way: each element contributes a 2-vector built from the two basis functions that touch it. For a general $f$ we evaluate the integral numerically — two-point Gauss-Legendre quadrature is exact for polynomials of degree up to three, more than enough resolution for the smooth forcings we'll throw at it.

def element_load(f, x_left, x_right):
    """Local 2-vector load via 2-point Gauss-Legendre quadrature.

    Exact for any polynomial f of degree <= 3. For sin / cos forcings it's
    not exact but the error is well below the discretization error of the
    method itself, so it doesn't show up in the plot.
    """
    h = x_right - x_left
    xi = np.array([-1 / np.sqrt(3),  1 / np.sqrt(3)])  # quadrature nodes on [-1, 1]
    w  = np.array([ 1.0, 1.0])                          # weights
    x_q = x_left + (xi + 1) * 0.5 * h                   # mapped to element

    phi_left  = (1 - xi) * 0.5
    phi_right = (1 + xi) * 0.5

    f_q = np.array([f(xq) for xq in x_q])
    # Jacobian h/2 included.
    return np.array([
        np.sum(w * f_q * phi_left  * (h * 0.5)),
        np.sum(w * f_q * phi_right * (h * 0.5)),
    ])

A note on the Jacobian. The quadrature rule $\int_{- 1}^{1} g (ξ) d ξ \approx \sum w_{q} g (ξ_{q})$ sits on the reference interval $[- 1, 1]$ ; mapping it to the physical element $[x_{e}, x_{e + 1}]$ introduces a factor of $h /2$ from the change of variables. That's the trailing h * 0.5 in the code.

Assembly

Assembly is the loop that turns local element matrices into a single global system. For each element, look up which two global nodes it touches, and add its 2×2 contribution into the matching 2×2 block of the global $K$ . Same for the load vector.

def assemble(N, f):
    """Assemble global K and b, ignoring boundary conditions for now."""
    nodes, elements = mesh_1d(N)
    n = N + 1
    K = np.zeros((n, n))
    b = np.zeros(n)

    for e in range(N):
        i, j = elements[e]
        h = nodes[j] - nodes[i]

        K_e = element_stiffness(h)
        K[np.ix_([i, j], [i, j])] += K_e

        b_e = element_load(f, nodes[i], nodes[j])
        b[[i, j]] += b_e

    return K, b, nodes

The np.ix_(...) trick is just numpy's way to write "the 2×2 block at rows [i, j] and columns [i, j]." Same pattern will show up in the 2D case, just with 3×3 blocks for triangles instead of 2×2 for line segments.

Boundary conditions

For Dirichlet BCs $u (0) = u (1) = 0$ , the cleanest move is to strip the boundary nodes from the unknown set: we know $c_{0} = c_{N} = 0$ already, so they don't need to be solved for. The implementation trick is to zero those rows and columns of $K$ , put a 1 on the diagonal, and zero the matching entries of $b$ . The resulting system has the boundary unknowns trivially equal to zero and leaves the interior block untouched:

def apply_dirichlet(K, b):
    """Pin the boundary nodes (0 and N) to zero by zeroing their rows/cols
    and putting a 1 on the diagonal — the standard "row replacement" trick."""
    n = K.shape[0]
    for i in (0, n - 1):
        K[i, :] = 0.0
        K[:, i] = 0.0
        K[i, i] = 1.0
        b[i] = 0.0

For non-zero Dirichlet conditions like $u (0) = u_{0}$ , you'd set b[i] = u_0 and subtract $u_{0} \cdot K [:, 0]$ from the original RHS before zeroing the column — a "lift the boundary data" operation. The rest of the machinery is unchanged. We'll stick to homogeneous BCs for the visualizations.

Solve and check

Tie it together. The full solver is just three steps — assemble, apply BCs, call np.linalg.solve:

def solve_1d_poisson(N, f):
    """End-to-end 1D Poisson solver: -u''(x) = f(x), u(0) = u(1) = 0."""
    K, b, nodes = assemble(N, f)
    apply_dirichlet(K, b)
    c = np.linalg.solve(K, b)
    return nodes, c


# Demo: f = pi^2 * sin(pi x), so u_exact(x) = sin(pi x).
import numpy as np
nodes, c = solve_1d_poisson(N=8, f=lambda x: np.pi**2 * np.sin(np.pi * x))

print("nodes:", nodes)
print("u_h:  ", c.round(4))
print("u_ex: ", np.sin(np.pi * nodes).round(4))

nodes: [0.    0.125 0.25  0.375 0.5   0.625 0.75  0.875 1.   ]
u_h:   [0.    0.3815 0.7053 0.9213 0.9991 0.9213 0.7053 0.3815 0.   ]
u_ex:  [0.    0.3827 0.7071 0.9239 1.     0.9239 0.7071 0.3827 0.   ]

Eight elements, and the FEM solution is already within $\sim 1 0^{- 3}$ of the analytic $sin (π x)$ at every node. There is in fact a small piece of magic going on here: for 1D Poisson with hat functions, the FEM solution is exact at the nodes (assuming exact integration of the load vector). The error you see is purely from quadrature, not from the basis. Between the nodes, however, $u_{h}$ is only piecewise linear, and that's what the visualization below makes obvious.

See it

The plot below shows the forcing $f (x)$ on top, the analytic solution $u (x)$ as a smooth line, and the FEM solution $u_{h} (x)$ as a polyline through the nodes. Switch problems and refine the mesh and watch the polyline snap onto the smooth curve.

problem

elementsmax nodal error ≈ 1.67e-5

What to read off the plot:

For the smooth single-bump $sin (π x)$ , even $N = 4$ hits the peak within a couple of percent. The faceted look is the piecewise-linear nature of $u_{h}$ — between the dots, FEM doesn't know any better than to draw a straight line.
For $sin (2 π x)$ , four elements are not enough — the mesh has only one node between the two extrema, so the polyline can't capture both peaks. Doubling to $N = 8$ gives the function room to dip and rise; by $N = 16$ the FEM solution is visually indistinguishable from the analytic one.
The max nodal error reported in the controls is dominated by quadrature error (very small, hovering near $1 0^{- 4}$ ) for the trig problems. For the parabolic $f = 1$ case it drops to numerical zero — the analytic $u (x) = x (1 - x) /2$ happens to live in our hat-function space pointwise at the nodes.

A clean theoretical statement, made visible: piecewise-linear FEM converges as $O (h^{2})$ in the $L^{2}$ norm. Empirically: each doubling of $N$ quarters the max between-node error of the polyline (the error you see). Try it — eyeball the worst-case gap at $N = 4$ , then look at it again at $N = 16$ .

Where this is going

The architecture established here — element stiffness, element load, assembly loop, boundary handling, solve — is the entire FEM playbook. Chapter 4 builds a 2D triangular mesh; chapter 5 swaps the 2-node line element for a 3-node triangle and solves Poisson on the unit square. The assembly routine in 2D will look almost identical to assemble above — just with 3×3 element matrices and a different geometry. Same idea, one dimension up.