The Central Idea
Consider an initial-value problem
\[ \boxed{ y'=f(x,y), \qquad y(x_0)=y_0 }. \]The differential equation tells us the slope of the solution at each point.
A numerical method uses these slopes to estimate the solution at a sequence of points:
\[ x_0,\ x_1,\ x_2,\ldots \]producing approximations
\[ y_0,\ y_1,\ y_2,\ldots \]to the exact values
\[ y(x_0),\ y(x_1),\ y(x_2),\ldots \]Why Numerical Methods Are Needed
Some differential equations have exact solutions that are difficult to derive.
Others have no convenient closed-form solution in elementary functions.
Numerical methods are useful when:
- an exact symbolic solution is unavailable;
- only numerical values are needed;
- the model is too large or nonlinear for direct analysis;
- the equation is part of a computer simulation;
- parameter studies require many repeated solutions.
The Step Size
Suppose we move through the independent variable using a fixed step size
\[ h. \]Then
\[ x_{n+1}=x_n+h. \]The choice of \(h\) controls how far the numerical method moves during each update.
Smaller steps usually improve accuracy but require more computation.
Euler's Method
Euler's method is the simplest widely used numerical method for an initial-value problem.
Starting from
\[ (x_n,y_n), \]the differential equation provides the slope
\[ f(x_n,y_n). \]Euler's update is
\[ \boxed{ y_{n+1} = y_n+h f(x_n,y_n) }. \]Where Euler's Method Comes From
The derivative gives the local linear approximation
\[ y(x+h) \approx y(x)+hy'(x). \]Since
\[ y'=f(x,y), \]we obtain
\[ y(x+h) \approx y(x)+hf(x,y). \]Euler's method repeatedly applies this tangent-line approximation.
Euler Update in Two Steps
Euler's method can be viewed as:
- find the current slope;
- move along that slope for one step.
In symbols,
\[ m_n=f(x_n,y_n), \]then
\[ y_{n+1}=y_n+hm_n. \]Example: One Euler Step
Consider
\[ y'=x+y, \qquad y(0)=1. \]Use
\[ h=0.1. \]At
\[ (x_0,y_0)=(0,1), \]the slope is
\[ f(0,1)=1. \]Therefore,
\[ y_1 = 1+0.1(1) = 1.1. \]Thus,
\[ \boxed{ y(0.1)\approx1.1 }. \]Example: A Second Euler Step
Continue from
\[ (0.1,1.1). \]The new slope is
\[ f(0.1,1.1)=1.2. \]Therefore,
\[ y_2 = 1.1+0.1(1.2) = 1.22. \]Thus,
\[ \boxed{ y(0.2)\approx1.22 }. \]Euler Table
Numerical work is often organized in a table.
| \(n\) | \(x_n\) | \(y_n\) | \(f(x_n,y_n)\) | \(y_{n+1}\) |
|---|---|---|---|---|
| 0 | \(0\) | \(1\) | \(1\) | \(1.1\) |
| 1 | \(0.1\) | \(1.1\) | \(1.2\) | \(1.22\) |
A table makes the iterative structure easy to verify.
Local Truncation Error
A numerical formula approximates the exact solution over each step.
The error introduced during one step, assuming the starting value is exact, is called the local truncation error.
Euler's method has local truncation error of order
\[ \boxed{O(h^2)}. \]Global Error
Errors accumulate over many steps.
The difference between the true solution and the numerical approximation after many updates is called the global error.
Under suitable smoothness and stability conditions, Euler's method has global error of order
\[ \boxed{O(h)}. \]Halving the step size therefore tends to reduce the global error by roughly a factor of two.
Order of a Numerical Method
A numerical method is said to have order \(p\) when its global error behaves approximately like
\[ Ch^p \]for sufficiently small \(h\).
Higher-order methods can achieve much greater accuracy for the same step size.
Improved Euler Method
Euler's method uses only the slope at the beginning of the step.
The improved Euler method, also called Heun's method in a common form, estimates a second slope and averages the two.
Predictor Step
First compute the Euler prediction:
\[ \boxed{ y^* = y_n+h f(x_n,y_n) }. \]This predicts the value at
\[ x_{n+1}=x_n+h. \]Corrector Step
Compute the slope at the predicted endpoint:
\[ f(x_{n+1},y^*). \]Average the initial and predicted slopes:
\[ \boxed{ y_{n+1} = y_n + \frac{h}{2} \left[ f(x_n,y_n) + f(x_{n+1},y^*) \right] }. \]Example: Improved Euler
Use one improved Euler step for
\[ y'=x+y, \qquad y(0)=1, \qquad h=0.1. \]The initial slope is
\[ k_1=f(0,1)=1. \]Euler predicts
\[ y^*=1+0.1(1)=1.1. \]The endpoint slope is
\[ k_2=f(0.1,1.1)=1.2. \]Therefore,
\[ y_1 = 1+\frac{0.1}{2}(1+1.2). \]Thus,
\[ \boxed{ y(0.1)\approx1.11 }. \]Why Improved Euler Is Better
Euler's method assumes the starting slope represents the entire interval.
Improved Euler uses information from both ends of the step.
This better captures curvature.
The method has global error of order
\[ \boxed{O(h^2)}. \]The Midpoint Method
Another second-order method estimates the slope near the middle of the step.
First compute
\[ k_1=f(x_n,y_n). \]Estimate the midpoint state:
\[ y_{\text{mid}} = y_n+\frac{h}{2}k_1. \]Then compute
\[ k_2 = f\left( x_n+\frac{h}{2}, y_{\text{mid}} \right). \]Update with
\[ \boxed{ y_{n+1}=y_n+hk_2 }. \]Runge-Kutta Methods
Runge-Kutta methods use several slope evaluations within each step to achieve higher accuracy without requiring higher derivatives of \(f\).
The most famous example is the classical fourth-order Runge-Kutta method, often called RK4.
The RK4 Slopes
For
\[ y'=f(x,y), \]define
\[ k_1=f(x_n,y_n), \] \[ k_2 = f\left( x_n+\frac{h}{2}, y_n+\frac{h}{2}k_1 \right), \] \[ k_3 = f\left( x_n+\frac{h}{2}, y_n+\frac{h}{2}k_2 \right), \] \[ k_4 = f(x_n+h,y_n+hk_3). \]The RK4 Update
Combine the four slopes using the weighted average
\[ \boxed{ y_{n+1} = y_n + \frac{h}{6} \left( k_1+2k_2+2k_3+k_4 \right) }. \]RK4 has global error of order
\[ \boxed{O(h^4)}. \]Why RK4 Is So Effective
RK4 samples the slope:
- at the beginning of the step;
- twice near the middle;
- at the end.
The weighted combination captures the changing slope much more accurately than a single tangent estimate.
Example: One RK4 Step
Consider
\[ y'=x+y, \qquad y(0)=1, \qquad h=0.1. \]First,
\[ k_1=f(0,1)=1. \]Next,
\[ k_2 = f(0.05,1.05) = 1.10. \]Then,
\[ k_3 = f(0.05,1.055) = 1.105. \]Finally,
\[ k_4 = f(0.1,1.1105) = 1.2105. \]Therefore,
\[ y_1 = 1 + \frac{0.1}{6} \left( 1+2(1.10)+2(1.105)+1.2105 \right). \]This gives approximately
\[ \boxed{ y(0.1)\approx1.11034 }. \]Comparing the Methods
For the same step size, a higher-order method usually gives a more accurate approximation.
A rough hierarchy is:
\[ \text{Euler} \quad\to\quad \text{second-order methods} \quad\to\quad \text{RK4}. \]However, higher accuracy per step comes at the cost of more slope evaluations.
Accuracy versus Computational Cost
A numerical method must balance:
- step size;
- order of accuracy;
- number of function evaluations;
- memory;
- stability;
- problem complexity.
The most efficient method is not always the method with the highest formal order.
Step-Size Refinement
A practical way to estimate numerical reliability is to compute the solution using one step size and then repeat with a smaller step size.
For example, compare results using
\[ h \]and
\[ \frac{h}{2}. \]If the results agree closely, the approximation may be approaching convergence.
Convergence
A numerical method is convergent if its approximations approach the exact solution as
\[ h\to0 \]under appropriate conditions.
Convergence is one of the central requirements of a reliable numerical method.
Consistency
A numerical method is consistent when its one-step formula approaches the original differential equation as the step size becomes small.
Informally, the discrete method must represent the correct continuous rate law in the small-step limit.
Stability
Stability concerns how numerical errors behave as the computation advances.
A stable method does not allow small perturbations or rounding errors to grow uncontrollably in situations where the true solution should remain well behaved.
A Standard Stability Test Equation
A common model for studying numerical stability is
\[ \boxed{ y'=\lambda y }. \]Its exact solution is
\[ y=y_0e^{\lambda t}. \]When
\[ \lambda<0, \]the exact solution decays.
Euler Stability for the Test Equation
Apply Euler's method:
\[ y_{n+1} = y_n+h\lambda y_n. \]Therefore,
\[ \boxed{ y_{n+1} = (1+h\lambda)y_n }. \]For the numerical values to decay, we need
\[ |1+h\lambda|<1. \]Euler Step-Size Restriction
If
\[ \lambda<0, \]the Euler stability condition becomes
\[ \boxed{ 0Example: Numerical Instability
Consider
\[ y'=-10y. \]The exact solution decays rapidly.
Euler's method gives
\[ y_{n+1}=(1-10h)y_n. \]Stability requires
\[ |1-10h|<1. \]Therefore,
\[ \boxed{0Stiff Differential Equations
A problem is called stiff when some components decay or vary on much faster time scales than others, forcing explicit methods to use very small step sizes for stability.
Stiffness is common in:
- chemical kinetics;
- electrical circuits;
- control systems;
- reaction networks;
- multiscale physical models.
Implicit Euler Method
One simple implicit method is
\[ \boxed{ y_{n+1} = y_n + h f(x_{n+1},y_{n+1}) }. \]Unlike ordinary Euler's method, the new value appears on both sides of the equation.
This may require solving an algebraic equation at each step.
Why Use Implicit Methods?
Implicit methods can have much better stability properties for stiff problems.
The tradeoff is additional computational work because each update may require solving a nonlinear or linear system.
Adaptive Step Sizes
A fixed step size is simple, but many solutions change at uneven rates.
An adaptive method changes \(h\) during the calculation.
It may:
- use smaller steps where the solution changes rapidly;
- use larger steps where the solution is smooth;
- estimate local error automatically;
- try to meet a requested tolerance efficiently.
Error Tolerances
Numerical software often uses two types of tolerance:
- absolute tolerance, controlling absolute error;
- relative tolerance, controlling error relative to the size of the solution.
A typical error scale combines both ideas.
Embedded Runge-Kutta Methods
Modern adaptive solvers often compute two approximations of different orders using largely shared slope evaluations.
Their difference provides an estimate of local error.
The step size is then increased or decreased automatically.
Numerical Methods for Systems
Numerical methods extend directly to vector-valued systems
\[ \boxed{ \mathbf{y}'=\mathbf{F}(t,\mathbf{y}) }. \]Euler's method becomes
\[ \boxed{ \mathbf{y}_{n+1} = \mathbf{y}_n + h\mathbf{F}(t_n,\mathbf{y}_n) }. \]Example: Euler Method for a System
Consider
\[ \begin{aligned} x'&=x+y,\\ y'&=x-y, \end{aligned} \]with
\[ x(0)=1, \qquad y(0)=2. \]Use
\[ h=0.1. \]At the initial state,
\[ x'=3, \qquad y'=-1. \]Therefore,
\[ x_1=1+0.1(3)=1.3 \]and
\[ y_1=2+0.1(-1)=1.9. \]Thus,
\[ \boxed{ (x(0.1),y(0.1)) \approx (1.3,1.9) }. \]Higher-Order Equations as Numerical Systems
Before applying a standard numerical solver, a higher-order differential equation is often rewritten as a first-order system.
For example,
\[ y''+3y'+2y=0 \]becomes
\[ \begin{aligned} y'&=v,\\ v'&=-3v-2y. \end{aligned} \]Numerical methods can then advance the state vector
\[ \begin{bmatrix} y\\ v \end{bmatrix}. \]Numerical Phase Portraits
Numerical solution of a system produces sequences
\[ (x_n,y_n). \]Plotting these state pairs creates an approximate phase trajectory.
Repeating the calculation for several initial states builds a numerical phase portrait.
Numerical Boundary-Value Problems
Boundary-value problems may require specialized numerical methods.
Common ideas include:
- shooting methods;
- finite-difference methods;
- collocation methods.
The shooting method converts a boundary problem into repeated initial-value problems.
Finite-Difference Idea
Derivatives can be approximated directly on a grid.
For example,
\[ y'(x_i) \approx \frac{y_{i+1}-y_i}{h} \]and
\[ \boxed{ y''(x_i) \approx \frac{y_{i+1}-2y_i+y_{i-1}}{h^2} }. \]Replacing derivatives with finite differences converts a differential equation into algebraic equations for the grid values.
Round-Off Error
Computers store numbers with finite precision.
This produces round-off error.
In many ordinary calculations, truncation error dominates.
But after extremely many steps or in poorly conditioned computations, round-off can become significant.
Truncation Error versus Round-Off Error
Truncation error comes from replacing an exact mathematical process with an approximation.
Round-off error comes from finite numerical precision.
Reducing the step size lowers truncation error, but it also increases the number of operations.
Verification with a Known Exact Solution
When an exact solution is known, a numerical method can be tested by comparing
\[ y_n \]with
\[ y(x_n). \]The absolute error is
\[ \boxed{ |y(x_n)-y_n| }. \]Relative Error
When the exact value is nonzero, relative error may be measured by
\[ \boxed{ \frac{|y(x_n)-y_n|}{|y(x_n)|} }. \]This scales the error relative to the size of the quantity being estimated.
Residual Checking
When the exact solution is unknown, one useful check is to ask how well the numerical approximation satisfies the differential equation.
This may involve estimating derivatives of the numerical solution and checking the residual
\[ y'-f(x,y). \]A small residual is useful evidence, though it does not by itself guarantee small solution error.
Conservation Checks
If the continuous model conserves a quantity such as
\[ E(x,y) \]or
\[ x+y, \]the numerical solution can be checked to see how well it preserves that invariant.
Large artificial drift may indicate that the step size or method is unsuitable.
Positivity and Physical Constraints
Some state variables must remain nonnegative.
Examples include:
- population;
- mass;
- concentration;
- probability;
- inventory.
A numerical method that produces physically impossible negative values may require a smaller step size or a method designed to preserve the constraint.
Event Detection
Numerical solvers often need to determine when a solution reaches a threshold.
Examples include:
- population reaching a target;
- position reaching a boundary;
- temperature crossing a critical value;
- concentration falling below a limit.
Modern solvers can locate such events between time steps.
Common Error: Using the New Slope Too Early in Euler's Method
Euler's method uses
\[ f(x_n,y_n), \]not
\[ f(x_{n+1},y_{n+1}). \]The latter corresponds to an implicit method, not ordinary Euler stepping.
Common Error: Forgetting to Update \(x\)
After each step,
\[ x_{n+1}=x_n+h. \]The next slope must be evaluated using the updated independent variable.
Common Error: Reusing the Old \(y\)-Value
Each step begins from the newest approximation.
After computing
\[ y_{n+1}, \]use it in the next slope calculation.
Common Error: Assuming Smaller \(h\) Always Fixes Everything
A smaller step often improves accuracy, but numerical difficulty can also come from:
- stiffness;
- poor conditioning;
- unstable algorithms;
- singular behavior;
- incorrect model formulation.
Step-size reduction is important, but it is not a universal remedy.
Common Error: Confusing Method Order with Differential-Equation Order
A fourth-order Runge-Kutta method does not mean the differential equation is fourth order.
The word order has different meanings:
- differential-equation order refers to the highest derivative;
- numerical-method order refers to how error scales with step size.
Common Error: Reporting Too Many Digits
A numerical approximation should not imply more accuracy than the method and step size justify.
Report enough digits to reflect the actual reliability of the calculation.
Common Error: Ignoring Stability
A method can be accurate in theory yet fail numerically when the step size is outside its stability region.
Always distinguish:
- accuracy;
- convergence;
- stability.
A Reliable Numerical Workflow
- Write the differential equation as a first-order system if necessary.
- Specify the initial state and interval.
- Choose a numerical method.
- Choose an initial step size.
- Compute the required slopes carefully.
- Advance the solution one step at a time.
- Repeat with a smaller step size when practical.
- Check stability and physical constraints.
- Compare against exact solutions or conserved quantities when available.
- Estimate or monitor error.
- Use adaptive stepping for difficult or uneven problems.
- Interpret the numerical solution in the context of the model.
Example: Full Euler Workflow
Approximate the solution of
\[ y'=2x-y, \qquad y(0)=1, \]at
\[ x=0.2 \]using Euler's method with
\[ h=0.1. \]Step 1:
\[ f(0,1)=-1. \]Therefore,
\[ y_1=1+0.1(-1)=0.9. \]Complete the Full Euler Workflow
Now
\[ (x_1,y_1)=(0.1,0.9). \]Compute
\[ f(0.1,0.9)=0.2-0.9=-0.7. \]Therefore,
\[ y_2 = 0.9+0.1(-0.7) = 0.83. \]Thus,
\[ \boxed{ y(0.2)\approx0.83 }. \]Example: Full Improved Euler Workflow
Use one improved Euler step for
\[ y'=2x-y, \qquad y(0)=1, \qquad h=0.1. \]Initial slope:
\[ k_1=-1. \]Predictor:
\[ y^*=1+0.1(-1)=0.9. \]Endpoint slope:
\[ k_2=f(0.1,0.9)=-0.7. \]Correct:
\[ y_1 = 1+\frac{0.1}{2}(-1-0.7). \]Therefore,
\[ \boxed{ y(0.1)\approx0.915 }. \]Practice 1: One Euler Step
For
\[ y'=x+2y, \qquad y(0)=1, \]use Euler's method with
\[ h=0.1 \]to approximate \(y(0.1)\).
Practice 2: Two Euler Steps
Continue Practice 1 to approximate
\[ y(0.2). \]Practice 3: Improved Euler
For
\[ y'=x+y, \qquad y(0)=1, \]use one improved Euler step with
\[ h=0.2. \]Practice 4: Euler Stability
For
\[ y'=-4y, \]what range of positive step sizes makes Euler's method stable for the test equation?
Practice 5: Method Order
What is the global order of Euler's method?
What is the global order of classical RK4?
Practice 6: Euler for a System
For
\[ \begin{aligned} x'&=x-y,\\ y'&=x+y, \end{aligned} \]with
\[ x(0)=2, \qquad y(0)=1, \]use one Euler step with
\[ h=0.1. \]Practice 7: Finite Difference
Write the centered finite-difference approximation for
\[ y''(x_i). \]Practice 8: Exact versus Numerical
If the exact value is \(2.500\) and the numerical approximation is \(2.470\), find the absolute error.
Practice 9: Step-Size Reasoning
If Euler's method has global error approximately proportional to \(h\), what happens approximately to the error when \(h\) is cut in half?
Practice 10: Stiffness
Why can a stiff differential equation force an explicit method to use a very small step size even when the solution appears smooth?
Practice Answers
Practice 1
At
\[ (0,1), \]the slope is
\[ f(0,1)=0+2=2. \]Therefore,
\[ y_1 = 1+0.1(2) = 1.2. \]Thus,
\[ \boxed{ y(0.1)\approx1.2 }. \]Practice 2
At
\[ (0.1,1.2), \]the slope is
\[ f(0.1,1.2) = 0.1+2.4 = 2.5. \]Therefore,
\[ y_2 = 1.2+0.1(2.5) = 1.45. \]Thus,
\[ \boxed{ y(0.2)\approx1.45 }. \]Practice 3
Initial slope:
\[ k_1=f(0,1)=1. \]Predictor:
\[ y^*=1+0.2(1)=1.2. \]Endpoint slope:
\[ k_2=f(0.2,1.2)=1.4. \]Therefore,
\[ y_1 = 1+\frac{0.2}{2}(1+1.4). \]Thus,
\[ \boxed{ y(0.2)\approx1.24 }. \]Practice 4
Euler stability requires
\[ |1+h\lambda|<1. \]Here
\[ \lambda=-4. \]Therefore,
\[ |1-4h|<1. \]This gives
\[ \boxed{ 0Euler's method has global order
\[ \boxed{1}. \]Classical RK4 has global order
\[ \boxed{4}. \]Practice 6
At
\[ (x,y)=(2,1), \]we have
\[ x'=2-1=1 \]and
\[ y'=2+1=3. \]Therefore,
\[ x_1=2+0.1(1)=2.1 \]and
\[ y_1=1+0.1(3)=1.3. \]Thus,
\[ \boxed{ (x(0.1),y(0.1)) \approx (2.1,1.3) }. \]Practice 7
The centered second-difference formula is
\[ \boxed{ y''(x_i) \approx \frac{y_{i+1}-2y_i+y_{i-1}}{h^2} }. \]Practice 8
The absolute error is
\[ |2.500-2.470| = 0.030. \]Therefore,
\[ \boxed{0.030}. \]Practice 9
Since Euler's global error is approximately proportional to \(h\), halving \(h\) approximately halves the error.
Therefore,
\[ \boxed{ \text{error is reduced by about a factor of }2 }. \]Practice 10
A stiff system contains rapidly decaying or rapidly varying modes that can destabilize an explicit method unless the step size is small.
The required step size may therefore be controlled by stability rather than by the visible smoothness of the dominant solution.
What to Remember
Numerical methods approximate differential-equation solutions by replacing continuous evolution with a sequence of discrete updates.
Euler's method is
\[ \boxed{ y_{n+1}=y_n+h f(x_n,y_n) }. \]- Euler's method is first-order accurate globally;
- improved Euler and midpoint methods use better slope estimates;
- RK4 uses four slope evaluations and is fourth-order accurate globally;
- smaller step sizes generally improve accuracy;
- stability can impose stricter step-size limits than accuracy;
- stiff equations often require implicit or specialized methods;
- adaptive methods adjust step size automatically;
- systems are solved numerically by updating vectors of state variables;
- verification should include error checks, physical constraints, and invariants.
The central computational idea is
\[ \boxed{ \text{current state} + \text{estimated change} \longrightarrow \text{next state} }. \]Where We Go Next
Ordinary differential equations describe change with respect to one independent variable. Many systems, however, depend on several independent variables at once, such as space and time.
The next lesson, Introductory Partial Differential Equations, introduces partial derivatives in dynamic models, the heat equation, wave equation, Laplace equation, boundary and initial data, and the basic idea of separating variables.