TMV035 Analysis and Linear Algebra B, 2005

LECTURE 1.1

In this lecture we present the integral. This covers AMBS Ch 27.

1. Introduction

Remember that the goal of ALA-A was to solve systems of algebraic equations:

f(x) = 0.

We introduced the following functions: the polynomials, the rational functions, and their inverse functions. The polynomial and rational functions can be computed exactly with rational numbers, but the inverse functions, for example, the square root \sqrt{x}, or more generally the power functions x^{p/q} with rational exponent, must be computed approximately by some iterative algorithm.

In this course the goal is to solve systems of ordinary differential equations (ODEs):

u′(x) = f(x, u(x)).

In this way we will introduce the functions ln(x), exp(x), sin(x), cos(x), and related functions. These are also based on algorithms for approximate computation. We will write Matlab programs for computing them.

We will also study linear algebra in n dimensions: matrix algebra, linear systems of equations, Gaussian elimination, orhtogonal matrices, least squares method.

We will do two applications projects together with the chemistry course: project 9, reaction kinetics, begins in week 4, examined by a written report; project 12, equilibrium equations, is done and examined in the studio class room in weeks 6–7.

2. The integral

We may consider two tasks in calculus:

1. Given u(x) determine u′(x). This was done in ALA-A.

2. Given u′(x) determine u(x). This is our first topic in ALA-B.

Example 1. We know from ALA-A that

Dx^{m} = mx^{m−1}, m ≥ 0, m ∈ Q.

Hence

u′(x) = D \frac{x^{m}}{m} = x^{m−1}, m > 0,

u′(x) = D \frac{x^{r+1}}{r + 1} = x^{r}, r = m − 1 > −1.

We may now guess that

u(x) = \frac{x^{r+1}}{r + 1} + c, r > −1,

where c is an arbitrary constant, because D(\frac{x^{r+1}}{r+1} + c) = x^{r}. We determine the constant:

u(0) = 0 ⟹ c = 0 ⟹ u(x) = \frac{x^{r+1}}{r + 1},

u(0) = u_{0} ⟹ c = u_{0} ⟹ u(x) = \frac{x^{r+1}}{r + 1} + u_{0}.

Date: October 23, 2005, Stig Larsson, Mathematical Sciences, Chalmers University of Technology.

Example 2. We also know from ALA-A that

Dx^{−m} = −mx^{−m−1}, m ≥ 0, m ∈ Q.

Hence

u′(x) = D \frac{x^{−m}}{−m} = x^{−m−1}, m > 0,

u′(x) = D \frac{x^{−r+1}}{−r + 1} = x^{−r}, r = m + 1 > 1,

We conclude that

u(x) = \frac{x^{−r+1}}{−r + 1} + c, r > 1,

where c is an arbitrary constant. We determine the constant:

u(1) = 0 ⟹ c = \frac{−1}{−r + 1} ⟹ u(x) = \frac{x^{−r+1} − 1}{−r + 1}

u(1) = u_{1} ⟹ c = \frac{−1}{−r + 1} + u_{1} ⟹ u(x) = \frac{x^{−r+1} − 1}{−r + 1} + u_{1}.

We now make the task more precise: Given an interval [a, b], a function f : [a, b] → R, and a number u_{a}, we want to find a unique function u such that

\{ \begin{matrix}u′(x) = f(x), x ∈ [a, b], \\ u(a) = u_{a}. \end{matrix} (1)

If we can find such a function u then we write

u(x) = u_{a} + ∫_{a}^{x} f(y) dy. (2)

The term ∫_{a}^{x} f(y) dy is called the integral of f over the interval [a, x]. The problem (1) is called an initial-value problem (“begynnelsevärdesproblem”). The first equation in (1) is a differential equation. The second equation is an initial condition and the number u_{a} is an initial value.

Example 1 means that

u(x) = u_{0} + ∫_{0}^{x} y^{r} dy = u_{0} + \frac{x^{r+1}}{r + 1} (r > 1)

is the solution of

\{ \begin{matrix}u′(x) = x^{r}, x ∈ [0, b], \\ u(0) = u_{0}, \end{matrix}

for any b > 0.

In a similar way we can integrate a constant function f(x) = K:

u(x) = u_{a} + ∫_{a}^{x} K dy = u_{a} + K(x − a), x ∈ [a, b]. (3)

This function satisfies

\{ \begin{matrix}u′(x) = K, x ∈ [a, b], \\ u(a) = u_{a}, \end{matrix}

for any b > a.

Example 2 means that

u(x) = u_{1} + ∫_{1}^{x} y^{−r} dy = u_{1} + \frac{x^{−r+1} − 1}{−r + 1} (r > 1)

is the solution of

\{ \begin{matrix}u′(x) = x^{−r}, x ∈ [1, b], \\ u(1) = u_{1}, \end{matrix}

for any b > 1.

In this way we can integrate x^{r} for all rational numbers r ≠ −1. What about x^{−1}? That is, ∫_{1}^{x} y^{−r} dy? This ought to be ln(x), but what is this? We will soon answer this question.

And what about an arbitrary function f(x)? That is, ∫_{1}^{x} f(y) dy, where f could be any function? The answer is that we must construct the function u(x) if we cannot guess it as we did in the

examples. We recall the constructive proof in four steps:

(1) An algoritm which produces a sequence.

(2) The sequence is a Cauchy sequence.

(3) The limit of the sequence solves the problem.

(4) The solution is unique.

Step 1. Algorithm. We begin by creating a mesh by bisecting the interval [a, b] repeatedly. When the number of bisections is n = 1 we have

x_{0}^{1} = a, x_{1}^{1} = \frac{a + b}{2}, x_{2}^{1} = b.

The points x_{i}^{1} are called meshpoints (or nodes) and the distance between them is the steplength h_{1} = (b − a)/2 = (b − a)2^{−1}. We continue with n = 2 and the steplength h_{2} = (b − a)2^{−2} and the meshpoints:

x_{0}^{2} = a, x_{1}^{2} = a + h_{2}, x_{2}^{2} = a + 2h_{2}, x_{3}^{2} = a + 3h_{2}, x_{4}^{2} = b.

Note that here the 2 in x_{i}^{2} is not an exponent but an index which indicates that we have bisected n = 2 times. With n = 3 the steplength is h_{3} = (b − a)2^{−3} and

x_{0}^{3} = a, x_{1}^{3} = a + h_{3}, x_{2}^{3} = a + 2h_{3}, x_{3}^{3} = a + 3h_{3}, x_{4}^{3} = a + 4h_{3},

x_{5}^{3} = a + 5h_{3}, x_{6}^{3} = a + 6h_{3}, x_{7}^{3} = a + 7h_{3}, x_{8}^{2} = b.

After n bisections we have h_{n} = (b − a)2^{−n} and

x_{0}^{n} = a, x_{1}^{n} = a + h_{n}, x_{2}^{n} = a + 2h_{n}, …, x_{i}^{n} = a + ih_{n}, …, x_{N}^{n} = b.

The meshpoints are {x_{i}^{n}}_{i=0}^{N}. Together they form a mesh on the interval [a, b]. The number of meshpoints is N + 1 and N = 2^{n}.

Then for each level of bisection n we create an approximate solution U^{n}(x) as follows. We first set

U^{n}(x_{0}^{n}) = U^{n}(a) = u_{a}.

Then we approximate the function f(x) on the interval [x_{0}^{n}, x_{1}^{n}] by its value f(x_{0}^{n}) in the left endpoint and integrate as in (3):

U^{n}(x) = U^{n}(x_{0}^{n}) + ∫_{x_{0}^{n}}^{x} f(x_{0}^{n}) dy
= U^{n}(x_{0}^{n}) + f(x_{0}^{n})(x − x_{0}^{n}), x ∈ [x_{0}^{n}, x_{1}^{n}].

At the next meshpoint we then have

U^{n}(x_{1}^{n}) = U^{n}(x_{0}^{n}) + f(x_{0}^{n})(x_{1}^{n} − x_{0}^{n})
= U^{n}(x_{0}^{n}) + f(x_{0}^{n})h_{n}.

We continue in the same way:

U^{n}(x) = U^{n}(x_{1}^{n}) + ∫_{x_{1}^{n}}^{x} f(x_{1}^{n}) dy
= U^{n}(x_{1}^{n}) + f(x_{1}^{n})(x − x_{1}^{n}), x ∈ [x_{1}^{n}, x_{2}^{n}],

and

U^{n}(x_{2}^{n}) = U^{n}(x_{1}^{n}) + f(x_{1}^{n})h_{n}.

In this way we construct U^{n}(x) in each mesh interval I_{i}^{n} = [x_{i−1}^{n}, x_{i}^{n}]. It is sufficient to compute the node values U^{n}(x_{i−1}^{n}) and U^{n}(x_{i}^{n}); then the intermediate values are obtained by

U^{n}(x) = U^{n}(x_{i−1}^{n}) + f(x_{i−1}^{n})(x − x_{i−1}^{n})

= U^{n}(x_{i−1}^{n}) + \frac{U^{n}(x_{i}^{n}) − U^{n}(x_{i−1}^{n})}{h_{n}}(x − x_{i−1}^{n}), x ∈ [x_{i−1}^{n}, x_{i}^{n}].

In other words: the graph of U^{n} is obtained by drawing a straight line between the node values. We can now formulate the algorithm in a compact form:

Algorithm. (Rectangle rule) First set the initial values:

\{ \begin{matrix}x_{0}^{n} = a, \\ U(x_{0}^{n}) = u_{a}. \end{matrix} (4)

Then set h_{n} = 2^{−n}(b − a), N = 2^{n}, and compute for n = 1, …, N:

\{ \begin{matrix}x_{i}^{n} = x_{i−1}^{n} + h_{n}, \\ U^{n}(x_{i}^{n}) = U^{n}(x_{i−1}^{n}) + f(x_{i−1}^{n})h_{n}. \end{matrix} (5)

This is very easy to program in Matlab. You will soon do this in the studio class.

Another motivation for the algorithm is that if a solution u exists, then it is differentiable with u′ = f so that

u(x_{i}^{n}) = u(x_{i−1}^{n}) + u′(x_{i−1}^{n})(x_{i}^{n} − x_{i−1}^{n}) + E_{u}(x_{i}^{n}, x_{i−1}^{n})
= u(x_{i−1}^{n}) + f(x_{i−1}^{n})h_{n} + E_{u}(x_{i}^{n}, x_{i−1}^{n}).

If we skip the remainder then we obtain (5).

It is important to note that U^{n}(x_{i}^{n}) can also be expressed as a sum:

U^{n}(x_{i}^{n}) = U^{n}(x_{i−1}^{n}) + f(x_{i−1}^{n})h_{n}
= U^{n}(x_{i−2}^{n}) + f(x_{i−2}^{n})h_{n} + f(x_{i−1}^{n})h_{n}
⋮

= U^{n}(x_{0}^{n}) + f(x_{0}^{n})h_{n} + ⋯ + f(x_{i−1}^{n})h_{n} (6)

= u_{a} + ∑_{j=1}^{i} f(x_{j−1}^{n})h_{n},

that is,

U^{n}(x_{i}^{n}) = u_{a} + ∑_{j=1}^{i} f(x_{j−1}^{n})h_{n}. (7)

Step 2. Cauchy sequence. We must show that {U^{n}(x)}_{n=1}^{∞} is a Cauchy sequence (for each fixed x), i.e.,

|U^{n}(x) − U^{m}(x)| → 0, as m, n → ∞.

The proof of this is rather long and complicated. It is written in detail in the book. In the proof it is important that f is Lipschitz continuous on [a, b]. You may skip this part of the proof if you like.

That {U^{n}(x)}_{n=1}^{∞} is a Cauchy sequence means that it generates a decimal expansion, and a decimal expansion is the same as a real number. Therefore we get a real number u(x) such that

u(x) = lim_{n→∞} U^{n}(x).

Using (7) this can also be expressed as

u(x) = u_{a} + lim_{n→∞} ∑_{j=1}^{i} f(x_{j−1}^{n})h_{n},

if i is related to n so that x = x_{i}^{n}, that is, i = \frac{x − a}{h_{n}} = \frac{x − a}{b − a}2^{n}.

(8)

We remark that this is possible (i.e., i is an integer) only if x is a meshpoint. But the set of all meshpoints x_{i}^{n} for all n and i is very large and densely spread out in [a, b], in fact all real numbers x can be approximated by a sequence of meshpoints. We don’t prove this.

We have now constructed a new function u(x).

Step 3. The limit solves the equation. We have to show that the new new function u(x) satisfies the initial-value problem (1). First of all, we have U^{n}(a) = u_{a} so that u(a) = lim_{n→∞} U^{n}(a) = u_{a}, so that the initial condition is satisfied. Then we must show that u is differentiable and that u′(x) = f(x). We will do this later, AMBS 28.8.

Step 4. Uniqueness. We must show that u is the only solution of (1). Assume then that v is another solution. This means that

\{ \begin{matrix}v′(x) = f(x), x ∈ [a, b], \\ v(a) = u_{a}. \end{matrix}

Since v is differentiable we have as in (6) and (7)

v(x) = v(x_{i}^{n}) = v(x_{i−1}^{n}) + v′(x_{i−1}^{n})(x_{i}^{n} − x_{i−1}^{n}) + E_{v}(x_{i}^{n}, x_{i−1}^{n})

= v(x_{i−1}^{n}) + f(x_{i−1}^{n})h_{n} + E_{v}(x_{i}^{n}, x_{i−1}^{n})

= u_{a} + ∑_{j=1}^{i} f(x_{j−1}^{n})h_{n} + ∑_{j=1}^{i} E_{v}(x_{j}^{n}, x_{j−1}^{n}),

if i is related to n so that x = x_{i}^{n}, that is, i = \frac{x − a}{b − a}2^{n} = \frac{x − a}{h_{n}}. We know from (8) that

u_{a} + ∑_{j=1}^{i} f(x_{j−1}^{n})h_{n} → u(x) as n → ∞.

For the remainder we have

| ∑_{j=1}^{i} E_{v}(x_{j}^{n}, x_{j−1}^{n})| ≤ ∑_{j=1}^{i} |E_{v}(x_{j}^{n}, x_{j−1}^{n})| ≤ ∑_{j=1}^{i} K_{v}|x_{j}^{n} − x_{j−1}^{n}|^{2} = ∑_{j=1}^{i} K_{v}h_{n}^{2}

= K_{v}h_{n}^{2}i = K_{v}(ih_{n})h_{n} = K_{v}(x − a)h_{n} \{because i = \frac{x − a}{h_{n}}\}

≤ K_{v}h_{n}(b − a) → 0.

Therefore

v(x) = v(x_{i}^{n}) = u_{a} + ∑_{j=1}^{i} f(x_{j−1}^{n})h_{n} + ∑_{j=1}^{i} E_{v}(x_{j}^{n}, x_{j−1}^{n}) → u(x) as n → ∞.

We conclude that v(x) = u(x) so there is only one solution.

The fundamental theorem of calculus. We can now formulate what we have done as a theorem.

Theorem. Assume that f : [a, b] → R is Lipschitz continuous and u_{a} is a number. Then there is a unique function u : [a, b] → R such that

\{ \begin{matrix}u′(x) = f(x), x ∈ [a, b], \\ u(a) = u_{a}. \end{matrix}

6 LECTURE 1.1

This function can be written u(x) = u_{a} + ∫_{a}^{x} f(y) dy, (9) where the integral is constructed as the limit

∫_{a}^{x} f(y) dy = lim_{n→∞} ∑_{j=1}^{i} f(x_{j−1}^{n})h_{n},

where h_{n} = (b − a)2^{−n}, x_{j}^{n} = a + jh_{n}, if i is related to n so that x = x_{i}^{n}, that is, i = \frac{x − a}{h_{n}} =

\frac{x − a}{b − a}2^{n}.

The function u is uniformly differentiable on [a, b] with constant K_{u} = \frac{1}{2}L_{f}.

The function u is a primitive function of f. A primitive function of f is a function F such that F′ = f. It is only determined up to a constant: F (x) + C is also a primitive function. It can be made unique by adding an initial condition. Example. f(x) = x^{2} and F (x) = \frac{1}{3}x^{3} is a primitive function. We set u(x) = F (x) + C = \frac{1}{3}x^{3} + C and u(1) = 1. Then C = \frac{2}{3} and u(x) = \frac{1}{3}x^{3} + \frac{2}{3}.

Sometimes we write (9) (with x = b) as

u(b) = u(a) + ∫_{a}^{b} f(x) dx = u(a) + ∫_{a}^{b} u′(x) dx,

which means that the integral of u′ can be expressed as

∫_{a}^{b} u′(x) dx = u(b) − u(a) = [u(x)]_{x=a}^{b}.

Here we introduced the convenient notation

[u(x)]_{x=a}^{b} = u(b) − u(a).

Example.

∫_{1}^{2} x^{2} dx = [\frac{x^{3}}{3} ]_{x=1}^{2} = \frac{8}{3} − \frac{1}{3} = \frac{7}{3}.

It is important to note that the integration variable y in (9) is a “dummy variable” which can be replaced by anything but x and a:

∫_{a}^{x} f(y) dy = ∫_{a}^{x} f(s) ds = ∫_{a}^{x} f(ö) dö.

Example.

∫_{1}^{x} y^{2} dy = [\frac{y^{3}}{3} ]_{y=1}^{x} = \frac{x^{3}}{3} − \frac{1}{3},

∫_{1}^{x} ö^{2} dö = [\frac{ö^{3}}{3} ]_{ö=1}^{x} = \frac{x^{3}}{3} − \frac{1}{3}.

The same thing holds for sums:

∑_{i=1}^{N} a_{i} = ∑_{n=1}^{N} a_{n} = a_{1} + ⋯ + a_{N}.

Interpretation as area. Assume now that f : [a, b] → R_{+} is a nonnegative function. Then the sum in (9) (with x = b = x_{N}^{n})

∑_{j=1}^{N} f(x_{j−1}^{n})h_{n}

is a sum of small rectangle areas f(x_{j−1}^{n})h_{n}, each rectangle with base h_{n} and height f(x_{j−1}^{n}). When n → ∞ and h_{n} → 0 this sum converges to the area under the graph y = f(x), x ∈ [a, b]. On the other hand we know that the limit is

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

We conclude that

A = ∫_{a}^{b} f(x) dx

is the area under the graph.

Example. The area under y = x^{2}, x ∈ [1, 2], is

A = ∫_{1}^{2} x^{2} dx = [\frac{x^{3}}{3} ]_{x=1}^{2} = \frac{8}{3} − \frac{1}{3} = \frac{7}{3}.

/stig