49 P2 Finite Elements

Quadratic Finite Elements

Published

April 25, 2026

using Plots, LaTeXStrings, SparseArrays, LinearAlgebra

We now look at piecewise quadratic basis functions instead of our piecewise linear functions we used previously.

On each element, the P1 basis functions looked like “half-hat” functions. These are scaled and translated versions of these functions:

# P1 element basis functions on [-1,1]
ψ₋(ξ) = (1-ξ)/2
ψ₊(ξ) = (1+ξ)/2

plot( ψ₋, -1, 1, title="P1 basis functions", label=L"\psi_-")
plot!( ψ₊, label=L"\psi_+")

When we consider P2 basis functions, we consider 3 basis functions per element:

ϕ₋(ξ) = ξ*(ξ-1)/2
ϕ₀(ξ)  = 1 - ξ^2
ϕ₊(ξ) = ξ*(ξ+1)/2

plot( ϕ₋, -1, 1, title="P2 basis functions",label=L"\varphi_-" )
plot!( ϕ₀, label=L"\varphi_0" )
plot!( ϕ₊, label=L"\varphi_+" )

"""
    fem(n, α, β, f; a=0, b=1)

Approximate the solution to 
    -(α u′)′ + β u = f on (a,b)
    u(a) = 0, u(b) = 0
using finite element method with P1 or P2 basis functions
"""
function fem(n, α, β, f; a=0, b=1, P=1)
    if P==1

        h = (b-a) / n
        x = a .+ h .* (0:n)

        α_av = (α.(x[1:n])+α.(x[2:n+1]))/2
        β_av = (β.(x[1:n])+β.(x[2:n+1]))/2
        f_av = (f.(x[1:n])+f.(x[2:n+1]))/2

        K_d  = [ (α_av[i]+α_av[i+1])/h for i ∈ 1:n-1 ]  
        K_dl = [ -α_av[i+1]/h for i ∈ 1:n-2 ]    
        K = SymTridiagonal(K_d, K_dl)

        M_d = [ h/3 *(β_av[i] + β_av[i+1]) for i ∈ 1:n-1 ]
        M_dl = [ h/6 *β_av[i+1] for i ∈ 1:n-2 ]
        M = SymTridiagonal(M_d, M_dl)
        
        F = [ (h/2) * ( f_av[i] + f_av[i+1] ) for i = 1:n-1 ]
    
        u = (K+M) \ F
        u = [0; u; 0]
        return x, u 

    elseif P==2 

        N = 2n + 1
        h = (b-a) / n
        x = range(a, b, length=N)
        
        # Local K & M on [-1,1], scaled
        Ke_local = (1/(3h)) * [7 -8 1; -8 16 -8; 1 -8 7]
        Me_local = (h/30) * [4 2 -1; 2 16 2; -1 2 4]

        I, J, V = Int[], Int[], Float64[]
        F = zeros(N)

        # for each element
        for el ∈ 1:n
            
            idx = [2el-1, 2el, 2el+1]

            # Midpoint of element for coefficient evaluation
            x_mid = x[2el]
            α_val = α(x_mid)
            β_val = β(x_mid)

            for i ∈ 1:3, j ∈ 1:3
                push!(I, idx[i]); 
                push!(J, idx[j])
                push!(V, α_val * Ke_local[i,j] + β_val * Me_local[i,j])
            end

            # Simpson's rule
            F[idx[1]] += (h/6) * f(x[2el-1])
            F[idx[2]] += (4h/6) * f(x[2el])
            F[idx[3]] += (h/6) * f(x[2el+1])
        end

        A = sparse(I, J, V)
        u = A[2:N-1,2:N-1] \ F[2:N-1]
        u = [0; u; 0]
        return x, u 
        
    else 
        @warn "P should be 1 or 2"
    end 
end


"""
    P2_interp(y, x, u, n)

Evaluates the piecewise quadratic interpolation of (x, u) at y
"""
function P2_interp(y, x, u, n)
    h = (x[end] - x[1]) / n
    
    # find el ∋ y
    el = ceil(Int, (y - x[1]) / h)
    el = clamp(el, 1, n)
    
    idx = [2el-1, 2el, 2el+1]
    x0 = x[idx[2]]
    
    # get local coordinate ξ ∈ [-1, 1]
    ξ = 2 * (y - x0) / h
    
    return u[idx[1]]*ϕ₋(ξ) + u[idx[2]]*ϕ₀(ξ) + u[idx[3]]*ϕ₊(ξ)
end

"""
    P1_interp(y, x, u)

Evaluates the piecewise linear interpolation of (x,u) at y
"""
function P1_interp(y, x, u)
    i = searchsortedfirst(x, y)
    
    # Boundary handling
    if i <= 1
        return u[1]
    elseif i > length(x)
        return u[end]
    end
    
    # local element parameters
    h = x[i] - x[i-1]
    x_mid = (x[i] + x[i-1]) / 2
    ξ = 2 * (y - x_mid) / h
       
    return u[i-1] * ψ₋(ξ) + u[i] * ψ₊(ξ)
end

P1_interp

α(x) = 1.
β(x) = 0.
f(x) = sin(π*x)

# exact solution
u_e(x) = sin(π*x)/π^2

plot(u_e, 0, 1;
     label="exact solution", linestyle=:dash, 
     color=:black)

# P1
m = 4
x, u = fem(m, α, β, f; P=1)
plot!(x, u; 
     label="P1 FEM (n=$m)", m=5,
     color=:red)

n = 2
x, u = fem(n, α, β, f; P=2)
scatter!( x,u;
     color=:blue, markersize=5, primary=false)
plot!(y -> P2_interp(y, x, u, n), 0, 1; 
     label="P2 FEM (n=$n)", color=:blue)

# we measure the error on the grid egrid:
egrid = range(0,1,500)
Δx = step(egrid)
u_e_vals = u_e.(egrid)

nn = [2^j for j∈2:12]
err_P1 = zeros(length(nn))
err_P2 = zeros(length(nn))

for (i,n) ∈ enumerate(nn)
    x, u = fem(n, α, β, f)
    u_vals = [P1_interp(y, x, u) for y ∈ egrid]
    err_P1[i] = sqrt(Δx) * norm(u_vals - u_e_vals, 2)

    x, u = fem(n, α, β, f; P=2)
    u_vals = [P2_interp(y, x, u, n) for y ∈ egrid]
    err_P2[i] = sqrt(Δx) * norm(u_vals - u_e_vals, 2) 
end

c1 = err_P1[5] * nn[5]^2
c2 = err_P2[5] * nn[5]^3

plot( nn, err_P1;
    xaxis=:log, yaxis=:log, m=2, color=:red,
    xlabel=L"n", ylabel="error", label="P1" )
plot!( nn, err_P2;
    label="P2", m=2, color=:blue )
plot!( nn, 2c1.*nn.^(-2);  
    color=:black, linestyle=:dot, label=L"O(h^2)")
plot!( nn, 2c2.*nn.^(-3);
    color=:black, linestyle=:dash, label=L"O(h^3)")

Other Links