TMV035 Analysis and Linear Algebra B, 2005

LECTURE 2.1

In this lecture we present the logarithm and numerical quadrature. This covers AMBS Ch 29, 30.

1. The logarithm

1.1. The construction. We fix a number b > 1 and consider the initial-value problem

\{ \begin{matrix}u′(x) = \frac{1}{x}, x ∈ [1, b], \\ u(1) = 0. \end{matrix} (1)

Recall that f(x) = 1/x is Lipschitz continuous on [δ, ∞) for any δ > 0. Therefore the Fundamental Theorem of Calculus gives a unique solution

u(x) = ∫_{1}^{x} \frac{1}{y} dy = lim_{n→∞} ∑_{j=1}^{i} \frac{1}{x_{j−1}^{n}} h_{n}, x ∈ [1, b]. (2)

where h_{n} = (b − 1)2^{−n}, x_{j}^{n} = 1 + jh_{n}, if i is related to n so that x = x_{i}^{n}, that is, i = (x − 1)/h_{n}. We have constructed a new function u(x). Since b is arbitrary it is defined for all x ≥ 1. It is called the natural logarithm and it is denoted

log(x) = ∫_{1}^{x} \frac{1}{y} dy, x ∈ [1, ∞). (3)

For 0 < x ≤ 1 we define it as the backward integral:

log(x) = ∫_{1}^{x} \frac{1}{y} dy = − ∫_{x}^{1} \frac{1}{y} dy, x ∈ (0, 1]. (4)

This function is denoted ln(x) in Swedish, French, German books. Thus, log(x) = ln(x). It is not the same as lg(x) = log_{10}(x), the base 10 logarithm, which we will introduce soon. In Matlab the natural logarithm is computed as

>> y=log(x)

To plot the graph on [1, 5]:

>> x=1:0.01:5;
>> y=log(x);
>> plot(x,y)

You can also compute and plot it with your own program my int.m:

>> [x,y]=my_int('funk',[1,5],0,1e-2)
>> plot(x,y)

where funk.m is

   function y=funk(x)
   y=1/x;

We now list the properties of log.

1.2. Domain of definition. It is clear that log(x) is defined for x > 0:

D_{log} = (0, ∞). (5)

Date: November 1, 2005, Stig Larsson, Mathematical Sciences, Chalmers University of Technology.

1.3. Derivative and initial value. It is clear from our construction, see (1), that log is uniformly differentiable on all intervals [a, b] with a > 0 and

D log(x) = \frac{1}{x}, log(1) = 0. (6)

Hence it is Lipschitz continuous on all such intervals [a, b].

1.4. Monotonicity. The logarithm is strictly increasing because D log(x) = \frac{1}{x} > 0 for x > 0. Hence log(x) > 0 for x > 1 and log(x) < 0 for 0 < x < 1.

1.5. Logarithm of product and quotient.

(7) log(ab) = log(a) + log(b), a, b > 0,

(8) log(a/b) = log(a) − log(b), a, b > 0.

Proof. By the chain rule and (6):

D log(ax) = \frac{1}{ax} a = \frac{1}{x},

which implies

log(ax) = log(x) + C.

We determine the constant by taking x = 1:

log(a) = log(1) + C = C, C = log(a).

Therefore

log(ax) = log(x) + log(a).

With x = b we get (7). The other formula (8) is obtained by using (7):

log(a) = log (\frac{a}{b} · b) = log (\frac{a}{b}) + log(b).
□

1.6. Logarithm of a power.

log(a^{n}) = n log(a), a > 0, n = 0, ±1, ±2, …, (9)

log(a^{r}) = r log(a), a > 0, r ∈ Q. (10)

Recall that the power function x^{r} with rational expontent r ∈ Q was defined for x ≥ 0 in AMBS 18.4, see Lecture 3.2 in ALA-A. We shall soon define it for r ∈ R.

Proof. We first consider r = n ∈ N. Repeated application of (7) with b = a gives:

log(a^{2}) = log(a) + log(a) = 2 log(a), log(a^{n}) = n log(a).

Also from (8):

log(a^{−1}) = log (\frac{1}{a}) = log(1) − log(a) = − log(a), log(a^{−n}) = −n log(a).

Since log(a^{0}) = log(1) = 0 we conclude (9). For the general case we make a substitution of variables:

log(a^{r}) = ∫_{1}^{a^{r}} \frac{1}{x} dx = \{ \begin{matrix}y = x^{1/r} \\ x = y^{r} \\ \frac{dx}{dy} = ry^{r−1} \\ dx = ry^{r−1} dy \end{matrix}\} = ∫_{1}^{a} \frac{1}{y^{r}} r y^{r−1}dy = r ∫_{1}^{a} \frac{1}{y} dy = r log(a).

(Recall the derivative \frac{dx^{r}}{dx} = rx^{r−1} from AMBS 24.10.) □

1.7. Asymptotic behavior.

lim_{x→∞} log(x) = ∞, (11)

lim_{x→0^{+}} log(x) = −∞. (12)

Here x → 0^{+} means that x → 0 with x > 0, i.e., x tends to zero from the right.

Proof. From (9), and since log(10) > 0, we get

log(10^{n}) = n log(10) → ∞ as n → ∞.

For any x > 1 we can find an integer n such that 10^{n} ≤ x < 10^{n+1}. Then

log(x) ≥ log(10^{n}) = n log(10).

We conclude (11) because x → ∞ implies n → ∞. For (12) we note that x^{−1} → ∞ as x → 0^{+} so that

log(x) = − log(x^{−1}) → −∞ as x → 0^{+}.
□

Thus log(x) tends to infinity as x tends to infinity. But the slope of the graph is 1/x which tends to zero. Therefore we expect that log(x) tends to infinity very slowly. In fact, it is possible to show that

(13)  lim_{x→∞} x^{−r} log(x) = 0 for any r > 0.

(We do not prove it here.) Here x^{−r} → 0 very slowly if r > 0 is small and log(x) → ∞, but nevertheless x^{−r} log(x) → 0. In other words: log(x) → ∞ more slowly than x^{r} for any r > 0.

1.8. Range. We have seen that the logarithm is strictly increasing from −∞ to ∞. Since it is Lipschitz continuous it cannot skip any values. We conclude that y = log(x) takes all values y between −∞ and ∞. Therefore

R_{log} = R. (14)

1.9. Graph. With the previous information we can now draw the graph of log. It is helpful to compute tangents, for example, for x = 0 the derivative is 1 so the tangent is y = x − 1. A Matlab plot can be seen in Figure 1.

1.10. An inequality.

log(1 + x) ≤ x, x > −1. (15)

The proof is left as an exercise, Problem 29.4. If we write the inequality in the equivalent form

log(x) ≤ x − 1, x > 0,

then Figure 1 indicates that it is true. This is not a proof, but it helps us to remember the inequality.

1.11. Logarithmic scale. Very large numbers and very small numbers are mapped to numbers of size near 1 by the logarithm. For example:

log(10^{10}) = 10 log(10) ≈ 30,
log(10^{−10}) = −10 log(10) ≈ −30.

This is useful when we want to graph data with the quotient of the largest to smallest number is very large. Instead of plotting y = f(x) we plot log(y) as a function of log(x) if both x and y span over large ranges. This is called a log-log plot. If only one of the variables in plotted in logarithmic scale then we have a semilog plot.

For example, the power function

y = ax^{r}, a, x > 0,

Figure 1. The graph y = log(x) together with the tangent y = x − 1.

becomes

log(y) = r log(x) + log(a),

which is a straight line with slope r in the log(x)-log(y) plane.

Example. We illustrate this by plotting the function y = 10x^{3} on the interval [10^{−3}, 10^{3}] in the original scale, see Figure 2, and in the log-log scale, see Figure 3. Note the large ranges of values spanned by both x and y: \frac{x_{max}}{x_{min}} = 10^{6} and \frac{y_{max}}{y_{min}} = 10^{18} . Note that the log-log plot has slope 3 and intersects the vertical axis at log(10) ≈ 3.

Figure 3 was created by

>> x=0.001:0.1:1000;
>> y=10*x.^3;
>> plot(log(x),log(y),'*-')
>> grid on
>> axis equal
>> xlabel('log(x)')
>> ylabel('log(y)')

Do this now with y = 10x^{−3}.

Matlab has functions for logarithmic plots: loglog, semilogx, semilogy. They are based on the base 10 logarithm lg = log_{10}. Figure 4 was created by

>> loglog(x,y)
>> grid on

2. Numerical quadrature

The word “quadrature” refers to computation of area. You may have heard of the classical problem of “the quadrature of the circle”, which is to find a square with the same area as a given circle. We know that computation of area is more or less the same as computation of integrals.

LECTURE 2.1 5

Figure 2. Plot of y = 10x^{3}. Figure 3. Log-log plot of y = 10x^{3}.

Therefore “numerical quadrature” means numerical computation of integrals. A “quadrature rule”

Figure 4. Log-log plot of y = 10x^{3}.

is an algorithm for the numerical computation of integrals. We have seen one such quadrature rule in the proof of the Fundamental Theorem of Calculus, namely the rectangle rule.

2.1. The rectangle rule. The Fundamental Theorem of Calculus constructs the integral as

∫_{a}^{b} f(x) dx = lim_{n→∞} ∑_{j=1}^{N} f(x_{j−1}^{n})h_{n}, h_{n} = (b − a)2^{−n}, x_{j}^{n} = a + jh_{n}, N = 2^{n}. (16)

We shall now estimate how fast the sum on the right converges to the integral as n → ∞, or equivalently, as h_{n} → 0. But first we replace the sum by a more general quadrature rule:

• variable steplength h_{j} = x_{j} − x_{j−1} and intervals I_{j} = [x_{j−1}, x_{j}];

• arbitrary interpolation point x̂_{j} ∈ I_{j}. Then we get

∫_{a}^{b} f(x) dx ≈ ∑_{j=1}^{N} f(x̂_{j})h_{j}. (17)

The sum on the right is called a Riemann sum. It can also be viewed as the general form of the rectangle rule. The sum in (16) is of this form with x̂_{j} = x_{j−1}^{n} and h_{j} = (b − a)2^{−n}.

We now consider the quadrature error:

Q_{h} = | ∫_{a}^{b} f(x) dx − ∑_{j=1}^{N} f(x̂_{j})h_{j}|. (18)

Theorem. (Error estimate for the general rectangle rule) Assume that f is differentiable on [a, b] with Lipschitz continuous and bounded derivative. Then

Q_{h} = | ∫_{a}^{b} f(x) dx − ∑_{j=1}^{N} f(x̂_{j})h_{j}| ≤ \frac{1}{2} max_{y∈[a,b]} |f′(y)|(b − a)h, (19)

where h = max_{1≤j≤N} h_{j}.

This means that the error converges to zero when the maximal steplength h tends to zero. The rate of convergence is

Q_{h} ≤ K_{1}h, K_{1} = \frac{1}{2} max_{y∈[a,b]} |f′(y)|(b − a). (20)

Thus, we gain one decimal of accuracy if we decrease the steplength by a factor 0.1. Another consequence of the theorem is that all Riemann sums

∑_{j=1}^{N} f(x̂_{j})h_{j}

converge to ∫_{a}^{b} f(x) dx as the maximal step h goes to zero, independently of the choice of interpolation point x̂_{j} and mesh intervals I_{j}. This proves that the integral is unique.

Proof. By writing ∫_{a}^{b} f(x) dx = ∑_{j=1}^{N} ∫_{x_{j−1}}^{x_{j}} f(x) dx = ∑_{j=1}^{N} ∫_{I_{j}} f(x) dx and using the triangle inequality for sums, we get

(21)  Q_{h} = | ∑_{j=1}^{N} ( ∫_{I_{j}} f(x) dx − f(x̂_{j})h_{j})|

≤ ∑_{j=1}^{N} | ∫_{I_{j}} f(x) dx − f(x̂_{j})h_{j}|. (22)

This expresses the global error as the sum of the local errors. We now look at the local error:

| ∫_{I_{j}} f(x) dx − f(x̂_{j}) \underbrace{h_{j}}_{=∫_{I_{j}} dx} | = | ∫_{I_{j}} f(x) dx − ∫_{I_{j}} f(x̂_{j}) dx|

= | ∫_{I_{j}} (f(x) − f(x̂_{j})) dx|
≤ ∫_{I_{j}} |f(x) − f(x̂_{j})| dx.

Here

|f(x) − f(x̂_{j})| = | ∫ x ; x̂_{j} f′(y) dy| ≤ max_{y∈I_{j}} |f′(y)||x − x̂_{j}|.

Hence

| ∫_{I_{j}} f(x) dx − f(x̂_{j})h_{j}| ≤ max_{y∈I_{j}} |f′(y)| ∫_{I_{j}} |x − x̂_{j}| dx ≤ max_{y∈I_{j}} |f′(y)|h_{j}^{2}.

Here we used the fact that |x − x̂_{j}| ≤ h_{j} so that ∫_{I_{j}} |x − x̂_{j}| dx ≤ h_{j}^{2}. A more detailed calculation shows that

\begin{aligned}∫_{I_{j}} |x − x̂_{j}| dx & = ∫ x̂_{j} ; x_{j−1}(x̂_{j} − x) dx + ∫ x_{j} ; x̂_{j}(x − x̂_{j}) dx = \frac{1}{2}((x̂_{j} − x_{j−1})^{2} + (x_{j} − x̂_{j})^{2}) \\ & = \{x̂_{j} = x_{j−1} + sh_{j}, 0 ≤ s ≤ 1\} = \frac{1}{2}(s^{2} + (1 − s)^{2})h_{j}^{2} ≤ \frac{1}{2}h_{j}^{2}. \end{aligned} (23)

With this slightly smaller estimate we get instead

| ∫_{I_{j}} f(x) dx − f(x̂_{j})h_{j}| ≤ max_{y∈I_{j}} |f′(y)| ∫_{I_{j}} |x − x̂_{j}| dx ≤ \frac{1}{2} max_{y∈I_{j}} |f′(y)|h_{j}^{2}. (24)

Finally we insert this bound of the local error into (21) to get

Q_{h} ≤ \frac{1}{2} ∑_{j=1}^{N} max_{y∈I_{j}} |f′(y)|h_{j}^{2} ≤ \frac{1}{2} max_{y∈[a,b]} |f′(y)|( ∑_{j=1}^{N} h_{j}) max_{1≤j≤N} h_{j} = \frac{1}{2} max_{y∈[a,b]} |f′(y)|(b − a)h,

because ∑_{j=1}^{N} h_{j} = b − a. □

2.2. The midpoint rule. The quantity in (23) is minimal (and = \frac{1}{4}h_{j}^{2}) if and only if s = \frac{1}{2}. This means that x̂_{j} = x_{j−1} + \frac{1}{2}h_{j} = \frac{1}{2}(x_{j−1} +x_{j}) is the midpoint of I_{j}. With this choice of interpolation point we get the so-called midpoint rectangle rule (or just the midpoint rule). It converges much faster than the general rectangle rule.

Theorem. (Error estimate for the midpoint rectangle rule) Assume that f is twice differentiable on [a, b] with Lipschitz continuous and bounded second derivative. Then, with x̂_{j} = x_{j−1} + \frac{1}{2}h_{j} = \frac{1}{2}(x_{j−1} + x_{j}), we have

Q_{h} = | ∫_{a}^{b} f(x) dx − ∑_{j=1}^{N} f(x̂_{j})h_{j}| ≤ \frac{1}{24}(b − a) max_{y∈[a,b]} |f′′(y)|h^{2}, (25)

where h = max_{1≤j≤N} h_{j}.

We skip the proof. It is based on Taylor’s formula.

The rate of convergence for the midpoint rule is

Q_{h} ≤ K_{2}h^{2}, K_{2} = \frac{1}{24}(b − a) max_{y∈[a,b]} |f′′(y)|. (26)

Thus, we gain two decimals of accuracy if we decrease the steplength by a factor 0.1.

2.3. The trapezoidal rule. The rectangle rule is based on integrating the piecewise constant function which is obtained by interpolating the function f in a point x̂_{j} ∈ I_{j}. This gives rise to

∫_{I_{j}} f(x) dx ≈ f(x̂_{j})h_{j}.

So the local integral is approximated by a small rectangle area. Another possibility is to interpolate f by the piecewise linear function:

f(x_{j−1})\frac{x_{j} − x}{h_{j}} + f(x_{j})\frac{x − x_{j−1}}{h_{j}} .

If we integrate this we get

∫_{I_{j}} f(x) dx ≈ \frac{1}{2}(f(x_{j−1}) + f(x_{j}))h_{j}.

This means that the local integral is approximated by the area of a small trapeze. For the global integral we get the trapezoidal rule:

∫_{a}^{b} f(x) dx = ∑_{j=1}^{N} ∫_{I_{j}} f(x) dx ≈ \frac{1}{2} ∑_{j=1}^{N} (f(x_{j−1}) + f(x_{j}))h_{j}

= \frac{1}{2}f(x_{0})h_{1} + \frac{1}{2} ∑_{j=1}^{N−1} f(x_{j})(h_{j} + h_{j+1}) + \frac{1}{2}f(x_{N})h_{N}.

It can be proved that it converges with the same rate as the midpoint rule:

Q_{h} ≤ Kh^{2}. (27)

2.4. Rate of convergence. We have seen that the quadrature error satisfies

Q_{h} ≤ Kh^{s}. (28)

where s = 1 for the general rectangle rule, and where s = 2 for the midpoint and trapezoidal rules. Such an inequality is best illustrated in a log-log scale:

(29) log(Q_{h}) ≤ log(K) + s log(h).

Here we used the fact that log is an increasing function so that the inequality is preserved. The exponent s can easily be estimated by computing the slope of the straight line.

2.5. Adaptivity. Recall the detailed estimate of the local error in (24):

| ∫_{I_{j}} f(x) dx − f(x̂_{j})h_{j}| ≤ \frac{1}{2} max_{y∈I_{j}} |f′(y)|h_{j}^{2}.

It shows that if the derivative is large on I_{j}, then this term will make a large contribution to the global error. This can be compensated by choosing h_{j} small. And vice versa: if max_{y∈I_{j}} |f′(y)| is small then we can choose h_{j} large. So we would take small steps where f changes rapidly and large steps where f changes slowly.

An algorithm which adjusts itself to the data of the problem in this way is called an adaptive algorithm.

For example, if we have an error tolerance TOL, then we might choose h_{j} so that

\frac{1}{2} max_{y∈I_{j}} |f′(y)|h_{j} ≈ \frac{TOL}{b − a}.

Then the global error becomes

Q_{h} ≤ \frac{1}{2} ∑_{j=1}^{N} max_{y∈I_{j}} |f′(y)|h_{j}^{2} = ∑_{j=1}^{N} \underbrace{\frac{1}{2} max_{y∈I_{j}} |f′(y)|h_{j}}_{≈ \frac{TOL}{b−a}} h_{j} ≈ \frac{TOL}{b − a}( ∑_{j=1}^{N} h_{j}) = TOL.

Further discussion of adaptivity can be found in AMBS 30.4. This is advanced, and can be skipped.

/stig