Thomas Algorithm (Tridiagonal Solver) Simulator

Parameters

Interior grid points n

Number of unknowns in the tridiagonal system (resolution)

Left boundary value u(0)

Dirichlet boundary condition (value of u at x=0)

Right boundary value u(1)

Dirichlet boundary condition (value of u at x=1)

Source magnitude f

Strength of the right-hand side of −u''=f

Source distribution

Spatial shape of the right-hand side f(x)

Results

—

Grid points n

—

Operation count (≈8n)

—

Central u value

—

Maximum u value

—

Max error (vs analytic, uniform)

—

Complexity class

—

Algorithm visualisation — forward sweep & back substitution

The solution u(x) is drawn while a highlight sweeps left→right (forward elimination) then right→left (back substitution). The small matrix at top-right shows only the three highlighted diagonals.

Numerical solution vs analytic solution

Operation count — Thomas O(n) vs Gaussian elimination O(n³)

Theory & Key Formulas

$$a_i\,u_{i-1} + b_i\,u_i + c_i\,u_{i+1} = d_i$$

Row i of the tridiagonal system. a: sub-diagonal, b: main diagonal, c: super-diagonal, d: right-hand side. For the central-difference discretisation a=c=−1, b=2.

$$c'_i=\frac{c_i}{b_i-a_i c'_{i-1}},\qquad d'_i=\frac{d_i-a_i d'_{i-1}}{b_i-a_i c'_{i-1}}$$

The forward sweep. Start at i=1 with c'_1=c_1/b_1, d'_1=d_1/b_1, then advance through i=2..n.

$$u_n=d'_n,\qquad u_i=d'_i-c'_i\,u_{i+1}$$

Back substitution, solving in reverse order i=n−1..1. The total cost is O(n), against O(n³) for general Gaussian elimination — dramatically faster.

What is the Thomas Algorithm?

🙋

I've never heard of the "Thomas algorithm". I get that it solves systems of equations, but how is it different from ordinary Gaussian elimination?

🎓

Roughly speaking, it's the "express version" of Gaussian elimination. The Thomas algorithm targets tridiagonal matrices — matrices where numbers appear only on the main diagonal and the diagonals just above and just below it, with zeros everywhere else. General Gaussian elimination solves any matrix, but its cost blows up like n³. With a tridiagonal matrix you never have to touch the zero parts, so it solves in just O(n), roughly 8n operations.

🙋

Only three diagonals... does such a convenient matrix actually show up in practice?

🎓

It shows up surprisingly often. The very problem this tool solves is one. When you discretise the 1D differential equation −u''=f with central differences — the standard approach — the equation at each node involves only three unknowns: "itself", "left neighbour" and "right neighbour". Write that as a matrix and it is automatically tridiagonal. When you change the grid count n on the left, the tool solves that n×n tridiagonal system with the Thomas algorithm each time.

🙋

I see! So how exactly is it solved? In the animation at top-right an arrow moves from left to right.

🎓

Two stages. First the "forward sweep" — from the leftmost equation onward, you eliminate the sub-diagonal term. At each row you compute c'_i = c_i/(b_i − a_i·c'_{i-1}) while moving right. That is the left→right sweep in the animation. By the last row, that row has only one unknown left and u_n is determined immediately. From there comes "back substitution", going right→left: using u_i = d'_i − c'_i·u_{i+1}, you use the value just found on the right and solve one node at a time toward the left. That is the right→left sweep.

🙋

How accurate is the answer? The "max error" is showing something tiny like 1e-13.

🎓

Good question. For the default uniform source this problem has an exact analytic solution, u(x) = 50·x·(1−x). The Thomas algorithm is not an approximation — it solves the discrete system exactly, so apart from rounding the error is essentially zero, only about 1e-13, the floating-point limit. The central u should be exactly 12.5 as theory predicts. But if you switch the source to "centred" or "linear", the analytic formula changes, so the error field shows "—". The uniform case is the one that confirms the Thomas algorithm itself is exact.

🙋

If it's that fast and accurate, I almost feel we could just use the Thomas algorithm for everything.

🎓

I understand the feeling, but you can only use it "when the matrix is tridiagonal". On top of that, the Thomas algorithm does no pivoting, so unless the matrix is diagonally dominant (the diagonal entry at least the sum of the off-diagonals), the denominator can approach zero mid-sweep and errors run wild. Luckily this tool's matrix has diagonal 2 and off-diagonals −1, so it is diagonally dominant and safe. In practice, once you have "tridiagonal plus diagonally dominant", reaching for the Thomas algorithm is the standard move.

Frequently Asked Questions

The Thomas algorithm is a dedicated method for solving a tridiagonal linear system A·u = d — one whose only non-zero entries lie on the main diagonal and the diagonals just above and below it. It is general Gaussian elimination specialised to this banded structure, and consists of two stages: a forward sweep and a back substitution. General Gaussian elimination needs about (2/3)n³ operations, whereas the Thomas algorithm needs only about 8n, i.e. O(n). Tridiagonal systems appear everywhere in engineering — 1D boundary-value problems, cubic spline interpolation, implicit time-stepping of PDEs — which makes it one of the most important basic algorithms in numerical computing.

In general Gaussian elimination each elimination step updates all rows and columns below it, so the cost grows like n³. In a tridiagonal matrix every row has at most three non-zero entries, and eliminating one row affects only its single neighbour. The forward sweep therefore costs a constant number of operations per row, giving a total proportional to n. Back substitution is also one subtraction and one multiplication per row, so it too is O(n). Altogether the work is roughly 8n floating-point operations, and the memory is just the three diagonals — O(n).

The Thomas algorithm does no pivoting, so it is numerically stable when the matrix is diagonally dominant (in each row the diagonal entry is at least as large in magnitude as the sum of the off-diagonals) or symmetric positive definite. The 1D Poisson matrix solved by this tool has diagonal 2 and off-diagonals −1, so it is diagonally dominant and solves stably. For a general tridiagonal matrix that is not diagonally dominant, the denominator m = b_i − a_i·c'_{i-1} can approach zero during the forward sweep and amplify errors; then a banded solver with partial pivoting is required.

The most typical case is the 1D boundary-value problem: discretising a second derivative with central differences, as in the −u''=f equation this tool solves, always produces a tridiagonal matrix. Others include determining cubic-spline coefficients, every step of an implicit time integration of heat or diffusion equations (such as Crank-Nicolson), the splitting in the ADI method for 2D problems, and tridiagonal structural systems. Whenever a tridiagonal system appears, the standard move is to check whether the Thomas algorithm can solve it in O(n).

Real-World Applications

Implicit time integration of PDEs: When you advance a heat or diffusion equation in time, using an implicit scheme (such as Crank-Nicolson) for stability means solving one tridiagonal system at every time step. A 10,000-step simulation means solving 10,000 times with the Thomas algorithm. If you used O(n³) Gaussian elimination here the run time would balloon by thousands of times, so the Thomas algorithm is what makes implicit schemes practical.

Cubic spline interpolation: Cubic splines, used in CAD, graphics and data visualisation to draw smooth curves, solve a system of equations to determine each segment's coefficients. The equations that arise from the spline's "continuity of the second derivative" condition are exactly tridiagonal. Even a smooth curve through hundreds of points has its coefficients fixed in an instant by the Thomas algorithm.

1D structural and heat-transfer analysis: Beam deflection, fin temperature distribution, conduction in the ground — 1D boundary-value problems are the Thomas algorithm's home turf. The −u''=f that this tool solves is exactly this archetype, the 1D version of the Poisson and Laplace equations. In CAE teaching and verification, this simple problem is used first to confirm that the discretisation and the solver are correct before moving on to 2D and 3D.

As a building block of large sparse solvers: Even for 2D and 3D problems, splitting the problem one direction at a time, as in the ADI (alternating-direction implicit) method, reduces to solving large numbers of tridiagonal systems internally. The idea of the Thomas algorithm is also used in block-tridiagonal solvers and in the preconditioning stage of preconditioned iterative methods. It works not only as a standalone solver but as the beating heart of large-scale solvers.

Common Misconceptions and Pitfalls

A common one first: "the Thomas algorithm safely solves any tridiagonal matrix". The Thomas algorithm does no pivoting (no row swapping) at all, so if the forward-sweep denominator m = b_i − a_i·c'_{i-1} approaches zero, errors are amplified at once. This is guaranteed not to happen only when the matrix is diagonally dominant or symmetric positive definite. This tool's 1D Poisson matrix (diagonal 2, off-diagonals −1) is diagonally dominant and therefore safe, but feeding in a general tridiagonal matrix as-is is risky. When it is not diagonally dominant, use a banded solver with partial pivoting (such as LAPACK's dgtsv).

Next, assuming that "because it is O(n), the error is also zero". The Thomas algorithm solves the discretised linear system exactly (apart from rounding), but that means the "answer to the discrete equation" is exact — which is not the same as the "answer to the original differential equation". The central-difference discretisation has an O(h²) truncation error, and increasing the grid count n brings you closer to the true solution. The reason the error is essentially zero in this tool's uniform-source case is a special circumstance: for this particular problem the discrete and analytic solutions happen to coincide. In general you must distinguish the "algorithm's rounding error" from the "discretisation's truncation error".

Finally, "if it is tridiagonal, the Thomas algorithm is always fastest" is not necessarily true. As a direct method the Thomas algorithm is optimal, but if the coefficient matrix is reused for multiple right-hand sides, computing an LU factorisation once (also O(n)) lets you repeat just the forward and back substitution. Also, a cyclic tridiagonal matrix — with entries in the "corners" from periodic boundary conditions — cannot be solved by the pure Thomas algorithm and needs a correction such as the Sherman-Morrison formula. Look closely at the problem's structure and choose between the Thomas algorithm, the block Thomas algorithm and the cyclic variant — that is the correct attitude in practice.

What is the Thomas Algorithm?

Frequently Asked Questions

Real-World Applications

Common Misconceptions and Pitfalls

How to Use

Worked Example

Practical Notes