Agricultural Statistics & Data Science · Foundations · also Numerical Analysis, Calculus of Finite Differences, Mathematics for Statistics (part)
Ordinary differential equations, interpolation, numerical differentiation and integration and root-finding: what AIEEA-PG Group 07 names explicitly, and what growth models, life tables and iterative estimation (Newton-Raphson for maximum likelihood) rely on.
In the exams: AIEEA-PG Group 07 Units II-III (heavy: differential equations, finite differences, interpolation, numerical integration, Newton-Raphson named); AICE-JRF Computer Application Unit 10 numerical analysis (heavy)
Statistics in agriculture constantly meets equations that cannot be solved by hand. Crop growth curves are solutions of differential equations; life tables, population projections and price cycles are difference equations; tables of reference values must be interpolated; areas under response curves and probabilities of non-standard distributions must be integrated numerically; and almost every maximum likelihood estimate beyond the simplest models is found by iterating Newton–Raphson or a relative of it. This subject teaches the classical tools for these tasks: ordinary differential equations and their analytical solutions, the calculus of finite differences, interpolation, numerical differentiation and integration, difference equations, root-finding, the solution of linear systems, numerical solution of differential equations, and, because it appears alongside them in many syllabi, the coordinate geometry of conics.
The emphasis is on understanding what each method does, when it works, how accurate it is and how its errors behave as the step size shrinks, so that the reader can choose a method sensibly and recognise when software output is unreliable. Each chapter includes short worked examples of the kind examinations set, and points to the statistical problems where the method is used: logistic and Gompertz growth curves, the delta and Newton–Raphson methods, Simpson's rule for probabilities, Gauss–Seidel iteration for large mixed-model equations, and Runge–Kutta solution of crop and disease models.
Master's admission examinations in agricultural statistics include differential equations, finite differences, interpolation, numerical integration and the Newton–Raphson method by name, usually as short numerical problems; doctoral examinations in computer applications for agriculture test numerical analysis more deeply, including error analysis and iterative methods; and graduate entrance examinations in statistics and mathematics include differential equations and basic numerical methods.
An ordinary differential equation (ODE) relates a function of one variable to its derivatives; its order is that of the highest derivative. A first-order equation dy/dx = f(x, y) is separable if f factors as g(x)h(y); then ∫dy/h(y) = ∫g(x)dx. Exponential growth dN/dt = rN separates to give N = N₀e^(rt), with doubling time ln 2/r, and decay with rate k has half-life ln 2/k. An equation M(x, y)dx + N(x, y)dy = 0 is exact if ∂M/∂y = ∂N/∂x; then a function F with ∂F/∂x = M and ∂F/∂y = N exists and the solution is F(x, y) = C. A non-exact equation can sometimes be made exact by an integrating factor: if (∂M/∂y − ∂N/∂x)/N depends on x alone, e^(∫ that dx) is one. The linear first-order equation dy/dx + P(x)y = Q(x) has integrating factor e^(∫P dx), and multiplying through gives d/dx[y e^(∫P dx)] = Q e^(∫P dx). For dy/dt + 2y = 4 with y(0) = 0, the integrating factor is e^(2t), and y = 2(1 − e^(−2t)), which approaches the equilibrium 2; at t = 1, y ≈ 1.73.
A Bernoulli equation dy/dx + P(x)y = Q(x)yⁿ, with n ≠ 0, 1, becomes linear in v = y^(1−n). The logistic growth model dN/dt = rN(1 − N/K) is of this type (n = 2) and is also separable; its solution is N(t) = K/(1 + Ae^(−rt)) with A = (K − N₀)/N₀. With K = 1,000, N₀ = 100 and r = 0.5, N(2) = 1,000/(1 + 9e^(−1)) ≈ 232. The logistic curve is sigmoid, with its inflection at N = K/2, where growth is fastest. The Gompertz model dN/dt = rN ln(K/N), asymmetric with inflection at K/e, and the monomolecular (Mitscherlich) model dN/dt = r(K − N), which has no inflection and describes diminishing returns to fertiliser, are other growth laws used in agriculture and forestry. Newton's law of cooling, dT/dt = −k(T − Tₐ), and first-order kinetics of drug elimination, pesticide decay and soil organic matter decomposition have the same linear form.
Not every initial value problem has a unique solution. The Picard–Lindelöf theorem states that if f(x, y) and ∂f/∂y are continuous near (x₀, y₀), the problem dy/dx = f(x, y), y(x₀) = y₀ has a unique solution near x₀. The equation dy/dx = 3y^(2/3) with y(0) = 0 has both y = 0 and y = x³ as solutions, because ∂f/∂y is unbounded at y = 0. Picard's method of successive approximations, yₙ₊₁(x) = y₀ + ∫ from x₀ to x of f(t, yₙ(t))dt, both proves the theorem and gives an approximation method.
Choosing among exponential, monomolecular, logistic and Gompertz models is partly a matter of biology and partly of data. Exponential growth fits only the early phase of growth, before resources become limiting. The monomolecular model suits responses that rise fastest at first, such as yield response to nutrients. The logistic model gives a symmetric S-shape and the Gompertz an asymmetric one with a longer approach to the asymptote, often better for animal and tree growth. Richards' generalisation adds a shape parameter that includes several of these as special cases. Fitted models should be compared by residual plots and information criteria, and the asymptote K should be biologically plausible rather than an extrapolation far beyond the data.
A linear equation of order n with constant coefficients, aₙy⁽ⁿ⁾ + … + a₁y′ + a₀y = F(x), has general solution y = y_c + y_p, the complementary function (general solution of the homogeneous equation F = 0) plus any particular integral. Writing D for d/dx, the homogeneous equation is f(D)y = 0, and trying y = e^(mx) gives the auxiliary (characteristic) equation f(m) = 0. If its roots are real and distinct, y_c = c₁e^(m₁x) + c₂e^(m₂x) + …; a repeated root m of multiplicity k contributes (c₁ + c₂x + … + c_k x^(k−1))e^(mx); and a complex pair α ± iβ contributes e^(αx)(c₁ cos βx + c₂ sin βx). For y″ − 3y′ + 2y = 0 the roots are 1 and 2, so y = c₁eˣ + c₂e^(2x); for y″ − 2y′ + y = 0 the double root 1 gives (c₁ + c₂x)eˣ; for y″ + 4y = 0 the roots ±2i give c₁ cos 2x + c₂ sin 2x, an oscillation with period π. Damped oscillations, α < 0, describe systems returning to equilibrium, such as predator–prey models linearised around equilibrium.
The particular integral can be found by undetermined coefficients, guessing a form like the right-hand side, or by the inverse operator 1/f(D). For F(x) = e^(ax) with f(a) ≠ 0, y_p = e^(ax)/f(a): for y″ − 3y′ + 2y = e^(3x), f(3) = 9 − 9 + 2 = 2, so y_p = e^(3x)/2. If f(a) = 0 (resonance), multiply by x: if a is a simple root, y_p = x e^(ax)/f′(a). For F = sin bx or cos bx, replace D² by −b² in f(D). For polynomial F, expand 1/f(D) as a power series in D and stop at the degree of F; for y″ + y = x², y_p = (1 − D² + …)x² = x² − 2. Variation of parameters handles general F. Initial conditions fix the constants after y_c + y_p is formed, not before.
Higher-order equations can be rewritten as systems of first-order equations, and systems of linear equations x′ = Ax are solved through the eigenvalues and eigenvectors of A: solutions combine terms e^(λt)v. Their stability depends on the eigenvalues' real parts: all negative means trajectories decay to equilibrium. Compartment models of nutrient flows in soil, pharmacokinetics of veterinary drugs and Lotka–Volterra models of pests and natural enemies are systems of this kind, and linearising a non-linear system near an equilibrium uses the Jacobian matrix.
For a function tabulated at equally spaced points x₀, x₀ + h, x₀ + 2h, …, the forward difference is Δf(x) = f(x + h) − f(x), the backward difference ∇f(x) = f(x) − f(x − h), and the central difference δf(x) = f(x + h/2) − f(x − h/2). The shift operator E is defined by Ef(x) = f(x + h), the averaging operator μf(x) = [f(x + h/2) + f(x − h/2)]/2, and D the derivative. The operators are related by Δ = E − 1, ∇ = 1 − E⁻¹, δ = E^(1/2) − E^(−1/2), Δ = E∇ = ∇E, and E = e^(hD), so hD = ln(1 + Δ) = Δ − Δ²/2 + Δ³/3 − …, which is the basis of numerical differentiation. Higher differences are formed by repeated application: Δ²f(x) = f(x + 2h) − 2f(x + h) + f(x), and in general Δⁿyₖ = Σ (−1)ʲ C(n, j) yₖ₊ₙ₋ⱼ. A difference table lists the values and successive columns of differences.
The n-th difference of a polynomial of degree n is constant, equal to aₙ n! hⁿ, where aₙ is the leading coefficient, and all higher differences vanish. With h = 1, Δ²(x²) = 2 and Δ³(x³) = 6. This fact lets one detect the degree of a polynomial from a table, check tables for errors (an error ε in one entry spreads through the k-th differences with binomial coefficients and alternating signs, ε, −kε, …), and estimate missing values: if four of five equally spaced values of a cubic are known, Δ⁴y₀ = 0 gives the fifth. For exponential functions, Δaˣ = aˣ(aʰ − 1), so differences of an exponential are again exponential. Factorial polynomials x⁽ⁿ⁾ = x(x − h)(x − 2h)…(x − (n − 1)h) play the role of powers: Δx⁽ⁿ⁾ = nh x⁽ⁿ⁻¹⁾, which makes summation of series by finite differences parallel to integration. These methods appear in actuarial and demographic work, where life tables are tabulated by age and graduated, and in the construction of orthogonal polynomials for equally spaced treatment levels.
Finite differences also sum series, in a way that parallels integration. If Δu(x) = f(x), then Σ from x = a to b of f(x) equals u(b + 1) − u(a) for h = 1, just as an integral equals the difference of an antiderivative. Because Δx⁽ⁿ⁾ = nx⁽ⁿ⁻¹⁾ for factorial polynomials with unit step, the 'anti-difference' of x⁽ⁿ⁾ is x⁽ⁿ⁺¹⁾/(n + 1). To sum ordinary powers, one first writes them in factorial form: x² = x⁽²⁾ + x⁽¹⁾, so Σ from x = 1 to n of x² = [x⁽³⁾/3 + x⁽²⁾/2] from 1 to n + 1, which simplifies to the familiar n(n + 1)(2n + 1)/6. The Euler–Maclaurin formula connects sums and integrals in general, Σf(x) ≈ ∫f(x)dx + [f(a) + f(b)]/2 + correction terms involving derivatives at the ends, and is the basis for refining numerical integration and for approximating sums that appear in likelihoods and in demographic tables.
Interpolation estimates a function between tabulated points by a polynomial passing through them. For equally spaced points, Newton's forward difference formula, with p = (x − x₀)/h, is y = y₀ + pΔy₀ + p(p − 1)Δ²y₀/2! + p(p − 1)(p − 2)Δ³y₀/3! + …, and is used near the beginning of a table. Newton's backward formula, with p = (x − xₙ)/h, uses backward differences at the end of the table: y = yₙ + p∇yₙ + p(p + 1)∇²yₙ/2! + …. Central difference formulae (Gauss, Stirling, Bessel) are used near the middle. For values 1, 3 and 7 at x = 0, 1, 2, Δy₀ = 2 and Δ²y₀ = 2, so y = 1 + 2p + p(p − 1) = p² + p + 1, giving 4.75 at x = 1.5. The error of interpolation by a polynomial of degree n through n + 1 points is f⁽ⁿ⁺¹⁾(ξ)(x − x₀)(x − x₁)…(x − xₙ)/(n + 1)! for some ξ in the range, so accuracy is best within the table and deteriorates in extrapolation.
For unequally spaced points, Lagrange's formula writes the interpolating polynomial as y(x) = Σ yᵢ Lᵢ(x), where Lᵢ(x) = Π over j ≠ i of (x − xⱼ)/(xᵢ − xⱼ) equals 1 at xᵢ and 0 at the other nodes. Through (1, 2), (2, 3) and (4, 11) it gives the quadratic x² − 2x + 3, whose value at x = 3 is 6. Lagrange's form is convenient for theory but must be recomputed when a point is added. Newton's divided difference formula avoids this: the first divided difference is f[x₀, x₁] = (f(x₁) − f(x₀))/(x₁ − x₀), higher ones are defined recursively, and y = f(x₀) + (x − x₀)f[x₀, x₁] + (x − x₀)(x − x₁)f[x₀, x₁, x₂] + …. Divided differences are symmetric in their arguments, and for equally spaced points reduce to ordinary differences divided by k! hᵏ. Inverse interpolation, finding x for a given y, interchanges the roles of x and y and is used to find, for example, the dose giving a stated response.
High-degree polynomial interpolation through many equally spaced points can oscillate wildly near the ends of the interval (Runge's phenomenon), so in practice low-degree piecewise polynomials, especially cubic splines, are preferred; splines also underlie smoothing in generalised additive models. Interpolation assumes the function is smooth; noisy data call for regression or smoothing instead, because forcing a curve through every point fits the noise.
Numerical differentiation estimates derivatives from tabulated values by differentiating an interpolating polynomial. From Newton's forward formula, at x = x₀, dy/dx ≈ (1/h)[Δy₀ − Δ²y₀/2 + Δ³y₀/3 − …] and d²y/dx² ≈ (1/h²)[Δ²y₀ − Δ³y₀ + (11/12)Δ⁴y₀ − …], which follow from hD = ln(1 + Δ). From the backward formula at the end of the table, dy/dx ≈ (1/h)[∇yₙ + ∇²yₙ/2 + ∇³yₙ/3 + …]. Simple finite-difference approximations are used most: the forward difference f′(x) ≈ [f(x + h) − f(x)]/h has error proportional to h; the central difference f′(x) ≈ [f(x + h) − f(x − h)]/(2h) has error proportional to h² (−h²f‴/6); and the second derivative f″(x) ≈ [f(x + h) − 2f(x) + f(x − h)]/h² also has error of order h². For f(x) = x² at x = 1 with h = 0.1, the forward difference gives 2.1 and the central difference gives exactly 2, since central differences are exact for quadratics.
Numerical differentiation is ill-conditioned: as h shrinks, truncation error falls but rounding error grows, because a small difference of nearly equal numbers is divided by a small h. The total error is smallest at an intermediate h, roughly the square root of machine precision for forward differences and its cube root for central differences, scaled to the size of x. Richardson extrapolation combines estimates at h and h/2 to cancel the leading error term: for the central difference, D = [4D(h/2) − D(h)]/3 has error of order h⁴. In statistics, numerical derivatives are used to approximate the Hessian of a log-likelihood when analytic second derivatives are not available, giving standard errors of estimates; poor step choices produce visibly unstable standard errors.
Numerical integration (quadrature) approximates ∫ₐᵇ f(x)dx from values at points x₀ = a, x₁, …, xₙ = b spaced h = (b − a)/n apart. The trapezoidal rule joins points by straight lines: ∫ ≈ (h/2)[y₀ + 2(y₁ + … + yₙ₋₁) + yₙ], with error −(b − a)h²f″(ξ)/12, so it is exact for linear functions and its error falls by a factor of four when h is halved. Simpson's one-third rule fits parabolas through successive pairs of intervals and needs an even number of intervals: ∫ ≈ (h/3)[y₀ + 4(y₁ + y₃ + …) + 2(y₂ + y₄ + …) + yₙ], with error −(b − a)h⁴f⁽⁴⁾(ξ)/180; it is exact for cubics, and halving h reduces the error sixteen-fold. Simpson's three-eighths rule, needing a multiple of three intervals, is (3h/8)[y₀ + 3y₁ + 3y₂ + 2y₃ + 3y₄ + 3y₅ + 2y₆ + … + yₙ]. Weddle's rule uses six intervals and is more accurate still for smooth functions.
For ∫₀¹ x² dx with h = 0.5, the trapezoidal rule gives (0.25)[0 + 2(0.25) + 1] = 0.375, while Simpson's rule gives (0.5/3)[0 + 4(0.25) + 1] = 1/3, the exact value. For ∫₀¹ dx/(1 + x) = ln 2 ≈ 0.693147 with h = 0.25, the ordinates are 1, 0.8, 0.6667, 0.5714 and 0.5; the trapezoidal rule gives 0.69702 and Simpson's rule 0.69325, errors of about 0.0039 and 0.0001. Romberg integration applies Richardson extrapolation to trapezoidal estimates with successively halved h, and Gaussian quadrature chooses both the points and weights optimally: the two-point Gauss–Legendre rule on [−1, 1], ∫ ≈ f(−1/√3) + f(1/√3), is exact for cubics with only two function values.
Numerical integration computes probabilities and quantiles of distributions without closed-form distribution functions, such as the normal, t and non-central distributions; areas under disease-progress curves, a standard summary in plant pathology computed by the trapezoidal rule from repeated severity assessments; areas under growth and dose-response curves; and marginal likelihoods in generalised linear mixed models, where adaptive Gauss–Hermite quadrature integrates out random effects. Where the dimension is high, quadrature becomes impractical, and Monte Carlo integration, whose error falls as 1/√n regardless of dimension, takes over.
A difference equation relates values of a sequence at successive steps, yₙ₊ₖ = f(yₙ, …, yₙ₊ₖ₋₁, n); its order is k. The first-order linear equation yₙ₊₁ = ayₙ + b has the equilibrium y* = b/(1 − a) when a ≠ 1, and the solution yₙ = y* + (y₀ − y*)aⁿ; it converges to equilibrium if |a| < 1, monotonically if 0 < a < 1 and with damped oscillation if −1 < a < 0. For yₙ₊₁ = 0.5yₙ + 10, the equilibrium is 20. A linear homogeneous equation with constant coefficients, such as yₙ₊₂ − 5yₙ₊₁ + 6yₙ = 0, is solved by trying yₙ = mⁿ, giving the characteristic equation m² − 5m + 6 = 0 with roots 2 and 3, so yₙ = A2ⁿ + B3ⁿ. A repeated root m gives (A + Bn)mⁿ, and complex roots re^(±iθ) give rⁿ(A cos nθ + B sin nθ). Non-homogeneous equations add a particular solution found by undetermined coefficients. The Fibonacci recurrence yₙ₊₂ = yₙ₊₁ + yₙ has characteristic roots (1 ± √5)/2, and its terms grow like the golden ratio 1.618 to the n.
Difference equations describe systems observed at discrete times. The cobweb model of agricultural markets, in which supply depends on last season's price and demand on this season's, gives a first-order equation whose price oscillations converge or explode depending on the ratio of supply to demand slopes. Discrete population models range from the geometric Nₜ₊₁ = λNₜ to the logistic map Nₜ₊₁ = rNₜ(1 − Nₜ), whose behaviour passes from stable equilibrium through cycles to chaos as r increases. The Leslie matrix model, nₜ₊₁ = Lnₜ, is a vector difference equation whose dominant eigenvalue gives the long-run growth rate. In time series, an autoregressive process yₜ = φ₁yₜ₋₁ + φ₂yₜ₋₂ + εₜ is a stochastic difference equation, and it is stationary exactly when the roots of the characteristic equation m² − φ₁m − φ₂ = 0 lie inside the unit circle. Interest, annuity and loan schedules in farm finance are also first-order difference equations.
Many equations f(x) = 0, such as likelihood equations, cannot be solved in closed form. The bisection method starts from an interval [a, b] with f(a) and f(b) of opposite sign, so that a root lies inside by the intermediate value theorem, and repeatedly halves the interval, keeping the half where the sign changes. After n steps the root is located within (b − a)/2ⁿ; to reach a tolerance of 0.001 from an interval of length 1 requires 2ⁿ > 1,000, so n = 10 steps. Bisection always converges but slowly (linear convergence, gaining one binary digit per step). The regula falsi (false position) method replaces the midpoint with the x-intercept of the chord between the end points; it keeps the bracket and usually converges faster, though one end point can stagnate.
The Newton–Raphson method uses the tangent at the current estimate: xₙ₊₁ = xₙ − f(xₙ)/f′(xₙ). Near a simple root it converges quadratically, roughly doubling the number of correct digits each step. For √2, solving x² − 2 = 0 from x₀ = 1 gives 1.5, 1.41667 and 1.414216. For x³ − x − 1 = 0 from x₀ = 1.5, the iterates are 1.3478 and 1.3252, approaching the root 1.3247. The method can fail when f′ is near zero, when the start is far from the root, or at multiple roots, where convergence becomes linear. The secant method replaces f′ with the slope through the last two iterates, avoiding derivatives at the cost of a slightly lower order of convergence, about 1.618. Fixed-point iteration rewrites the equation as x = g(x) and iterates xₙ₊₁ = g(xₙ); it converges near a fixed point when |g′(x)| < 1 there.
Maximum likelihood estimation solves the score equation U(θ) = ∂ℓ/∂θ = 0. Newton–Raphson iterates θₙ₊₁ = θₙ − U(θₙ)/ℓ″(θₙ), or in several parameters θₙ₊₁ = θₙ − H⁻¹U, where H is the Hessian of the log-likelihood. Fisher scoring replaces the observed Hessian by its expectation, the negative Fisher information, which is more stable and, for generalised linear models with canonical links, gives the iteratively reweighted least squares algorithm used by standard software. Good starting values, such as method-of-moments estimates, and step-halving when the likelihood decreases, make these algorithms reliable; failure to converge often signals a problem with the model or data, such as complete separation in logistic regression.
Gaussian elimination reduces a system Ax = b to upper triangular form by row operations and then solves by back substitution. It requires about 2n³/3 arithmetic operations for an n × n system, so doubling n multiplies the work by eight. Partial pivoting, choosing as pivot the largest available entry in absolute value in the current column, controls the growth of rounding errors and is standard. Elimination is equivalent to factorising A = LU into lower and upper triangular matrices (with row permutations if pivoting is used); once the factorisation is done, systems with new right-hand sides cost only about 2n² operations each. Gauss–Jordan elimination continues to reduced form and gives the inverse, at greater cost, but explicit inverses are rarely needed: solving the system directly is cheaper and more accurate. For symmetric positive definite matrices, such as X′X and covariance matrices, the Cholesky factorisation A = LL′ takes half the work of LU and is numerically stable without pivoting. The condition number of A bounds how much relative errors in b or A can be magnified in x; ill-conditioned normal equations are why regression software uses QR factorisation of X instead of forming X′X.
For large sparse systems, iterative methods can be much cheaper than elimination. The Jacobi method solves each equation for its diagonal unknown using the previous iterate's values for the others: xᵢ⁽ᵏ⁺¹⁾ = (bᵢ − Σⱼ≠ᵢ aᵢⱼxⱼ⁽ᵏ⁾)/aᵢᵢ. The Gauss–Seidel method uses each new value as soon as it is computed, and usually converges about twice as fast. Both converge for any starting vector if A is strictly diagonally dominant (each diagonal entry exceeds in absolute value the sum of the others in its row), and Gauss–Seidel also converges for symmetric positive definite A. For 10x + y = 11 and x + 10y = 11, starting from zero, the first Jacobi iterate is (1.1, 1.1), already close to the solution (1, 1). Successive over-relaxation accelerates Gauss–Seidel with a relaxation factor between 1 and 2, and conjugate-gradient methods are standard for large symmetric positive definite systems.
In animal and plant breeding, the mixed-model equations for predicting breeding values can involve millions of unknowns, one for each animal or tree in a pedigree, with a very sparse coefficient matrix; they have been solved for decades by Gauss–Seidel and preconditioned conjugate-gradient iteration rather than by inversion.
The condition number κ(A) = ‖A‖‖A⁻¹‖, which for the usual two-norm equals the ratio of the largest to the smallest singular value, measures how sensitive the solution of Ax = b is to small changes in A or b: relative errors in the data can be magnified by a factor of up to κ in the solution. A rule of thumb is that one loses about log₁₀κ significant digits; with κ = 10⁸ and double-precision arithmetic of about 16 digits, only about eight digits of the solution can be trusted. Ill-conditioning is a property of the problem, not of the algorithm, so no method can rescue a badly conditioned system, but unstable algorithms, such as elimination without pivoting, can make matters worse. In regression, near-collinear explanatory variables make X ill-conditioned, so coefficient estimates become unstable and sensitive to rounding and to small changes in the data; centring and scaling variables, removing redundant ones or using ridge regression reduce the problem. The residual of a computed solution, b − Ax, can be small even when the solution is inaccurate, so a small residual alone does not prove accuracy for ill-conditioned systems.
For the initial value problem y′ = f(x, y), y(x₀) = y₀, numerical methods compute approximations yₙ at points xₙ = x₀ + nh. Euler's method follows the tangent: yₙ₊₁ = yₙ + hf(xₙ, yₙ). Its local error per step is of order h² and its global error of order h, so it is a first-order method: halving h roughly halves the error. For y′ = y, y(0) = 1, one Euler step of h = 0.1 gives 1.1, against the exact e^0.1 ≈ 1.105171. The modified Euler (Heun) method averages the slopes at the start and at the Euler-predicted end: yₙ₊₁ = yₙ + (h/2)[f(xₙ, yₙ) + f(xₙ₊₁, yₙ + hf(xₙ, yₙ))]; for the same problem it gives 1 + 0.05(1 + 1.1) = 1.105, a second-order method. Taylor series methods use higher derivatives of y computed from the equation, and are accurate but need those derivatives.
Runge–Kutta methods achieve high order using only evaluations of f. The classical fourth-order method computes k₁ = f(xₙ, yₙ), k₂ = f(xₙ + h/2, yₙ + hk₁/2), k₃ = f(xₙ + h/2, yₙ + hk₂/2) and k₄ = f(xₙ + h, yₙ + hk₃), and sets yₙ₊₁ = yₙ + (h/6)(k₁ + 2k₂ + 2k₃ + k₄). Its global error is of order h⁴. For y′ = y with h = 0.1, k₁ = 1, k₂ = 1.05, k₃ = 1.0525 and k₄ = 1.10525, giving y₁ = 1.1051708, correct to seven significant figures. Adaptive methods pair Runge–Kutta formulae of different orders to estimate the local error and adjust h automatically, and are the default in scientific software. Multistep methods (Adams–Bashforth and Adams–Moulton predictor–corrector) reuse earlier values to save evaluations.
A method is stable for a problem if errors do not grow from step to step. For the test equation y′ = −λy with λ > 0, Euler's method is stable only if h < 2/λ; with a larger step the numerical solution oscillates and explodes although the true solution decays. Stiff systems, with processes on very different time scales, such as fast and slow pools in soil carbon models, force explicit methods to tiny steps, and implicit methods, such as backward Euler, are used instead. Crop growth simulation models, pest and disease epidemic models and pharmacokinetic models are systems of ODEs solved numerically; their outputs depend on the time step, so modellers check that results do not change materially when the step is reduced.
Many agricultural and ecological models are differential equations whose parameters must be estimated from observations, such as growth rates and carrying capacities from repeated weighings, decay rates of pesticide residues from soil samples, or transmission rates of a plant disease from field surveys. When an analytical solution exists, as for the logistic or monomolecular models, the solution curve is fitted directly by non-linear least squares or maximum likelihood. When it does not, each evaluation of the fitting criterion requires solving the equations numerically for the current parameter values, so the ODE solver sits inside the optimisation loop. The solver's tolerance must then be tight enough that numerical error does not distort the estimated parameters or their standard errors, and sensitivities of the solution to the parameters, needed for gradient-based optimisers, can be computed by solving augmented sensitivity equations alongside the model. Identifiability matters: some parameter combinations may not be estimable from the available measurements, which shows up as very flat likelihoods and huge standard errors.
Coordinate geometry describes curves by equations. The distance between (x₁, y₁) and (x₂, y₂) is √[(x₂ − x₁)² + (y₂ − y₁)²]; a line has slope m = (y₂ − y₁)/(x₂ − x₁) and equation y − y₁ = m(x − x₁); two lines are perpendicular if the product of their slopes is −1; and the perpendicular distance from (x₁, y₁) to ax + by + c = 0 is |ax₁ + by₁ + c|/√(a² + b²). The circle x² + y² + 2gx + 2fy + c = 0 has centre (−g, −f) and radius √(g² + f² − c); for x² + y² − 4x + 6y − 3 = 0 the centre is (2, −3) and the radius √(4 + 9 + 3) = 4. The general second-degree equation ax² + 2hxy + by² + 2gx + 2fy + c = 0 represents a conic. When it is non-degenerate, it is an ellipse if h² < ab (a circle if, in addition, a = b and h = 0), a parabola if h² = ab, and a hyperbola if h² > ab; the cross term 2hxy indicates axes rotated relative to the coordinate axes, removed by rotating through an angle θ with tan 2θ = 2h/(a − b).
A conic is the locus of points whose distance from a fixed point (focus) is e times their distance from a fixed line (directrix); e is the eccentricity: e = 0 for a circle, 0 < e < 1 for an ellipse, e = 1 for a parabola and e > 1 for a hyperbola. The parabola y² = 4ax has focus (a, 0), directrix x = −a and latus rectum of length 4a; y² = 12x has focus (3, 0). The ellipse x²/a² + y²/b² = 1 with a > b has eccentricity e = √(1 − b²/a²), foci (±ae, 0) and directrices x = ±a/e; for x²/25 + y²/9 = 1, e = 4/5 and the foci are (±4, 0). The hyperbola x²/a² − y²/b² = 1 has e = √(1 + b²/a²), so x²/16 − y²/9 = 1 has e = 5/4, and asymptotes y = ±(b/a)x; a rectangular hyperbola, a = b, has e = √2, and xy = c² is its form with the asymptotes as axes.
Conics appear wherever quadratic forms do. The contours of a bivariate normal density are ellipses, (x − μ)′Σ⁻¹(x − μ) = constant, centred at the mean and with axes along the eigenvectors of Σ; their size for a given probability content comes from the chi-square distribution with two degrees of freedom. Joint confidence regions for two regression coefficients are ellipses, and so are prediction ellipses used to detect outliers in two variables. A quadratic response surface y = β₀ + β₁x₁ + β₂x₂ + β₁₁x₁² + β₂₂x₂² + β₁₂x₁x₂, fitted in fertiliser or irrigation experiments, has elliptical contours around a maximum or minimum when the quadratic part is definite and hyperbolic contours around a saddle point when it is indefinite, which canonical analysis identifies. Rectangular hyperbolas describe Michaelis–Menten and other saturating responses, and parabolas the quadratic yield–dose curves common in agronomy.
A student uses this subject to solve the differential equations behind growth curves, to build difference tables and interpolate values, to integrate response curves by Simpson's rule, and to understand how software finds maximum likelihood estimates by Newton–Raphson. In examinations, the key skills are recognising the type of a differential equation and its integrating factor, writing the complementary function from the roots of the auxiliary equation, applying Newton's and Lagrange's formulae correctly, choosing the right integration rule for the number of intervals, and carrying out two or three Newton–Raphson or Runge–Kutta steps accurately.
A statistician five years into practice uses these methods inside software and must recognise when they go wrong: a non-converging likelihood, a numerical Hessian that gives implausible standard errors, an interpolated reference value outside the table, or a growth model whose fitted curve depends on the solver's step size. They compute areas under disease-progress curves, fit logistic and Gompertz curves to crop and livestock growth data, and interpolate missing values in weather or price series with appropriate caution.
A senior statistician with twenty-five years of experience chooses the numerical approach for large institutional analyses, such as national breeding evaluations solved iteratively, crop simulation models run across regions, and population projections, and sets standards for checking convergence, accuracy and reproducibility. They advise scientists on when a mechanistic differential-equation model is justified and when a statistical curve fit is enough, and review the numerical methods described in reports and papers. They also decide when an analysis needs a specialist numerical solver rather than default settings.
A researcher with a PhD and long experience develops numerical methods tailored to statistical problems: algorithms for fitting non-linear mixed models, quadrature and Laplace approximations for complex likelihoods, solvers for stiff ecological and physiological models with parameters estimated from data, and efficient iterative methods for very large sparse systems in genomics. Such researchers teach numerical analysis for statisticians, assess the numerical correctness of new software, and examine candidates on both the mathematics and its computation.
This manual covers classical analytical and numerical methods to the level of master's admission and doctoral computer-application examinations. It does not cover partial differential equations, boundary value problems, eigenvalue algorithms, optimisation algorithms beyond Newton–Raphson, or rigorous error analysis of floating-point arithmetic. Operation counts and convergence rates are stated in their standard asymptotic form; real performance depends on implementation.
Original practice questions written to the syllabus and pattern of the exams named in each question's examTags. They are not past papers and do not reproduce any board's questions.
These notes are free to read here, with no account.
What the Laisha app adds:
Free for ten days, then a subscription. We would rather write that here than have you find it later.