Skip to content
Get started

3. Weak forms and finite elements

How can straight segments solve a curved problem?

Section titled “How can straight segments solve a curved problem?”

A finite-element solution may be piecewise linear even when its differential equation contains a second derivative. Within each segment that second derivative is zero. How can it represent a nonzero source?

The answer is to ask the equation to hold in an integrated sense. We will derive that statement, construct one element, and run a mesh comparison. For physical context, uu can represent a scaled temperature difference in steady conduction.

Use a unit interval and dimensionless coordinate ξ=x/L\xi=x/L, with L=1mL=1\,\mathrm m. The scaled problem is

u(ξ)=f(ξ)=π2sin(πξ),u(0)=u(1)=0.-u''(\xi)=f(\xi)=\pi^2\sin(\pi\xi),\qquad u(0)=u(1)=0.

Here primes denote differentiation with respect to ξ\xi. Direct substitution gives u(ξ)=sin(πξ)u(\xi)=\sin(\pi\xi), so the field is symmetric, positive inside the interval, and has maximum 1 at the center. Its endpoint derivatives are u(0)=πu'(0)=\pi and u(1)=πu'(1)=-\pi.

The .eqi source below uses dimensional xx: its wave number is π/L\pi/L and source scale is π2/L2\pi^2/L^2. Keeping these scales explicit makes the sine’s argument dimensionless.

Multiply by a test function vv vanishing at both endpoints and integrate:

01(u)vdξ=[uv]01+01uvdξ=01fvdξ.\int_0^1(-u'')v\,d\xi =-[u'v]_0^1+\int_0^1u'v'\,d\xi =\int_0^1fv\,d\xi.

The boundary term vanishes because vv vanishes there. The weak problem is therefore: find a function satisfying the endpoint values such that

a(u,v)=(v)for all admissible v,a(u,v)=01uvdξ,(v)=01fvdξ.a(u,v)=\ell(v)\quad\text{for all admissible }v, \qquad a(u,v)=\int_0^1u'v'\,d\xi,\quad \ell(v)=\int_0^1fv\,d\xi.

Only first derivatives now appear. Continuous piecewise linear functions have square-integrable first derivatives even though their slopes jump at nodes. That makes them admissible trial functions. A formal account of the spaces and boundary traces is in Schöberl’s weak formulation of Poisson’s equation.

On an element [ξi,ξi+1][\xi_i,\xi_{i+1}] of length hh, take

N1=(ξi+1ξ)/h,N2=(ξξi)/h,uh=N1Ui+N2Ui+1.N_1=(\xi_{i+1}-\xi)/h,\qquad N_2=(\xi-\xi_i)/h, \qquad u_h=N_1U_i+N_2U_{i+1}.

Their derivatives are 1/h-1/h and 1/h1/h. Integrating each product of derivatives gives the element matrix

K(e)=1h(1111),Fa(e)=ξiξi+1fNadξ.K^{(e)}=\frac1h \begin{pmatrix}1&-1\\-1&1\end{pmatrix},\qquad F^{(e)}_a=\int_{\xi_i}^{\xi_{i+1}}fN_a\,d\xi.

Notice two structural facts before looking at any output. The matrix is symmetric, and each row sums to zero. A constant local field therefore creates no diffusive flux. Essential endpoint values remove the global constant nullspace in this problem.

Two neighboring elements contribute to the equation at their common node. On a uniform interior mesh their stiffness contributions produce

Ui1+2UiUi+1h=Fi.\frac{-U_{i-1}+2U_i-U_{i+1}}h=F_i.

For constant ff, Fi=fhF_i=fh, giving the familiar centered second difference after division by hh. For a variable source, keep the integrated load: replacing it by hf(ξi)hf(\xi_i) is an additional approximation.

As a hand calculation, use f=2f=2 and two elements of length 1/21/2. The sole interior equation is 4U1=14U_1=1, so U1=1/4U_1=1/4. The exact solution u=ξ(1ξ)u=\xi(1-\xi) also has value 1/41/4 there. Between the nodes the numerical field is still linear, while the exact field is curved. Exact nodal values do not imply an exact field.

poisson.eqi
model manufactured_poisson() {
domain interval = box(0, 1);
domain lower_end = boundary(interval, axis = 0, side = lower);
domain upper_end = boundary(interval, axis = 0, side = upper);
variable potential: 1 on interval;
parameter wave_number: 1 / m = 3.141592653589793;
parameter source_scale: 1 / m ^ 2 = 9.869604401089358;
relation balance on interval {
-div(grad(potential)) - source_scale * math.sin(wave_number * coordinate(0)) = 0;
}
relation lower_value on lower_end {
trace(potential) = 0;
}
relation upper_value on upper_end {
trace(potential) = 0;
}
}

From the source checkout created in Get started:

Terminal windowbash
cargo run -p eqiora-numerics --example poisson_convergence

Read fem_l2 first. It measures the error in the continuous piecewise linear field, integrated over the interval. The comparison program runs meshes with 8, 16, 32, 64, and 128 cells. Predict a decreasing error before reading it; the refinement lesson derives the order calculation.

The complete run program reads this equation and compares finite elements with finite volumes. The choice of approximation belongs to the computation, while the source states the mathematical problem.

For Python applications, the corresponding choice is eqiora.fem.Q1(). In one dimension it gives the linear construction just derived. On Cartesian rectangles, Q1 uses products of one-dimensional linear basis functions; the resulting fields are bilinear within each cell. This product structure also explains why mesh geometry and quadrature matter.

We imposed endpoint values by restricting the trial field and taking zero test values there. These are essential conditions. If a test function is allowed to be nonzero at an endpoint, the integration-by-parts boundary term remains and carries derivative or flux data: this is the origin of a natural condition. Always derive its sign using the outward normal. Do not silently drop that term on a boundary where the field value is free.

  1. Assemble the matrix for three equal elements with zero endpoint values. Which entries receive contributions from two elements?
  2. For the two-element constant-source example, evaluate the numerical and exact fields at ξ=1/4\xi=1/4. Explain their difference despite the exact midpoint value.
  3. Show that a(v,v)>0a(v,v)>0 for a nonzero piecewise linear function with zero endpoint values. Why does removing both endpoint constraints change this?
  4. Derive the boundary term when the right endpoint derivative is prescribed. State separately the outward diffusive flux qn=unq_n=-u'n.
  5. In two dimensions, write the four products of linear basis functions on a unit square. Check that they sum to one.

Previous: Time integration · Book map · Next: Finite-volume balance