Galerkin P1 on the log-transformed PDE

The weak form

The PDE is the same one the finite-difference module solves — in x = ln S, forward in τ = T - t:

with a = σ²/2 and b = r - q - σ²/2. The finite element method does not discretise this directly. It multiplies by a test function φ and integrates the second derivative by parts:

where (f, g) denotes the inner product over the domain. The second derivative is gone; what remains is first derivatives of both the solution and the test function, which a piecewise-linear basis can represent exactly.

P1 hat basis

The solution is expanded in P1 (piecewise-linear, continuous) hat functions. Each basis function is 1 at its own node, 0 at every other node, and linear in between — a triangular hat. On a uniform mesh of size h there are n+1 nodes and n elements, and an interior row of the assembled system couples only the node and its two neighbours.

Writing V(x) = Σ vⱼ φⱼ(x) and taking each hat as a test function gives one equation per interior node. The system is:

M is the mass matrix, K the stiffness matrix, and B the convection matrix. All three are tridiagonal, so the θ-time- stepping — exactly the same one as the finite-difference module — builds a tridiagonal system and solves it with the same Thomas algorithm.

Element matrices

For a uniform mesh with element size h, the per-row contributions are:

The system matrix is assembled row by row as A_row = -a·K_row + b·B_row - r·M_row. Because the mesh is uniform, every interior row is the same, so the lower, diagonal, and upper vectors are filled with constants before the time loop:

const double a_low = -a*(-1/h) + b*(-0.5) - r*m_off;
const double a_dia = -a*( 2/h) + b*( 0.0) - r*m_dia;
const double a_upp = -a*(-1/h) + b*( 0.5) - r*m_off;

for (size_t k = 0; k < m; ++k) {
    lower[k] = m_off - θ·dt·a_low;
    diag[k]  = m_dia - θ·dt·a_dia;
    upper[k] = m_off - θ·dt·a_upp;
}

The θ·dt term is what puts the time-stepping into the matrix. For θ = 0 (explicit) the matrix is just the mass matrix row; for θ = 1 (implicit) it is the full M - dt·A row; for θ = ½ it is the Crank–Nicolson blend.

Mass matrix: consistent vs lumped

The finite-difference module uses the identity mass matrix implicitly — each node is its own degree of freedom with no coupling to its neighbours in the time term. The finite-element module offers two choices:

The consistent mass matrix is the correct one for a Galerkin method — it is the L² projection of the solution onto the hat basis. The lumped mass is a common approximation that is cheaper (diagonal, no coupling) and, for this problem, happens to reproduce the finite-difference scheme to within 2 × 10⁻⁶. The implementation offers both so that the equivalence can be demonstrated rather than asserted.

The finite-difference equivalence

Lumping the mass matrix collapses the off-diagonal coupling in the time term to zero and sets the diagonal to h. The system that results is:

The -a/h K + b/h B - r I row is, element by element, the same as the finite-difference operator row α, β, γ — the stiffness gives -a/h · (-1, 2, -1), the convection gives b · (-½, 0, ½), and the reaction gives -r on the diagonal. So the lumped FE scheme is the finite-difference scheme with a particular (and correct) identification of the step size.

Time-stepping

The time-stepping is identical in structure to the finite-difference module, because the matrices are tridiagonal in both cases and the Thomas solver is shared (finite_element.hpp includes finite_difference.hpp for solve_tridiagonal, FdScheme, and theta_of).

Each step builds the right-hand side as (M + (1-θ)·dt·A) Vⁿ with the boundary terms moved across, solves the tridiagonal system, and applies American projection if requested:

for (size_t k = 0; k < m; ++k) {
    const size_t i = k + 1;
    rhs[k] = (m_off*v[i-1] + m_dia*v[i] + m_off*v[i+1])
           + w*(a_low*v[i-1] + a_dia*v[i] + a_upp*v[i+1]);
}
rhs[0] -= (m_off - θ·dt·a_low) * v_lo;
rhs[m-1] -= (m_off - θ·dt·a_upp) * v_hi;

solve_tridiagonal(lower, diag, upper, rhs, scratch);

The boundary conditions are the same asymptotics as the finite-difference module: deep out-of-the-money worthless, deep in-the-money call = discounted forward, put = discounted strike less discount spot, with the early-exercise maximum for an American contract.

Boundaries and American exercise

The two boundary nodes are imposed before each step from the same asymptotics as the finite-difference module. For a European contract they are the discounted forward (call) or the discounted strike-minus-spot (put). For an American contract the boundary value is the maximum of that asymptotic value and the payoff.

American exercise is handled the same way as in the finite-difference module — explicit projection onto the payoff after each step:

if (opt.exercise == Exercise::American)
    for (size_t i = 0; i <= n; ++i)
        v[i] = std::fmax(v[i], opt.payoff(exp(x[i])));

The same limitation applies: this is explicit projection, first order in time near the free boundary even under Crank–Nicolson, and stated for the same reason it is stated in the finite-difference module.

Source