Numerical Analysis

Newton divided difference interpolation

Newton's divided difference form writes the interpolating polynomial through nodes x_0,\dots,x_n as P(x)=f[x_0]+f[x_0,x_1](x-x_0)+\cdots+f[x_0,\dots,x_n]\prod_{j=0}^{n-1}(x-x_j). The divided differences are defined recursively by f[x_i]=f(x_i) and f[x_i,\dots,x_{i+k}]=\dfrac{f[x_{i+1},\dots,x_{i+k}]-f[x_i,\dots,x_{i+k-1}]}{x_{i+k}-x_i}. Because the coefficients are the top diagonal of the difference table, the Newton form lets you read the coefficients of the polynomial directly, and it is convenient when nodes are added one at a time.

Key facts to remember:

• f[x_0,\dots,x_k] is symmetric in its nodes, so the order of the nodes does not change its value.
• If f is a polynomial of degree n, then the divided difference of order n+1 is zero, and the divided difference of order n equals the leading coefficient.
• For a polynomial q of degree at most 4 with leading coefficient c, one has q^{(4)}(x)=24c, and c=f[x_0,\dots,x_4]. Therefore the fourth derivative of an interpolant is 4! times the fourth divided difference.
• The interpolation error is f(x)-P(x)=\dfrac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{j}(x-x_j) for smooth f. When f is not smooth, compute the error function on each subinterval where it is smooth, find its critical points there, and compare the values at those points with the endpoint values to get the supremum.

Example: for nodes 0,1,2 with values 1,3,7, we get f[0]=1, f[0,1]=2 and f[0,1,2]=1, so P(x)=1+2x+x(x-1)=x^2+x+1. Checking P(2)=7 confirms the result.

Quadrature rule exactness

A quadrature rule \int_a^b f(x)\,dx\approx\sum_i w_i f(x_i) has degree of exactness m if it integrates every polynomial of degree at most m exactly. Since the rule is linear, it suffices to test the monomials 1,x,\dots,x^m. Each test gives one linear equation in the unknown weights (and sometimes nodes), so exactness conditions are used to solve for unknown coefficients.

• Trapezoidal rule — \frac{b-a}{2}[f(a)+f(b)], exact for polynomials of degree at most 1. Its composite form with step h is \frac{h}{2}[f_0+2f_1+\cdots+2f_{n-1}+f_n]. For a single panel of width h, the rule is \frac{h}{2}[f(a)+f(b)].
• Simpson's rule — exact for cubic polynomials, because the symmetric placement of the nodes cancels the odd-degree error term.
• Piecewise behavior — a composite rule is exact for a piecewise polynomial only if every kink or break of f lies at a node. A kink inside a panel destroys exactness of that panel's rule.
• Derivative data — rules may use values of f' at the nodes. Testing the monomials then requires the derivatives of the test functions, for example \frac{d}{dx}x^k=kx^{k-1}, and this yields additional equations for the weights.

Example: to find the rule \int_0^1 f\,dx\approx w_0f(0)+w_1f(1) exact for degree 1, testing f=1 gives w_0+w_1=1, and testing f=x gives w_1=\frac12. Hence w_0=w_1=\frac12, which is the trapezoidal rule.

Jacobi and Gauss-Seidel iteration

Write A=D+L+U, where D is the diagonal, L the strictly lower part and U the strictly upper part. Jacobi iteration is x^{(k+1)}=D^{-1}\big(b-(L+U)x^{(k)}\big) with iteration matrix T_J=-D^{-1}(L+U). Gauss-Seidel iteration is x^{(k+1)}=(D+L)^{-1}\big(b-Ux^{(k)}\big), which uses the newest components as soon as they are computed. In a Jacobi step, every component uses only the old values, so each new component is the old right-hand side minus off-diagonal contributions, divided by the diagonal entry.

Convergence for every initial vector holds if and only if the spectral radius of the iteration matrix is less than 1. A convenient sufficient condition is strict diagonal dominance, namely |a_{ii}|>\sum_{j\ne i}|a_{ij}| for every row. In that case the row sum norm of T_J is less than 1, so the spectral radius is less than 1, and Gauss-Seidel also converges for every starting vector.

Adding \varepsilon I to a matrix only increases the diagonal entries. Thus if a matrix is only weakly dominant, adding a positive multiple of the identity can make it strictly diagonally dominant, and then both methods converge from any initial guess.

If the spectral radius of the Jacobi matrix exceeds 1, then iterates diverge for suitable starting vectors. For example, if the error x_0-x^* is an eigenvector of T_J with eigenvalue of modulus greater than 1, then the error grows at every step. A one-step Jacobi computation is also a common exam task: identify D^{-1}b, compute -(D^{-1}A-I)x_0, and add.

Order of convergence of fixed-point iteration

For the fixed-point iteration x_{n+1}=g(x_n) with fixed point p=g(p), let e_n=x_n-p. Taylor expansion gives e_{n+1}=g'(p)e_n+\frac{g''(p)}{2}e_n^2+\cdots. The iteration has order m if |e_{n+1}|\approx C|e_n|^m with C\ne 0.

• If 0<|g'(p)|<1, the convergence is linear (order 1). If |g'(p)|>1, the fixed point is repelling and iterates do not converge to it from nearby points.
• If g^{(k)}(p)=0 for k=1,\dots,m-1 and g^{(m)}(p)\ne 0, then the order is m. In particular, g'(p)=0 and g''(p)=0 with g'''(p)\ne 0 gives order 3, because e_{n+1}\approx \frac{g'''(p)}{6}e_n^3.

A useful pattern: if g(x)=p+c(x-p)^m with m\ge 2, then the error satisfies e_{n+1}=c\,e_n^m and the order is m. Monotonicity also matters. If g is increasing, the iterates form a monotone sequence, and the direction of motion is determined by whether g(x)<x or g(x)>x on the interval. Between two fixed points where g(x)<x, iterates move down toward the lower fixed point, which can be attracting.

For a contraction, suppose g maps [a,b] into itself and |g'(x)|\le L<1 on [a,b]. Then g has a unique fixed point there, and the iteration converges from every starting point in the interval, with |x_n-p|\le \dfrac{L^n}{1-L}|x_1-x_0|. To verify the hypothesis, check both the range condition and the derivative bound.

Order of convergence of Newton-Raphson

Newton-Raphson iteration is x_{n+1}=x_n-\dfrac{f(x_n)}{f'(x_n)}, which is fixed-point iteration with g(x)=x-\dfrac{f(x)}{f'(x)}. For a simple root \alpha, g'(\alpha)=0, and generically g''(\alpha)\ne 0, so the convergence is quadratic.

For a root of multiplicity m, one has g'(\alpha)=\dfrac{m-1}{m}\ne 0, so Newton's method converges only linearly, with asymptotic error ratio \dfrac{m-1}{m}. The quadratic rate is restored by the modified iteration x_{n+1}=x_n-m\dfrac{f(x_n)}{f'(x_n)}, which is a standard remedy when the multiplicity is known.