Taylor Series, Linear Systems by Elimination and Iteration, Linear Regression, and the One-Dimensional Diffusivity and Laplace Equations
1. Taylor series of a function of one variable
If f has derivatives of every order at x = a, its Taylor series about a is f(x) = f(a) + f′(a)(x − a) + f″(a)(x − a)²/2! + f‴(a)(x − a)³/3! + …, and about a = 0 it is called the Maclaurin series. Truncating after the term in (x − a)ⁿ leaves a remainder that, in Lagrange’s form, is Rₙ = f⁽ⁿ⁺¹⁾(ξ)(x − a)ⁿ⁺¹/(n + 1)! for some ξ between a and x — so the error of a truncated series is controlled by the first neglected term. Five Maclaurin series are worth knowing cold: eˣ = 1 + x + x²/2! + x³/3! + … (all x); sin x = x − x³/3! + x⁵/5! − … and cos x = 1 − x²/2! + x⁴/4! − … (all x); ln(1 + x) = x − x²/2 + x³/3 − … (−1 < x ≤ 1); and 1/(1 − x) = 1 + x + x² + … (|x| < 1).
Taylor’s series is the source of every finite-difference formula. Writing f(x + h) and f(x − h) about x and subtracting gives the central difference f′(x) ≈ [f(x + h) − f(x − h)]/(2h), with truncation error −h²f‴(x)/6, of order h²; keeping only f(x + h) gives the forward difference [f(x + h) − f(x)]/h, of order h. Adding the two expansions gives the second derivative, f″(x) ≈ [f(x + h) − 2f(x) + f(x − h)]/h², also of order h². That last formula is exactly how a reservoir simulator discretises the ∂²p/∂x² of the diffusivity equation on a grid, which is why the numerical-methods and reservoir-simulation sub-items of this paper are the same piece of mathematics.
| Formula | Approximates | Truncation error |
|---|---|---|
| [f(x + h) − f(x)]/h | f′(x), forward | O(h) |
| [f(x) − f(x − h)]/h | f′(x), backward | O(h) |
| [f(x + h) − f(x − h)]/(2h) | f′(x), central | O(h²) |
| [f(x + h) − 2f(x) + f(x − h)]/h² | f″(x), central | O(h²) |
2. Numerical solution of linear algebraic equations
A system Ax = b can be solved directly or iteratively. Gauss elimination is direct: subtract multiples of each pivot row from the rows below it until A is upper triangular, then back-substitute from the last unknown upward. It costs about n³/3 multiplications for n equations, and it fails outright if a pivot is zero — and loses accuracy if a pivot is merely small — so practical codes use partial pivoting, swapping in the row with the largest available pivot before each elimination step. A tridiagonal system, which is what a one-dimensional finite-difference grid produces, is solved by the Thomas algorithm, Gauss elimination that touches only the three diagonals and costs of order n.
Jacobi iteration rewrites each equation for its diagonal unknown and computes every new value from the previous iterate: xᵢ⁽ᵏ⁺¹⁾ = [bᵢ − Σⱼ≠ᵢ aᵢⱼxⱼ⁽ᵏ⁾]/aᵢᵢ. Gauss-Seidel is the same formula except that each new value is used as soon as it is computed, so within one sweep x₂ already uses the new x₁; it usually converges about twice as fast as Jacobi. Both are guaranteed to converge from any starting guess when A is strictly diagonally dominant (|aᵢᵢ| > Σⱼ≠ᵢ |aᵢⱼ| in every row), and Gauss-Seidel also converges for any symmetric positive-definite A. Reordering the equations so that the largest coefficients sit on the diagonal is often all that a non-converging system needs.
3. Linear regression analysis
Given n pairs (xᵢ, yᵢ), the least-squares line y = a + bx is the one that minimises the sum of squared vertical deviations Σ(yᵢ − a − bxᵢ)². Setting the two partial derivatives to zero gives the normal equations Σy = na + bΣx and Σxy = aΣx + bΣx², whose solution is b = Sxy/Sxx and a = ȳ − b x̄, where Sxy = Σ(xᵢ − x̄)(yᵢ − ȳ) and Sxx = Σ(xᵢ − x̄)². The fitted line always passes through the point of means (x̄, ȳ). The correlation coefficient r = Sxy/√(Sxx·Syy) lies between −1 and 1 and has the sign of the slope; its square, the coefficient of determination R², is the fraction of the variance of y that the line explains. The regression of y on x and that of x on y are different lines unless |r| = 1, and they cross at (x̄, ȳ).
Most petroleum data are made straight before they are fitted. An exponential production decline q = qᵢe^(−Dt) becomes the line ln q = ln qᵢ − Dt, so the decline rate is minus the slope of ln q against t. A pressure build-up plotted against the logarithm of the Horner time ratio is a line whose slope gives permeability, and a drawdown plotted against log t is a line whose slope does the same. In each case the regression is ordinary least squares on transformed variables, and the physical constant is read from the slope, not the intercept.
| Quantity | Value |
|---|---|
| x̄, ȳ | 2.5, 4.75 |
| Sxx = Σ(x − x̄)² | 5 |
| Sxy = Σ(x − x̄)(y − ȳ) | 9.5 |
| Syy = Σ(y − ȳ)² | 18.75 |
| Slope b = Sxy/Sxx; intercept a = ȳ − b x̄ | 1.9; 0 |
| Correlation coefficient r | 0.981 |
4. Solving the one-dimensional diffusivity equation and the Laplace equation
Slightly compressible single-phase flow through a porous medium obeys the diffusivity equation. In one linear dimension it is ∂²p/∂x² = (1/η)∂p/∂t, with hydraulic diffusivity η = k/(φμcₜ); in radial form it is (1/r)∂/∂r(r ∂p/∂r) = (1/η)∂p/∂t. Mathematically it is the heat equation: parabolic, diffusive, smoothing its initial data. Two solution methods matter. Separation of variables for a bounded domain 0 < x < L with p held at zero at both ends gives p(x, t) = Σ Bₙ sin(nπx/L) exp(−n²π²ηt/L²), with the Bₙ the Fourier sine coefficients of the initial profile; the higher modes die fastest, so at late time only the first survives. A similarity solution for a semi-infinite domain whose face at x = 0 is suddenly changed from pᵢ to p₀ gives (p − p₀)/(pᵢ − p₀) = erf(x/(2√(ηt))): the disturbance travels as √(ηt), which is why a pressure pulse needs four times as long to reach twice the distance.
When nothing changes with time the right-hand side vanishes and the diffusivity equation becomes the Laplace equation, the steady state. In one linear dimension d²p/dx² = 0, so p is linear in x between its boundary values. In radial flow d/dr(r dp/dr) = 0 integrates to p = C₁ ln r + C₂, and with p = p_w at the wellbore radius r_w and p = p_e at the outer radius r_e, p(r) = p_w + (p_e − p_w) ln(r/r_w)/ln(r_e/r_w). The pressure is logarithmic in radius, so most of the drawdown is spent close to the well — the reason a damaged zone a few feet thick (a skin) matters so much. Substituting this profile into Darcy’s law gives the steady-state radial inflow equation of the Reservoir Engineering chapter.
- Initial-value versus boundary-value. The diffusivity equation needs an initial condition and two boundary conditions; the Laplace equation needs boundary conditions only.
- Constant-rate inner boundary. A well producing at constant rate fixes the gradient r ∂p/∂r at r_w (a Neumann condition), not the pressure; the line-source solution of the Oil and Gas Well Testing chapter is the answer to exactly that problem.
- Diffusivity sets the time scale. A time scale L²/η: higher permeability shortens it, higher viscosity, porosity or compressibility lengthens it.
Key takeaways
- A truncated Taylor series errs by roughly its first neglected term; central differences are O(h²), forward and backward differences O(h).
- Gauss elimination is direct and needs pivoting; Jacobi and Gauss-Seidel are iterative and converge for strictly diagonally dominant systems, Gauss-Seidel faster because it uses new values at once.
- Least squares: b = Sxy/Sxx, a = ȳ − b x̄, the line passes through the means, and r = Sxy/√(Sxx·Syy).
- Linearise first: ln q against t for exponential decline, pressure against log time for transient tests; the physics is in the slope.
- The diffusivity equation has η = k/(φμcₜ) and spreads as √(ηt); at steady state it becomes Laplace’s equation, linear in x and logarithmic in r.
Practice questions (13)
Attempt each one before opening the answer. Every explanation names the tempting wrong option as well as the right one, because that is where marks are lost.
Using the first three terms of the Maclaurin series of eˣ, the approximate value of e^0.2 is ______ (to two decimal places).
Numerical answer — type the value.
Show answer
Answer: 1.22
e^x ≈ 1 + x + x²/2 = 1 + 0.2 + 0.04/2 = 1 + 0.2 + 0.02 = 1.22. The true value is 1.2214; the first neglected term x³/6 = 0.0013 accounts for almost all of the difference, as the Lagrange remainder says it should.The Maclaurin series of ln(1 + x), valid for −1 < x ≤ 1, is
Show answer
Answer: A — x − x²/2 + x³/3 − x⁴/4 + …
Differentiate: d/dx ln(1 + x) = 1/(1 + x) = 1 − x + x² − …, and integrating term by term from 0 gives x − x²/2 + x³/3 − …. The all-plus series is −ln(1 − x), the third is eˣ and the fourth is sin x.The derivative of f(x) = x³ at x = 1 is estimated by the central difference [f(x + h) − f(x − h)]/(2h) with h = 0.1. The estimate is ______ (to two decimal places).
Numerical answer — type the value.
Show answer
Answer: 3.01
f(1.1) = 1.331 and f(0.9) = 0.729, so the estimate is (1.331 − 0.729)/0.2 = 0.602/0.2 = 3.01. Check by the error formula: the exact derivative is 3 and the truncation error is h²f‴/6 = 0.01 × 6/6 = 0.01, giving 3.01 again.If the step size h of a central-difference estimate of f′(x) is halved, its truncation error becomes approximately
Show answer
Answer: A — one quarter of what it was
The central difference has truncation error −h²f‴/6, of order h², so halving h divides the error by 2² = 4. One half would be the answer for a forward or backward difference, which is only first order.The system 4x + y = 9, x + 3y = 7 is solved by Gauss-Seidel iteration starting from x = 0, y = 0, updating x first. The value of y after the first iteration is ______ (to two decimal places).
Numerical answer — type the value.
Show answer
Answer: 1.58
x₁ = (9 − y₀)/4 = 9/4 = 2.25. Gauss-Seidel uses this new x at once: y₁ = (7 − x₁)/3 = 4.75/3 = 1.583, which is 1.58. Jacobi would have used x₀ = 0 and given 2.33, which is the distractor to avoid.Which of the following guarantee that Gauss-Seidel iteration converges for Ax = b from any starting vector? (More than one option may be correct.)
Show answer
Answer: A — A is strictly diagonally dominant; B — A is symmetric and positive definite
Both strict diagonal dominance and symmetric positive-definiteness are sufficient conditions for Gauss-Seidel. Non-singularity guarantees that a solution exists, not that the iteration finds it, and a zero diagonal entry makes the update formula divide by zero.Partial pivoting is used in Gauss elimination mainly to
Show answer
Answer: A — avoid division by zero or by very small pivots and so limit round-off error
Swapping the row with the largest available pivot into place keeps the multipliers at most 1 in magnitude, which stops round-off from being amplified and removes the zero-pivot failure. It does not change the order-n³ cost, the symmetry of A, or the direct nature of the method.For the data x = 1, 2, 3, 4 and y = 2, 4, 5, 8, the slope of the least-squares regression line of y on x is ______ (to one decimal place).
Numerical answer — type the value.
Show answer
Answer: 1.9
x̄ = 2.5 and ȳ = 19/4 = 4.75. Sxy = (−1.5)(−2.75) + (−0.5)(−0.75) + (0.5)(0.25) + (1.5)(3.25) = 4.125 + 0.375 + 0.125 + 4.875 = 9.5 and Sxx = 2.25 + 0.25 + 0.25 + 2.25 = 5, so b = 9.5/5 = 1.9. Check with the normal equations: Σx = 10, Σx² = 30, Σy = 19, Σxy = 2 + 8 + 15 + 32 = 57; b = (4·57 − 10·19)/(4·30 − 100) = 38/20 = 1.9.For the same data (x = 1, 2, 3, 4; y = 2, 4, 5, 8), the correlation coefficient r is closest to
Show answer
Answer: A — 0.98
Syy = 7.5625 + 0.5625 + 0.0625 + 10.5625 = 18.75, so r = 9.5/√(5 × 18.75) = 9.5/√93.75 = 9.5/9.682 = 0.981. The sign is positive because the slope is positive; R² = 0.963, so the line explains about 96% of the variance of y.Production rate data believed to follow q = qᵢe^(−Dt) are best checked for a straight line on a plot of
Show answer
Answer: A — ln q against t
Taking logarithms gives ln q = ln qᵢ − Dt, which is linear in t with slope −D and intercept ln qᵢ. A plot of q against t is curved, and 1/q against t is the straight line of a harmonic decline instead.Steady radial flow obeys the Laplace equation. With p = 1000 psi at r_w = 0.5 ft and p = 3000 psi at r_e = 500 ft, the pressure at r = 5 ft is ______ psi (to the nearest whole number).
Numerical answer — type the value.
Show answer
Answer: 1667
p(r) = p_w + (p_e − p_w) ln(r/r_w)/ln(r_e/r_w). Here r/r_w = 10 and r_e/r_w = 1000, so the fraction is ln 10/ln 1000 = 1/3 and p = 1000 + 2000/3 = 1666.7, which rounds to 1667 psi. A third of the whole drawdown is recovered within 4.5 ft of the well — the logarithmic profile at work.For the one-dimensional diffusivity equation ∂²p/∂x² = (1/η)∂p/∂t with η = k/(φμcₜ), which statements are correct? (More than one option may be correct.)
Show answer
Answer: A — It is a parabolic partial differential equation; B — A pressure disturbance penetrates a distance proportional to √(ηt); D — At steady state it reduces to d²p/dx² = 0, so p is linear in x
With one time derivative against two space derivatives the equation is parabolic (B² − 4AC = 0). The similarity variable x/(2√(ηt)) makes penetration scale as √(ηt). Setting ∂p/∂t = 0 leaves d²p/dx² = 0, a linear profile. Viscosity sits in the denominator of η, so doubling μ halves the diffusivity.In the separation-of-variables solution of ∂p/∂t = η ∂²p/∂x² on 0 < x < L with p = 0 at both ends, the n-th Fourier mode decays with time as
Show answer
Answer: A — exp(−n²π²ηt/L²)
Substituting X(x)T(t) gives X″ + λX = 0 with X(0) = X(L) = 0, so λₙ = (nπ/L)², and then T′ = −ηλₙT, so T = exp(−n²π²ηt/L²). A cosine in time would belong to the wave equation, and η must multiply t, not divide it, because a more diffusive medium relaxes faster.