Numerical methods

While second-degree polynomials can be solved through factoring, completing the square, or applying the quadratic formula, straightforward analytical methods for cubic and higher-order equations exist only in a narrow set of cases. However, such higher-order equations are typical in engineering applications, so finding approximate solutions to the equations would be useful.

This module examines three fundamental areas of numerical methods: root finding using Newton’s method, numerical integration techniques, and numerical solutions of ordinary differential equations.

Root extraction: Newton’s method

Newton’s method is a particular form of fixed-point iteration. The goal is to find a value $x$ such that $f(x)=0$.

All fixed-point techniques require a starting point. Preferably, the starting point will be close to the actual root. And, while Newton’s method converges quickly, it requires the function to be continuously differentiable.

Given an initial guess of $x_0$, the next estimate for the solution to $f(x)=0$ is:

$$x_{i+1}=x_i-{\left[\frac{f(x)}{\frac{df(x)}{dx}}\right]}_{x=x_i}$$

or equivalently:

$$x_{i+1}=x_i-\frac{f(x_i)}{f^{\prime }(x_i)}$$

Newton’s method can be understood geometrically as follows:

1. Start with an initial guess $x_0$
2. Draw the tangent line to the curve $y=f(x)$ at the point $(x_0,f(x_0))$
3. Find where this tangent line intersects the x-axis; this intersection is $x_1$
4. Repeat the process using $x_1$ as the new starting point

The tangent line at $x_i$ has slope $f^{\prime }(x_i)$ and passes through $(x_i,f(x_i))$. Setting $y=0$ and solving for $x$ gives the next approximation $x_{i+1}$.

At each iteration ($i=0,1,2$, etc.), the method estimates the root. The maximum error is determined by looking at how much the estimate changes after each iteration. If the change between the previous and current estimates (representing the magnitude of error in the estimate) is too large, the current estimate is used as the independent variable for the subsequent iteration.

Convergence criterion:

$$|x_{i+1}-x_i|<ϵ$$

where $ϵ$ is a predetermined tolerance.

Properties of Newton’s method:

• Quadratic convergence: when it converges, Newton’s method typically doubles the number of correct digits with each iteration.
• Requires differentiability: the function must be continuously differentiable.
• Sensitive to starting point: poor initial guesses may lead to divergence or convergence to an unintended root.

Numerical integration

Numerical integration (also called numerical quadrature) is the process of approximating the value of a definite integral when an analytical solution is difficult or impossible to obtain. We present three common methods: Euler’s rule, the trapezoidal rule, and Simpson’s rule.

General setup

For all numerical integration methods, we divide the interval $[a,b]$ into $n$ subintervals of equal width:

$$\Delta x=\frac{b-a}{n}$$

The partition points are:

$$x_k=a+k\Delta x,k=0,1,2,\dots ,n$$

Euler’s rule (left endpoint rule)

Euler’s rule is the simplest numerical integration technique. It approximates the integral by summing rectangles whose heights are determined by the function value at the left endpoint of each subinterval.

$${\int }_a^{b}f(x)dx\approx \Delta x\sum _{k=0}^{n-1}f(a+k\Delta x)$$

Geometric interpretation: Each rectangle has width $\Delta x$ and height $f(x_k)$, where $x_k$ is the left endpoint of the $k$th subinterval.

Characteristics:

• Simplest method to implement.
• Generally less accurate than trapezoidal or Simpson’s rule.
• Error is proportional to $\Delta x$ (first-order accuracy).

Trapezoidal rule

The trapezoidal rule improves upon Euler’s rule by approximating the area under the curve using trapezoids instead of rectangles.

Formula for $n=1$ (single trapezoid):

$${\int }_a^{b}f(x)dx\approx \Delta x\left[\frac{f(a)+f(b)}{2}\right]$$

Formula for $n>1$ (multiple trapezoids):

$${\int }_a^{b}f(x)dx\approx \frac{\Delta x}{2}\left[f(a)+2\sum _{k=1}^{n-1}f(a+k\Delta x)+f(b)\right]$$

Alternative notation:

$$A=\frac{d}{2}\left(h_0+2\sum _{i=1}^{n-1}{h}_i+h_n\right)$$

where $d$ is the width of each cell and $h_i$ represents the function values (heights) at each partition point.

Geometric interpretation: Each trapezoid has parallel sides of lengths $f(x_k)$ and $f(x_{k+1})$, with width $\Delta x$. The area of each trapezoid is $\frac{\Delta x}{2}[f(x_k)+f(x_{k+1})]$.

Sidenote

When to use the trapezoidal rule

Areas of sections with irregular boundaries cannot be determined precisely, so approximation methods must be used. If the irregular side can be divided into a series of $n$ cells of equal width, and if the irregular side of each cell is fairly straight, the trapezoidal rule is appropriate.

Characteristics:

• More accurate than Euler’s rule.
• Error is proportional to $(\Delta x{)}^2$ (second-order accuracy).
• Exact for linear functions.

Simpson’s rule (parabolic rule)

Simpson’s rule achieves higher accuracy by approximating the function with parabolic arcs instead of straight lines.

Formula for $n=2$ (single parabola):

$${\int }_a^{b}f(x)dx\approx \left(\frac{b-a}{6}\right)\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]$$

Formula for general $n$ (must be even):

$${\int }_a^{b}f(x)dx\approx \frac{\Delta x}{3}\left[f(a)+2\sum _{\begin{array}{c}k=2,4,6,\dots \end{array}}^{n-2}f(a+k\Delta x)+4\sum _{\begin{array}{c}k=1,3,5,\dots \end{array}}^{n-1}f(a+k\Delta x)+f(b)\right]$$

where:

$$\Delta x=\frac{b-a}{n}$$

Alternative notation:

$$A=\frac{d}{3}\left(h_0+2\sum _{\begin{array}{c}i\text{ even}\\ i=2\end{array}}^{n-2}{h}_i+4\sum _{\begin{array}{c}i\text{ odd}\\ i=1\end{array}}^{n-1}{h}_i+h_n\right)$$

Important requirement: the value of $n$ must be even to use Simpson’s rule.

Geometric interpretation: Simpson’s rule fits a parabola through each consecutive set of three points and integrates the parabola exactly.

Sidenote

When to use Simpson's rule

If the irregular side of each cell is curved (parabolic), Simpson’s rule should be used instead of the trapezoidal rule.

Characteristics:

• More accurate than both Euler’s rule and the trapezoidal rule.
• Error is proportional to $(\Delta x{)}^4$ (fourth-order accuracy).
• Exact for polynomials of degree 3 or less.

Numerical solutions of ordinary differential equations

Many differential equations arising in engineering and science cannot be solved analytically. In such cases, numerical methods provide approximate solutions at discrete points.

Consider a first-order ordinary differential equation of the form:

$$\frac{dx}{dt}=f(x,t),x(0)=x_0$$

The goal is to find approximate values of $x(t)$ at discrete time points $t_k=k\Delta t$.

Euler’s approximation

Euler’s approximation (also called Euler’s method or forward Euler method) is the simplest method for numerically solving ordinary differential equations. It estimates the value of a function given the value and slope of the function at an adjacent location.

Given the differential equation

$$\frac{d}{dt}x(t)=f(x,t),$$

various points along the solution $x(t)$ may be approximated as:

$$x[(k+1)\Delta t]\cong x(k\Delta t)+\Delta t\cdot f[x(k\Delta t),k\Delta t]$$

$$x[(k+1)\Delta t]\cong x(k\Delta t)+\Delta t\cdot f[x(k\Delta t)]$$

$$x_{k+1}=x_k+\Delta t\cdot \frac{d{x}_k}{dt}$$

$$x_{k+1}=x_k+\Delta t\cdot f(x_k,t_k)$$

Geometric interpretation: The simplicity of the concept is illustrated by writing Euler’s approximation in terms of the traditional two-dimensional $x$-$y$ coordinate system:

$$y(x_2)=y(x_1)+(x_2-x_1)y^{\prime }(x_1)$$

This formula says: starting from a known point $(x_1,y(x_1))$, move along the tangent line (with slope $y^{\prime }(x_1)$) to estimate the function value at a nearby point $x_2$.

Applicability: Euler’s approximation applies to a differential equation of the form:

$$f(x,t)=\frac{dx}{dt},x(0)=x_0$$

The model applies to a general time $k\Delta t$, and it applies when $f(x)=dx/dt$ and can be expressed recursively.

As long as the derivative can be evaluated, Euler’s method can be used to predict the value of any function whose values are limited to discrete, sequential points in time or space (e.g., a difference equation).

Euler’s approximation for ordinary differential equations provides a straightforward method for obtaining numerical solutions when analytical methods fail. While it is exact for linear systems, care must be taken with nonlinear systems where error accumulates.

These numerical methods form the foundation for more advanced computational techniques used throughout engineering, physics, and applied mathematics. Modern software packages implement sophisticated versions of these algorithms, but understanding the basic principles remains essential for proper application and interpretation of results.