Numerical Differential Equations

The Central Idea

Consider an initial-value problem

A numerical method constructs approximations

at grid points

Ideally,

Why Numerical ODE Methods Are Needed

Many differential equations cannot be solved conveniently in closed form. Numerical methods are used when:

• the equation is nonlinear;
• the system has many interacting variables;
• coefficients vary with time or state;
• forcing terms are complicated;
• the analytic solution exists but is impractical to evaluate;
• the differential equation is part of a larger simulation.

From Derivative to Step

The derivative tells us the local slope:

Over a small interval of length \(h\),

Substituting the differential equation gives the simplest numerical update.

Euler's Method

Euler's method follows the tangent line at the beginning of each step.

Geometric Interpretation of Euler's Method

At \((t_n,y_n)\), the differential equation gives the slope

Euler's method moves a horizontal distance \(h\) along the tangent-line approximation.

The new point becomes the starting point for the next tangent step.

Example: One Euler Step

Suppose

With \(h=0.1\),

The exact value is

Euler as a Taylor Method

Taylor's theorem gives

Euler's method keeps only the first two terms:

Local Truncation Error

The local truncation error measures the error introduced in one step assuming the step begins from the exact solution value.

For Euler's method,

Global Error

The global error at step \(n\) is

Because many local errors accumulate, Euler's global error is typically

Euler's method is therefore first-order convergent.

Local and Global Error Are Different

A method may have a local error of one order higher than its global error.

For Euler:

The difference comes from accumulation over roughly \(1/h\) steps on a fixed time interval.

Smaller Steps

Reducing \(h\) usually improves accuracy because each step follows the solution over a shorter interval.

But smaller steps also require more function evaluations and may eventually increase floating-point accumulation.

Explicit Methods

Euler's method is explicit because

is computed directly from already known quantities.

Implicit Methods

An implicit method places the unknown future value inside the equation.

The backward Euler method is

Each step may require solving a nonlinear equation.

Why Use an Implicit Method?

Implicit methods cost more per step, but they can have much stronger stability properties.

This is especially important for stiff differential equations.

The Improved Euler Idea

Euler uses only the slope at the start of a step.

A better approximation uses additional slope information to represent the average behavior across the interval.

Heun's Method

First predict with Euler:

Then average the starting and predicted slopes:

Heun's Method Is Second Order

Under suitable smoothness assumptions,

and

The Explicit Midpoint Method

Estimate the midpoint state:

Then use the midpoint slope:

Runge-Kutta Methods

Runge-Kutta methods combine several slope evaluations within one step.

Their goal is to obtain higher-order accuracy without explicitly computing higher derivatives of \(y\).

Classical Fourth-Order Runge-Kutta

Define

The RK4 Update

Classical RK4 has global error

Why RK4 Is Popular

• high accuracy for smooth nonstiff problems;
• no analytic derivatives beyond \(f\) are required;
• simple fixed-step implementation;
• good balance of cost and accuracy for many introductory problems.

Cost per Step

Euler requires one evaluation of \(f\) per step.

Classical RK4 requires four.

A higher-cost step can still be more efficient overall if it permits a much larger step size for the required accuracy.

Order of a One-Step Method

If the global error behaves like

the method has order \(p\).

Halving \(h\) should reduce the leading error by approximately

Convergence Study

Compute a solution using

Compare the approximations at common times.

The pattern of differences can reveal whether the method is operating in its expected convergence regime.

Observed Order

If the exact solution is known for testing,

This same convergence-rate idea appears throughout numerical analysis.

Consistency

A method is consistent when its discrete step approaches the differential equation as

Informally, the numerical update must become a faithful local representation of the ODE at sufficiently small step size.

Stability

Stability describes how perturbations and numerical errors behave as the method advances.

A convergent formula must not allow harmless local errors to grow uncontrollably under ordinary conditions.

The Linear Test Equation

A standard stability model is

The exact solution is

If \(\operatorname{Re}(\lambda)<0\), the true solution decays.

Euler Applied to the Test Equation

Euler gives

Therefore,

Euler Stability Condition

For a decaying exact solution, the discrete sequence decays only if

Thus a mathematically stable differential equation can become numerically unstable if the step size is too large.

Absolute Stability Region

For a one-step method applied to

the update often has the form

The absolute stability region is

Stability Function

For forward Euler,

For backward Euler,

These different stability functions explain their very different behavior on stiff decay problems.

Backward Euler Stability

Backward Euler applied to \(y'=\lambda y\) gives

Hence

For \(\operatorname{Re}(\lambda)<0\), backward Euler remains stable over a very broad range of step sizes.

A-Stability

A method is called A-stable if its stability region contains the entire left half of the complex plane:

Backward Euler is A-stable.

Stiff Differential Equations

A problem is called stiff when stability requirements force an explicit method to use a much smaller step than accuracy alone would require.

Stiffness often appears when the solution contains multiple time scales, some of which decay very rapidly.

A Simple Stiffness Picture

Suppose one component decays like

while another changes on a time scale near

The fast mode may become physically negligible quickly, yet an explicit method may still need tiny steps to remain stable.

Implicit Methods and Stiffness

Implicit methods can remain stable for much larger steps on stiff problems.

The tradeoff is that each step may require solving

Root-finding and numerical linear algebra therefore reappear inside ODE solvers.

Newton's Method Inside an Implicit Solver

Backward Euler requires solving

Define

Newton iteration can then solve \(G(z)=0\).

The Trapezoidal ODE Method

Integrating

over one step gives

Applying the trapezoidal rule gives

Integration View of ODE Solvers

The exact solution satisfies

Many ODE methods can be understood as numerical integration formulas applied to this identity.

Systems of Differential Equations

Let

Euler becomes

Runge-Kutta methods extend in the same vector form.

Second-Order Equations as First-Order Systems

Consider

Introduce

Then

This converts the second-order equation into a first-order system.

Example: Harmonic Oscillator

The equation

becomes

Standard system solvers can then advance both position and velocity.

Multistep Methods

Runge-Kutta methods use several slope evaluations within one step.

Multistep methods reuse information from previous time levels.

This can reduce the number of new function evaluations per step.

Two-Step Adams-Bashforth Method

where

This is an explicit second-order multistep method.

Adams-Moulton Methods

Adams-Moulton formulas include the future slope

making them implicit.

They are often paired with explicit predictors in predictor-corrector schemes.

Predictor-Corrector Strategy

1. Use an explicit method to predict \(y_{n+1}\).
2. Evaluate the future slope using the prediction.
3. Use an implicit-style corrector formula.
4. Repeat the correction if needed.

Backward Differentiation Formulas

Backward differentiation formulas, or BDF methods, approximate derivatives using several past solution values and evaluate the ODE implicitly at the new time.

They are especially important in stiff ODE solvers.

Starting a Multistep Method

A \(k\)-step method needs several previous values before it can begin.

A one-step method such as Runge-Kutta is often used to generate the required starting values.

Variable Step Size

A fixed \(h\) is simple but may be inefficient.

A variable-step solver uses

and changes \(h_n\) according to estimated local error and solution behavior.

Embedded Runge-Kutta Methods

An embedded pair computes two approximations of different orders using largely the same slope evaluations.

Their difference estimates local error.

Adaptive Step-Size Control

If the estimated error is too large, reduce the next step.

If it is far below tolerance, increase the next step.

A common scaling idea is

Safety Factors

Practical adaptive solvers multiply the theoretical step update by a safety factor smaller than \(1\).

This reduces the chance that the next step barely misses the requested tolerance.

Absolute and Relative Error Control

A componentwise scaled error test may use

This prevents large and small solution components from being judged by the same absolute scale.

Step Rejection

Adaptive solvers do not have to accept every attempted step.

If the estimated error exceeds tolerance:

1. reject the trial step;
2. reduce \(h\);
3. recompute the step.

Dense Output

Sometimes values are needed between accepted time steps.

A solver can construct a local interpolant over each step:

This is called dense output.

Event Detection

Many simulations must stop or react when an event function satisfies

Examples include:

• impact with a surface;
• crossing a threshold;
• switching a controller;
• reaching a target state.

Event location combines ODE stepping with root-finding.

Conservation Properties

Some differential equations preserve quantities such as energy, mass, momentum, or probability.

A generic high-order method may approximate the solution accurately while slowly drifting in an invariant.

Specialized geometric integrators can preserve important structure better over long simulations.

Long-Time Integration

Small local errors can accumulate over a long simulation.

For oscillatory or conservative systems, phase error and structural drift may matter more than pointwise local accuracy.

Stability Versus Accuracy

A step can be accurate in a local Taylor sense yet unstable over repeated use.

Conversely, a stable method can be too inaccurate if the step size is overly large.

Reliable solvers must satisfy both requirements.

Roundoff Accumulation

Smaller steps increase the number of arithmetic operations.

Eventually, further reduction in \(h\) can stop improving the solution because floating-point error accumulates.

This creates another accuracy-versus-cost tradeoff.

Problem Conditioning

If nearby initial conditions

produce rapidly diverging exact solutions, the IVP itself is sensitive.

No numerical method can remove instability that belongs to the mathematical problem.

Numerical Stability Is Different from Physical Instability

A physical model may genuinely contain growing modes.

A numerical method should reproduce that growth accurately without introducing additional artificial instability.

A Reliable Euler Workflow

1. Write the IVP in first-order form.
2. Choose a starting step size.
3. Evaluate \(f(t_n,y_n)\).
4. Update with \(y_{n+1}=y_n+hf(t_n,y_n)\).
5. Repeat on the desired interval.
6. Rerun with \(h/2\).
7. Compare results to assess first-order convergence.

A Reliable Runge-Kutta Workflow

1. Confirm the problem is appropriate for an explicit solver.
2. Choose a step size or embedded adaptive pair.
3. Compute the required internal stage slopes.
4. Form the weighted step update.
5. Estimate local error when an embedded pair is available.
6. Accept or reject the step.
7. Adjust the next step size.
8. Verify convergence on representative problems.

A Reliable Stiff-Solver Workflow

1. Look for multiple time scales or explicit stability restrictions.
2. Select an implicit or stiffness-aware method.
3. Form the nonlinear step equation.
4. Solve each step with Newton or a related method.
5. Use Jacobian information or approximations when needed.
6. Control local error adaptively.
7. Monitor nonlinear-solver convergence and rejected steps.

A Reliable Verification Workflow

1. Test against an ODE with a known exact solution.
2. Compute with several step sizes.
3. Measure global error at common times.
4. Estimate observed convergence order.
5. Check invariants when the model has them.
6. Compare independent solvers when practical.
7. Distinguish modeling error from discretization error.

Common Error: Using Euler with a Large Step

Euler's method is simple but has only first-order global accuracy and a limited stability region.

A large step can produce both large truncation error and artificial instability.

Common Error: Assuming a Higher-Order Method Is Always Better

Higher order improves asymptotic accuracy, but method choice also depends on:

• stability;
• stiffness;
• function-evaluation cost;
• memory;
• event handling;
• required structure preservation.

Common Error: Ignoring Stability Restrictions

A step size may be small enough for accuracy yet too large for stability.

The test equation

provides a standard way to analyze this behavior.

Common Error: Treating Stiffness as Ordinary Inaccuracy

Repeatedly shrinking the step of an explicit solver may work, but it can be extremely inefficient.

A stiffness-aware implicit method may be the correct remedy.

Common Error: Using Only One Step Size

One numerical solution provides little evidence of convergence.

Repeat the calculation on finer grids and examine the error or solution differences.

Common Error: Ignoring the Nonlinear Solve in Implicit Methods

An implicit formula is not complete until the step equation has been solved to sufficient accuracy.

Poor nonlinear convergence can dominate the nominal discretization error.

Common Error: Overly Tight Tolerances

Extremely small requested tolerances can force excessive steps and function evaluations without producing meaningful additional accuracy.

Tolerances should reflect floating-point limits, model accuracy, and application requirements.

Common Error: Ignoring Scale Across Components

In systems, one variable may be near \(10^{-6}\) while another is near \(10^4\).

Combined absolute and relative tolerance scaling prevents one component from dominating the error test unfairly.

Example: Two Euler Steps

Solve approximately

using \(h=0.1\).

First step:

Second step:

The exact value at \(t=0.2\) is approximately \(1.22140\).

Example: Euler Stability

Consider

Euler gives

Stability requires

Thus,