Crank–Nicolson Method: Derivation, Implementation, and Examples

When you model how heat moves in a metal bar, how dye spreads in water, or how an option price changes over time, you often end up with an equation that evolves in time and space. To compute a good answer on a computer, you must turn that continuous change into small, careful steps. The Crank-Nicolson method is a popular way to do that. It blends two simple ideas to get a result that is both stable and accurate for many practical cases.

You may already know about “explicit” time stepping, which is simple but can blow up if you take steps that are too large, and “implicit” time stepping, which is more stable but can be more diffusive and sometimes less accurate for the same step size. The Crank–Nicolson method sits in the middle. It takes a balanced view of the present and the future state within each step. This balance is the key to its good behavior.

In this article, we explain the method in clear words, without formulas. We show how to build it, how to code it, and how to test it. We also include two tables you can use as quick guides: one to compare common time-stepping schemes and one to help you debug common issues. By the end, you will know when to pick the Crank–Nicolson method, how to implement it in one or more dimensions, and how to avoid the traps that often slow people down.

What the Crank–Nicolson Method Is (and Why It Helps)

At its core, the Crank–Nicolson method is a time-stepping rule for problems that change over time and also vary across space. Think of heat spreading along a rod, pollution moving in soil, or a smooth quantity that diffuses and may also move slightly with a flow. You divide space into points and time into small steps. You then predict the next time level from the current one.

There are three well-known ways to do this:

• Explicit step (forward in time): Use only values at the current time to predict the next time. It is fast per step because you do not need to solve a large system. But it can be unstable unless the time step is very small.

• Implicit step (backward in time): Use values at the next time to define the update. You must solve a system at each step. But it is stable for large time steps in many cases, though it may be a bit more diffusive.

• Crank–Nicolson: Take the average between the explicit and implicit updates across the same time step. This produces a second-order step in time (that is, it is often more accurate for the same time step) while still being stable for many diffusion-like problems.

The name “Crank–Nicolson” comes from John Crank and Phyllis Nicolson, who described this balanced time stepping for heat flow problems. Today, the method is used across physics, engineering, finance, and earth sciences. It is especially well-liked when you want a good blend of accuracy and stability without taking tiny time steps.

Derivation in Plain Words

We will outline the logic behind the method without symbols. The goal is to show the idea, not to prove it in a strict mathematical way.

1. Start from a Continuous Model

You have a quantity that changes in time and varies in space. For example, temperature along a rod. The rule that describes this change has two parts: how fast time moves the value, and how space differences push or spread the value.

2. Lay a Grid Over Space

Choose a set of points along the rod: left end, right end, and points in between. The distance between neighbor points is your “space step.” At each point, you will store a number that approximates the true value.

3. Lay Steps Over Time

Choose a small time step. You will update all the grid values from the current time to the next time. Repeat this many times to simulate an interval.

4. Approximate Space Changes by Central Differences

To estimate how the value curves in space at a point, you compare it with its left and right neighbors. This is a symmetric idea: it looks both ways and produces a fair approximation.

5. Approximate Time Change by an Average Slope

The central idea of Crank–Nicolson is to treat the time change across one step as the average between the “now” slope and the “next” slope.

• The “now” slope is what an explicit method would use.
• The “next” slope is what an implicit method would use.
• Taking their average gives a balanced update.

6. Turn the Rule Into a Linear System

Because the “next” slope depends on unknown next-time values, you cannot just compute the new values in one pass as with an explicit method. Instead, you collect all relations for all grid points into a set of linear equations. In one dimension, this system has a special structure: each row refers to a point and its two neighbors. This makes the matrix “tridiagonal,” which is fast to solve.

7. Solve the system at each step.

For each time step, you build a right-hand side using current values, and then solve the tridiagonal system to get the next values. You repeat until you reach your final time.

From this flow, we can sum up the idea: the Crank–Nicolson method averages the spatial action between the current and next time, which often brings second-order accuracy in time and good stability, while the central difference in space often brings second-order accuracy in space. The cost is that you must solve a linear system at each step.

Also Read: Top 10 Crypto Treasury Management Firms to Check Out in 2025

Implementation Guide You Can Follow

This section walks through a clean, code-ready plan. We keep the words simple and avoid symbols.

Define the Problem Data

• Domain and grid: Choose the size of the spatial domain. Pick the number of internal points. Compute the space step as domain length divided by the number of intervals.

• Time range and step: Choose the final time you want to reach. Pick the time step. Compute the number of steps as final time divided by time step (rounded as needed).

• Material or model coefficients: Set any constants in your model (for example, a diffusion rate). In more advanced cases, these can vary in space and time.

• Initial state: Create an array with the value at each grid point at time zero. This could be from a known function, from data, or just zeros with a spike at one location.

• Boundary conditions: Decide what happens at the edges. Common choices: fixed value at each boundary, fixed slope (no flux), or a mixed rule. We will discuss how to handle these later.

Build the Tri-Diagonal System

For a one-dimensional diffusion-type problem with central differences in space and a Crank–Nicolson average in time, the system each step links each interior point to its two neighbors. You can store the three diagonals as three arrays:

• Left diagonal (a): links to the left neighbor.
• Center diagonal (b): links to the point itself.
• Right diagonal (c): links to the right neighbor.

These arrays have fixed values when coefficients are constant. If coefficients vary in space, you compute each entry based on local values. If coefficients vary in time, you update them at each step. The system does not depend on the current solution values except through the right-hand side vector, so you can precompute the matrix if the coefficients are fixed.

Build the Right-Hand Side (RHS)

The Crank–Nicolson average mixes current and next levels. The “next” part moves to the left side of the system, while the “current” part remains on the right side. At each step:

• Start with the current solution array.
• Apply the spatial operator using the current array (the same central differences you would use in an explicit step).
• Combine it with the current array itself to form the RHS vector.
• Adjust the first and last entries in the RHS to reflect boundary conditions.

Solve the Tri-Diagonal System

Use the Thomas algorithm (a standard direct solver for tridiagonal systems). It runs in linear time with very small memory cost. The steps are:

1. Forward sweep: modify the center and RHS to eliminate the left diagonal.

2. Backward substitution: solve from the right end to the left.

Many numerical libraries offer a tridiagonal solver. In higher-level languages, you can also use a sparse direct solver. If you have many steps and a fixed matrix, it may be worth factoring the matrix once and reusing the factorization.

Apply Boundary Conditions

• Fixed value at a boundary (Dirichlet): Set the boundary point to the prescribed value at each step. In the system, you either remove that row/column and adjust RHS, or you force the center diagonal to one and the RHS to the value.

• No-flux or fixed slope at a boundary (Neumann): Replace the missing neighbor value with a mirror value or use a one-sided spatial difference that keeps the same order of accuracy. Adjust the first or last row accordingly.

• Mixed boundary rule: Combine value and slope information in the boundary row. This leads to a simple change in the first or last equation of the system.

Take care to keep the same level of accuracy at the boundary as in the interior. A sloppy boundary can spoil the whole solution.

March in Time

Repeat for each step:

1. Build the RHS from the current solution.
2. Adjust for boundary conditions.
3. Solve the tridiagonal system to get the next solution.
4. Swap “current” and “next.”

Store or plot results at selected times as needed.

Comparing Time-Stepping Schemes

Method	Time Accuracy	Space Accuracy (with central differences)	Stability behavior	Work per Step	Notes
Explicit (Forward)	First order	Second order (typical)	Often limited by step size	Very low	Simple to code; small steps often needed
Implicit (Backward)	First order	Second order (typical)	Strong for many diffusion problems	Higher (solve a system)	More diffusive feel; large steps possible
Crank–Nicolson	Second order	Second order (typical)	Good for many diffusion problems	Higher (solve a system)	Balanced update; good accuracy per step

Notes: “Order” here is a rough measure of how fast error falls when you cut the step size. Second order in time means that halving the time step tends to cut the time error by about four, all else equal. This is a rule of thumb; actual results vary with the model, the boundary rules, and the smoothness of the solution.

Examples and Use Cases

This section shows where the Crank–Nicolson method works well and how you would set each case up in practice, all without formulas.

Heat flow in a 1-D rod

Scenario: A long, thin rod with fixed temperature at both ends. The rod starts warmer in the middle. Over time, heat spreads toward the ends until the rod reaches a steady state.

Setup:

• Grid: choose a set of evenly spaced points along the rod.
• Initial state: higher in the middle, lower near the ends.
• Boundaries: left end fixed to one value, right end fixed to another (or both equal).
• Coefficients: constant diffusion rate.

With Crank–Nicolson:

Each step averages the spatial action between now and the next moment. You will see the peak in the middle flatten and spread out. With a reasonable time step, the method is stable and accurate. If you double the number of points in space and cut the time step to keep a similar balance, the solution should get sharper and closer to the expected shape at each time.

What to watch:

• Make sure the first and last rows in your system enforce the fixed values.
• If the rod is very long or needs very fine detail, use a tridiagonal solver rather than a dense solver to keep run time low.

Groundwater Diffusion in 1-D

Scenario: A simple model of water head along a straight soil column. The head evolves over time due to diffusion. One end is held at a fixed head. The other end allows no flux.

Setup:

• Grid: points from left boundary to right boundary.
• Initial state: a gentle slope in the head or a pulse.
• Left boundary: fixed head.
• Right boundary: no-flux (slope zero).
• Coefficients: can vary in space to mimic changing soil properties.

With Crank–Nicolson:

You will include the no-flux boundary by mirroring or by a one-sided space difference. The method handles spatial changes in soil properties well. If the coefficient jumps sharply, keep the grid aligned with the jump so that the property change lands on a grid point. This improves accuracy and avoids small oscillations.

What to watch:

• When coefficients vary, build the tridiagonal arrays from the local values at each interior point.
• Be careful with the boundary row when mixing fixed value on one side and no-flux on the other.

Option Pricing on a Uniform Grid

Scenario: Compute a fair value for a European option. The model can be recast as a diffusion-type problem in a transformed space coordinate and time.

Setup:

• Grid: choose a range of the transformed asset variable. Use more points near the expected strike region for better detail.
• Initial state: payoff at maturity.
• Boundaries: choose stable conditions at the far left and far right of the grid that match known behavior (for example, value goes to zero or grows linearly).
• Coefficients: can be constant or depend on the asset variable.

With Crank–Nicolson:

You step “backward” in time from maturity to the present. The balanced update helps reduce numerical damping while keeping things stable. Many classic option pricing codes use Crank–Nicolson as a solid default. If you see small ripples near the strike, try a slightly smaller time step, add a mild smoother after each step, or use a small blend toward the fully implicit method.

What to watch:

• Boundary choice can drive accuracy here. Study the expected behavior at far ends to set those rows correctly.
• If drift terms are present, the system is no longer symmetric; take care with the sign and with upwinding if needed.

2-D Heat on a Plate

You can extend the same idea to two space dimensions. The spatial operator now links each point to its four neighbors (up, down, left, right). The full matrix is no longer tri-diagonal, but it is still sparse. There are two common paths:

• Direct sparse solve: Build the sparse matrix and call a sparse solver.
• Alternating direction steps (ADI): Split each full step into two half steps, where you solve tridiagonal systems along rows in the first half and along columns in the second half. This keeps each sub-solve cheap and is easy to implement.

Crank–Nicolson with ADI is a classic pair for 2-D diffusion problems.

Also Read: Understanding Quant Hedge Funds: Strategies, Trends & AI

Practical Steps, Edge Cases, and Performance

This section gathers hands-on tips that save time and reduce error. It is based on what often goes wrong in first attempts.

Picking Time and Space Steps

• Balance matters: If you take a very large time step with very coarse space steps, you can still get a stable result, but details may be lost. Try to pick steps so that the physical pattern moves a modest fraction of a grid cell per time step.

• Refinement test: Run the same case with finer space and time steps. If the answer changes a lot, your original steps were too large. If it changes only a little, you are in a good range.

• Second-order in time: The method tends to improve quickly as you reduce the time step. If you cut the time step by half, you often see the error drop by roughly four, when the solution is smooth and the boundaries are accurate.

Coefficients that Vary in Space or Time

• Space-varying values: Compute each row of the matrix from the local coefficient. If a value changes sharply across the domain, align the grid with that change or use local averaging that respects the jump.

• Time-varying values: Rebuild the matrix or at least its diagonals at each step if needed. If rebuilding is too slow, consider a small lag (use the previous step’s matrix) and check if the error is acceptable.

Boundaries Without Pain

• Fixed value: Replace the boundary row with a simple rule that sets the boundary point to the desired value. Adjust the RHS so that interior rows do not “see” the old boundary value.

• No-flux or fixed slope: Use a one-sided difference that keeps the same accuracy. This means the first interior point and the boundary point work together to match the desired slope. You can also mirror the interior value to the outside “ghost” point.

• Mixed boundary: Combine parts of value and slope. This is common in heat transfer at a surface in contact with a fluid. The row weights both the point’s value and its slope. Keep units consistent if you work from physical data.

Avoiding Small Oscillations

Crank–Nicolson, because it averages present and future, can create small wiggles when the time step is large and the solution has sharp fronts. If that happens:

• Take a slightly smaller time step.
• Add a very small blend toward the implicit method (for example, average with a small extra weight on the future).
• Use a short, gentle smoother after each step.
• Refine the grid near the sharp front.

Speed and Memory

• Tri-diagonal wins in 1-D: The Thomas algorithm is fast and uses little memory. Precompute constant parts.

• Sparse matrices in 2-D and 3-D: Use sparse storage and solvers. Avoid dense matrices, which grow too big.

• Reuse factorization: If the matrix does not change with time, factor once and reuse factors at each step.

• Vectorization and blocking: In high-level languages, write loops so they work on whole arrays when possible. In lower-level languages, keep data in contiguous blocks to help the cache.

Testing and Validation

• Manufactured tests: Pick a simple pattern for the initial state and boundaries where you know how the solution should behave (for example, it should flatten smoothly and remain within certain bounds).

• Convergence study: Run with coarse, medium, and fine grids and time steps. Check that the results settle.

• Conservation checks: For closed systems with no flux at boundaries, the total “mass” should stay the same. Track the sum across the grid. If it drifts, you may have a boundary error or a bug in the solver.

Common Pitfalls and How to Fix Them

Symptom	Likely Cause	Fix in Practice
Solution blows up or becomes NaN	Boundary row coded wrong; wrong sign in system	Re-check boundary rows; verify diagonal signs; test a shorter step
Small ripples near sharp fronts	Time step too large for this case	Use a smaller step; add a tiny bias toward implicit; apply a mild smoother
Solution too diffusive	Space grid too coarse; first-order boundary rule	Refine grid; switch to a boundary rule with the same accuracy as interior
Slow runtime in 1-D	Dense solver used for a tri-diagonal matrix	Use a tri-diagonal (Thomas) solver; precompute coefficients
Memory spikes in 2-D or 3-D	Dense storage for large sparse matrix	Use sparse storage; consider ADI splitting in 2-D
Total “mass” not conserved (closed system)	Boundary treated like open; arithmetic error	Enforce no-flux correctly; check index ranges; compare sums at each step
Weird patterns at the edges	Ghost point rule wrong or missing	Use a correct one-sided difference or mirror value approach
Convergence stalls when refining	Mixed orders of accuracy across domain	Match boundary accuracy to interior; refine where coefficients jump
Solver does not converge (iterative method)	Poor preconditioner or bad matrix scaling	Use a direct solver if small; otherwise add a simple preconditioner and scale rows

Conclusion

The crank nicolson method gives you a stable and accurate path to simulate time-dependent processes that spread or smooth out in space. It blends present and future information within each step, which is why it performs so well for many diffusion-like problems. The price you pay is the need to solve a linear system at each step, but in one dimension this is fast, and in higher dimensions you can use sparse solvers or splitting methods.

If you are building your first solver, start with a one-dimensional case with fixed coefficients and simple boundaries. Use a tri-diagonal solver, check conservation if it applies, and do a small convergence study. Once that works, extend in small steps: variable coefficients, better boundaries, two dimensions, and so on. Keep your code clear and test each new feature before you add the next.

Finally, remember that no single method is perfect for all problems. The Crank–Nicolson method is a great default for many linear diffusion-type models, and it also works as a base that you can adapt. With the guidelines in this article, the two quick-look tables, and a careful test plan, you can implement and trust your solver in real projects—without getting lost in heavy math or complex formulas.