TMV035 Analysis and Linear Algebra B

Linearization. Jacobi matrix. Newton’s method.

Studio 1.2

Läraren går igenom avsnitt 1.1–1.4 nedan. Du jobbar med Problem 1.1–1.5. (Huvudsakligen handräkning.)

Studio 2.1

Du skriver ett program jacobi.m med anropet A=jacobi(f,x) som beräknar jacobimatrisen till funktionen f i punkten x. Du skall använda programskalet \underline{jacobi.m}. Testa programmet på Problem 1.2–1.5.

Studio 2.2

Läraren går igenom avsnitt 1.5 nedan. Du skriver ett program newton.m med anropet x=newton(f,x0,tol) som löser ekvationssystemet f(x) = 0. Du skall använda programskalet \underline{newton.m}. Testa programmet på Problem 1.6–1.7.

Introduction

The fixed point iteration (and hence also Newton’s method) works equally well for systems of equations. For example,

x_{2}(1 − x_{1}^{2}) = 0,
2 − x_{1}x_{2} = 0,

is a system of two equations in two unknowns. See Problem 1.6 below. If we define two functions

f_{1}(x_{1}, x_{2}) = x_{2}(1 − x_{1}^{2}),
f_{2}(x_{1}, x_{2}) = 2 − x_{1}x_{2},

then the equations may be written

f_{1}(x_{1}, x_{2}) = 0,
f_{2}(x_{1}, x_{2}) = 0.

With f = (f_{1}, f_{2}), x = (x_{1}, x_{2}), and 0 = (0, 0), we note that f : R^{2} → R^{2} and we can write the equations in the compact form

f(x) = 0.

In this lecture we will see how Newton’s method can be applied to such systems of equations.

Note that the bisection algorithm can only be used for a single equation, but not for a system of several equations. This is because it relies on the fact the the graph of a Lipschitz continuous function f : R → R must pass the value zero if it is positive in one point and negative in another point. This has no counterpart for functions f : R^{2} → R^{2}.

Before we discuss Newton’s method we need to define derivatives of such functions, namely, two functions of two variables, and more generally several functions of several variables.

1.1 Function of one variable, f : R → R

A function f : R → R of one variable is differentiable at a if the following limit exists:

f′(a) = lim_{x→a} \frac{f(x) − f(a)}{x − a}.

We write this in an equivalent form:

lim_{x→a} \frac{f(x) − f(a) − f′(a)(x − a)}{x − a} = 0.

Therefore we can say that a function f : R → R of one variable is differentiable at a if there is a function m(a), such that

lim_{x→a} \frac{f(x) − f(a) − m(a)(x − a)}{x − a} = 0.  (1)

Of course the function m(a) is the derivative of f at a:

m(a) = f′(a) = Df(a) = \frac{df}{dx}(a).

This formulation will be useful when we define the derivative of a function of two variables later. We also obtain the linearization formula

f(x) = f(a) + f′(a)(x − a) + E_{f}(x, a),  (2)

where the linearization error E_{f} is smaller than the second term on the right side when x is close to a.

It is convenient to use the abbreviation h = x − a, so that x = a + h and (1) becomes

lim_{h→0} \frac{f(a + h) − f(a) − m(a)h}{h} = 0,  (3)

and (2) becomes

f(x) = f(a + h) = f(a) + f′(a)h + E_{f}(x, a).  (4)

Note that the first term on the right side, f(a), is constant with respect to x. The second term,

f′(a)h = f′(a)(x − a),  (5)

is a linear function of the increment h = x − a. These two terms form the linearization of f at a,

L(x) = f(a) + f′(a)(x − a).  (6)

The straight line y = L(x) is the tangent to the curve y = f(x) at a.

Example 1. Let f(x) = x^{2}. Then f′(x) = 2x and the linearization at a = 3 is

L(x) = 9 + 6(x − 3).

Numerical computation of the derivative

In a previous studio lecture \underline{Numerical computation of derivatives} we learnt how to compute the derivative numerically. We quickly repeat it here. If we divide (4) by h, then we get

\frac{f(a + h) − f(a)}{h} = f′(a) + E_{f}(x, a)/h.  (7)

Here the remainder E_{f}(x, a)/h → 0 when h → 0. This suggests that we can approximate the derivative by the difference quotient

f′(a) ≈ \frac{f(a + h) − f(a)}{h}.  (8)

A better approximation is obtained by the symmetric difference quotient:

f′(a) ≈ \frac{f(a + h) − f(a − h)}{2h}.  (9)

The difference quotients in (8) and (9) are of the form ”small number divided by small number”. If this is computed with round-off error on a computer, then the total error will be large if the step h is very small. Therefore we must choose the step “moderately small” here, see \underline{Numerical computation of derivatives}. It can be shown that in Matlab a good choice for (8) is h = 10^{−8} and for (9) h = 10^{−5}.

1.2 Function of two variables, f : R^{2} → R

Let f(x_{1}, x_{2}) be a function of two variables, i.e., f : R^{2} → R. We want to imitate the formula in (3). We write x = (x_{1}, x_{2}) and f(x) = f(x_{1}, x_{2}) and introduce the increment vector h = (h_{1}, h_{2}) = (x_{1} − a_{1}, x_{2} − a_{2}) and its length |h| = \sqrt{h_{1}^{2} + h_{2}^{2}}.

We now say that function f is differentiable at a = (a_{1}, a_{2}), if there are functions m_{1}(a), m_{2}(a), such that

lim_{|h|→0} \frac{f(a + h) − f(a) − m_{1}(a)h_{1} − m_{2}(a)h_{2}}{|h|} = 0.  (10)

The corresponding linearization formula is

f(x) = f(a + h) = f(a) + m_{1}(a)h_{1} + m_{2}(a)h_{2} + E_{f}(x, a),  (11)

where the linearization error E_{f} is smaller than the second and third terms on the right side, more precisely, E_{f}(x, a)/|h| → 0 as |h| → 0.

If we take h = (h_{1}, 0), then we get

f(x_{1}, a_{2}) = f(a_{1} + h_{1}, a_{2}) = f(a) + m_{1}(a)h_{1} + E_{f}(x, a).

By comparison with (4) we see that this means that m_{1}(a) is the derivative of the one-variable function f̂(x_{1}) = f(x_{1}, a_{2}), obtained from f by keeping x_{2} = a_{2} fixed. By taking h = (0, h_{2}) we see in a similar way that m_{2}(a) is the derivative of the one-variable function, which is obtained from f by keeping x_{1} = a_{1} fixed. The functions m_{1}(a), m_{2}(a) are called the partial derivatives of f at a and we denote them by

m_{1}(a) = f′_{x_{1}}(a) = \frac{∂f}{∂x_{1}} (a), m_{2}(a) = f′_{x_{2}}(a) = \frac{∂f}{∂x_{2}} (a).  (12)

Now (11) may be written

f(x) = f(a + h) = f(a) + f′_{x_{1}}(a)h_{1} + f′_{x_{2}}(a)h_{2} + E_{f}(x, a), h = x − a.  (13)

It is convenient to write this formula by means of matrix notation. Let

a = [a_{1}, a_{2}] , b = [\begin{matrix}b_{1} \\ b_{2} \end{matrix}] .

We say that a is a row matrix of type 1 × 2 (one by two) and that b is a column matrix of type 2 × 1 (two by one). Their product is defined by

ab = [a_{1}, a_{2}] [\begin{matrix}b_{1} \\ b_{2} \end{matrix}] = a_{1}b_{1} + a_{2}b_{2}.

The result is a matrix of type 1×1 (one real number), according to the rule: (1×2)(2×1) = 1×1. Going back to (13) we define

f′(a) = Df(a) = [f′_{x_{1}}(a) f′_{x_{2}}(a)] , h = [\begin{matrix}h_{1} \\ h_{2} \end{matrix}] .

The row matrix f′(a) = Df(a) = [f′_{x_{1}}(a) f′_{x_{2}}(a)] is called the derivative (or Jacobi matrix) of f at a. Then (13) may be written

\begin{aligned}f(x) = f(a + h) & = f(a) + [f′_{x_{1}}(a) f′_{x_{2}}(a)] [\begin{matrix}h_{1} \\ h_{2} \end{matrix}] + E_{f}(x, a) \\ & = f(a) + f′(a)h + E_{f}(x, a), h = x − a. \end{aligned}  (14)

Note that the first term on the right side, f(a), is constant with respect to x. The second term,

f′(a)h = f′(a)(x − a),  (15)

is a linear function of the increment h = x − a. These terms are called the linearization of f at a,

L(x) = f(a) + f′(a)(x − a).  (16)

The plane x_{3} = L(x_{1}, x_{2}) is the tangent plane to the surface x_{3} = f(x_{1}, x_{2}) at a.

Example 2. Let f(x) = x_{1}^{2}x_{2}^{5}. Then

\frac{∂f}{∂x_{1}} (x) = \frac{∂f}{∂x_{1}} (x_{1}^{2}x_{2}^{5}) = 2x_{1}x_{2}^{5}, \frac{∂f}{∂x_{2}} (x) = \frac{∂f}{∂x_{2}} (x_{1}^{2}x_{2}^{5}) = 5x_{1}^{2}x_{2}^{4},

so that f′(x) = [2x_{1}x_{2}^{5} 5x_{1}^{2}x_{2}^{4}] and the linearization at a = (3, 1) is

L(x) = 9 + [6 45] [\begin{matrix}x_{1} − 3 \\ x_{2} − 1 \end{matrix}] .

1.3 Two functions of two variables, f : R^{2} → R^{2}

Let f_{1}(x_{1}, x_{2}), f_{2}(x_{1}, x_{2}) be two functions of two variables. We write x = (x_{1}, x_{2}) and f(x) = (f_{1}(x_{1}, x_{2}), f_{2}(x_{1}, x_{2})), i.e., f : R^{2} → R^{2}. The function f is differentiable at a = (a_{1}, a_{2}), if there are constants m_{11}(a), m_{12}(a), m_{21}(a), m_{22}(a), and K_{f}(a) such that

\begin{aligned}f_{1}(x) & = f_{1}(a + h) = f_{1}(a) + m_{11}(a)h_{1} + m_{12}(a)h_{2} + E_{f_{1}}(x, a), \\ f_{2}(x) & = f_{2}(a + h) = f_{2}(a) + m_{21}(a)h_{1} + m_{22}(a)h_{2} + E_{f_{2}}(x, a), \end{aligned}  (17)

where h = x−a and the linearization errors E_{f_{j}} satisfy E_{f_{j}}(x, a)/|h| → 0 when |h| → 0. As before |h| = \sqrt{h_{1}^{2} + h_{2}^{2}} denotes the norm (length) of the increment vector h = (h_{1}, h_{2}) = (x_{1}−a_{1}, x_{2}−a_{2}). From the previous subsection we recognize that the constants m_{ij}(a) are the partial derivatives of the functions f_{i} at a and we denote them by

\begin{aligned}m_{11}(a) & = f′_{1,x_{1}}(a) = \frac{∂f_{1}}{∂x_{1}} (a), m_{12}(a) = f′_{1,x_{2}}(a) = \frac{∂f_{1}}{∂x_{2}} (a), \\ m_{21}(a) & = f′_{2,x_{1}}(a) = \frac{∂f_{2}}{∂x_{1}} (a), m_{22}(a) = f′_{2,x_{2}}(a) = \frac{∂f_{2}}{∂x_{2}} (a). \end{aligned}

It is convenient to use matrix notation. Let

A = [\begin{matrix}a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix}] , b = [\begin{matrix}b_{1} \\ b_{2} \end{matrix}] .

We say that A is a matrix of type 2 × 2 (two by two) and that b is a column matrix of type 2 × 1 (two by one). Their product is defined by

Ab = [\begin{matrix}a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix}] [\begin{matrix}b_{1} \\ b_{2} \end{matrix}] = [\begin{matrix}a_{11}b_{1} + a_{12}b_{2} \\ a_{21}b_{1} + a_{22}b_{2} \end{matrix}] .

The result is a matrix of type 2 × 1 (column matrix), according to the rule: (2 × 2)(2 × 1) = 2 × 1.

Going back to (17) we define

f(x) = [\begin{matrix}f_{1}(x) \\ f_{2}(x) \end{matrix}] , f′(a) = Df(a) = [\begin{matrix}\frac{∂f_{1}}{∂x_{1}} (a) & \frac{∂f_{1}}{∂x_{2}} (a) \\ \frac{∂f_{2}}{∂x_{1}} (a) & \frac{∂f_{2}}{∂x_{2}} (a) \end{matrix}] , h = [\begin{matrix}h_{1} \\ h_{2} \end{matrix}] .  (18)

The matrix f′(a) = Df(a) is called the derivative (or Jacobi matrix) of f at a. Then (17) may be written

[\begin{matrix}f_{1}(x) \\ f_{2}(x) \end{matrix}] = [\begin{matrix}f_{1}(a + h) \\ f_{2}(a + h) \end{matrix}] = [\begin{matrix}f_{1}(a) \\ f_{2}(a) \end{matrix}] + [\begin{matrix}\frac{∂f_{1}}{∂x_{1}} (a) & \frac{∂f_{1}}{∂x_{2}} (a) \\ \frac{∂f_{2}}{∂x_{1}} (a) & \frac{∂f_{2}}{∂x_{2}} (a) \end{matrix}] [\begin{matrix}h_{1} \\ h_{2} \end{matrix}] + [\begin{matrix}E_{f_{1}}(x, a) \\ E_{f_{2}}(x, a) \end{matrix}] ,  (19)

or in more compact form

f(x) = f(a + h) = f(a) + f′(a)h + E_{f}(x, a), h = x − a.  (20)

Note that it is important that f, x, a, h are written as column vectors here.

Note that the first term on the right side, f(a), is constant with respect to x. The second term,

f′(a)h = f′(a)(x − a),  (21)

is a linear function of the increment h = x − a. These terms are called the linearization of f at a,

L(x) = f(a) + f′(a)(x − a).  (22)

Example 3. Let f(x) = [\begin{matrix}x_{1}^{2}x_{2}^{5} \\ x_{2}^{3} \end{matrix}]. Then

f′(x) = Df(x) = [\begin{matrix}\frac{∂f_{1}}{∂x_{1}} (x) & \frac{∂f_{1}}{∂x_{2}} (x) \\ \frac{∂f_{2}}{∂x_{1}} (x) & \frac{∂f_{2}}{∂x_{2}} (x) \end{matrix}] = [\begin{matrix}2x_{1}x_{2}^{5} & 5x_{1}^{2}x_{2}^{4} \\ 0 & 3x_{2}^{2} \end{matrix}]

and the linearization at a = (3, 1) is

L(x) = [\begin{matrix}9 \\ 1 \end{matrix}] + [\begin{matrix}6 & 45 \\ 0 & 3 \end{matrix}] [\begin{matrix}x_{1} − 3 \\ x_{2} − 1 \end{matrix}] .

1.4 Several functions of several variables, f : R^{n} → R^{m}

It is now easy to generalize to any number of functions in any number of variables. Let f_{i} be m functions of n variables x_{j}, i.e., f : R^{n} → R^{m}. As in (18) we define

x = [x_{1} ; ⋮ ; x_{n}] , h = [h_{1} ; ⋮ ; h_{n}] = [x_{1} − a_{1} ; ⋮ ; x_{n} − a_{n}] ,

f(x) = [ f_{1}(x_{1}, …, x_{n}) ; ⋮ ; f_{m}(x_{1}, …, x_{n})] , f′(a) = Df(a) = [\begin{matrix}\frac{∂f_{1}}{∂x_{1}} (a) & … & \frac{∂f_{1}}{∂x_{n}} (a) \\ ⋮ & & ⋮ \\ \frac{∂f_{m}}{∂x_{1}} (a) & … & \frac{∂f_{m}}{∂x_{n}} (a) \end{matrix}] .

The m × n matrix f′(a) = Df(a) is called the derivative (or Jacobi matrix) of f at a. In a similar way to (20) we get

f(x) = f(a + h) = f(a) + f′(a)h + E_{f}(x, a), h = x − a.  (23)

The linearization of f at a is

L(x) = f(a) + f′(a)(x − a).  (24)

Numerical computation of the derivative. In order to compute the j-th column \frac{∂f}{∂x_{j}} (a) of

the Jacobi matrix, we choose the increment h such that h_{j} = δ and h_{i} = 0 for i ≠ j, i.e.,

h = [0 ; ⋮ ; 0 ; δ ; 0 ; ⋮ ; 0] = δ [0 ⋮ 0 1 0 ⋮ 0] = δe_{j}, e_{j} = [0 ⋮ 0 1 0 ⋮ 0] ← row no. j.

Here the steplength δ is a small positive number and e_{j} is the j-th standard basis vector. If we use this increment in a symmetric difference quotient, see (9), we get

\frac{∂f}{∂x_{j}} (a) ≈ \frac{f(a + δe_{j}) − f(a − δe_{j})}{2δ}.  (25)

Note that f is a column so the result is a column: the j^{th} column of the matrix f′(a). Remember that the steplength δ should be small, but not too small.

1.5 Newton’s method for f(x) = 0

Consider a system of n equations with n unknowns:

f_{1}(x_{1}, …, x_{n}) = 0,
⋮
f_{n}(x_{1}, …, x_{n}) = 0.

If we define

x = [x_{1} ; ⋮ ; x_{n}] , f = [f_{1} ; ⋮ ; f_{n}] , 0 = [0 ⋮ 0] ,

then f : R^{n} → R^{n}, and we can write our system of equations in the compact form

f(x) = 0.  (26)

Suppose that we have found an approximate solution a. We want to find a better approximation x = a + h. Instead of solving (26) directly, which is usually impossible, we solve the linearized equation at a:

L(a + h) = f(a) + f′(a)h = 0.  (27)

We must solve for the increment h. Rearranging the terms we get

f′(a)h = −f(a).  (28)

Remember that the Jacobi matrix is of type n × n and the increment is of type n × 1. Therefore we have to solve a linear system of n equations with n unknowns to get the increment h. It is of the form Ah = b with A = f′(a) and b = −f(a). Then we set x = a + h.

In algorithmic form Newton’s method can be formulated:

while |h|>tol
   evaluate the residual    b=-f(x)
   evaluate the Jacobian    A=f'(x)
   solve the linear system  Ah=b
   update                   x=x+h
end

You will implement this algorithm in the studio exercises. You will use the Matlab command

h=A\b

to solve the system. But later in this course we will study linear systems of equations of the form Ah = b and we will answer important questions such as:

• Is there a unique solution h for each b?

• How do you compute the solution?

These questions can be answered for linear systems Ah = b, but not for the more general nonlinear systems f(x) = 0. Thus, Newton’s method transforms the task of solving a difficult equation to the task of solving an easier equation many times. The study of systems of linear equations is an important part of the subject “linear algebra” which we study in ALA-B.

90 Problems

Problem 1.1. Let

a = [1 2] , b = [\begin{matrix}1 \\ 2 \end{matrix}] , A = [\begin{matrix}1 & 2 \\ 3 & 4 \end{matrix}] .

Compute the products ab, Ab, Aa.

Problem 1.2. Compute the Jacobi matrix f′(x) (also denoted Df(x)). Compute the linearization of f at x̄.

(a) f(x) = [\begin{matrix}sin(x_{1}) + cos(x_{2}) \\ cos(x_{1}) + sin(x_{2}) \end{matrix}] , x̄ = 0; (b) f(x) = [\begin{matrix}1 \\ 1 + x_{1} \\ 1 + x_{1}e^{x_{2}} \end{matrix} ] , x̄ = [\begin{matrix}1 \\ 1 \end{matrix}] .

Problem 1.3. Compute the gradient vector ∇f(x) (also denoted f′(x) = Df(x)). Compute the linearization of f at x̄.

(a) f(x) = e^{−x_{1}} sin(x_{2}), x̄ = 0; (b) f(x) = |x|^{2} = x_{1}^{2} + x_{2}^{2} + x_{3}^{2}, x ∈ R^{3}, x̄ = [\begin{matrix}1 \\ 1 \\ 1 \end{matrix}] .

Problem 1.4. Here f : R → R^{2}. Compute f′(t). Compute the linearization of f at t̄.

(a) f(t) = [\begin{matrix}cos(t) \\ sin(t) \end{matrix}] , t̄= π/2; (b) f(t) = [\begin{matrix}t \\ 1 + t^{2} \end{matrix}] , t̄= 0.

Problem 1.5. Compute the linearization of f at x̄:

f(x) = [\begin{matrix}x_{1} − x_{1}x_{2} \\ −x_{2} + x_{1}x_{2} \end{matrix}] , x̄ = [\begin{matrix}1 \\ 1 \end{matrix}] .

Problem 1.6. (a) Write the system

u_{2}(1 − u_{1}^{2}) = 0,
2 − u_{1}u_{2} = 0

in the form f(u) = 0. Find the all the solutions by hand calculation. (b) Compute the Jacobi matrix Df(u).

(c) Perform the first step of Newton’s method for the equation f(u) = 0 with initial vector

u^{(0)} = [\begin{matrix}1 \\ 1 \end{matrix}].

(d) Solve the equation f(u) with your Matlab program newton.m.

Problem 1.7. (a) Write the system

u_{1}(1 − u_{2}) = 0,
u_{2}(1 − u_{1}) = 0,

in the form f(u) = 0. Find the all the solutions by hand calculation. (b) Compute the Jacobi matrix Df(u).

(c) Perform the first step of Newton’s method for the equation f(u) = 0 with initial vector

u^{(0)} = [\begin{matrix}2 \\ 2 \end{matrix}].

(d) Solve the equation f(u) with your Matlab program newton.m.

Answers and solutions

1.1. Use Matlab to check your answers.

ab = 5, Ab = [\begin{matrix}5 \\ 11 \end{matrix}] , Aa = not defined.

1.2. (a)

f′(x) = [\begin{matrix}cos(x_{1}) & − sin(x_{2}) \\ − sin(x_{1}) & cos(x_{2}) \end{matrix}] , L(x) = f(x̄) + f′(x̄)(x − x̄) = [\begin{matrix}1 \\ 1 \end{matrix}] + [\begin{matrix}1 & 0 \\ 0 & 1 \end{matrix}] [\begin{matrix}x_{1} \\ x_{2} \end{matrix}] .

(b)

f′(x) = [\begin{matrix}0 & 0 \\ 1 & 0 \\ e^{x_{2}} & x_{1}e^{x_{2}} \end{matrix} ] , L(x) = f(x̄) + f′(x̄)(x − x̄) = [\begin{matrix}1 \\ 2 \\ 1 + e \end{matrix}] + [\begin{matrix}0 & 0 \\ 1 & 0 \\ e & e \end{matrix}] [\begin{matrix}x_{1} − 1 \\ x_{2} − 1 \end{matrix}] .

1.3. (a)

∇f(x) = [−e^{−x_{1}} sin(x_{2}), e^{−x_{1}} cos(x_{2})] ,

L(x) = f(x̄) + f′(x̄)(x − x̄) = 0 + [0 1] [\begin{matrix}x_{1} \\ x_{2} \end{matrix}] = x_{2}.

(b) ∇f(x) = [2x_{1} 2x_{3} 2x_{3}] ,

L(x) = f(x̄) + f′(x̄)(x − x̄) = 3 + [2 2 2] [\begin{matrix}x_{1} − 1 \\ x_{2} − 1 \\ x_{3} − 1 \end{matrix}] = −3 + 2x_{1} + 2x_{2} + 2x_{3}.

1.4. (a)

f′(t) = [\begin{matrix}− sin(t) \\ cos(t) \end{matrix}] ,

L(t) = f(t̄) + f′(t̄)(t − t̄) = [\begin{matrix}0 \\ 1 \end{matrix}] + [\begin{matrix}−1 \\ 0 \end{matrix}] (t − π/2).

(b)

f′(t) = [\begin{matrix}1 \\ 2t \end{matrix}] ,

L(t) = f(t̄) + f′(t̄)(t − t̄) = [\begin{matrix}0 \\ 1 \end{matrix}] + [\begin{matrix}1 \\ 0 \end{matrix}] t = [\begin{matrix}t \\ 1 \end{matrix}] .

1.5.

f′(x) = [\begin{matrix}1 − x_{2} & −x_{1} \\ x_{2} & −1 + x_{2} \end{matrix}] ,

L(x) = f(x̄) + f′(x̄)(x − x̄) = [\begin{matrix}0 \\ 0 \end{matrix}] + [\begin{matrix}0 & −1 \\ 1 & 0 \end{matrix}] [\begin{matrix}x_{1} − 1 \\ x_{2} − 1 \end{matrix}]

1.6. (a) The solutions are given by

f(u) = [\begin{matrix}u_{2}(1 − u_{1}^{2}) \\ 2 − u_{1}u_{2} \end{matrix}] = [\begin{matrix}0 \\ 0 \end{matrix}] .

We find two solutions ū = [\begin{matrix}1 \\ 2 \end{matrix}] and ū = [\begin{matrix}−1 \\ −2 \end{matrix}].

(b) The Jacobian is

Df(u) = [\begin{matrix}−2u_{1}u_{2} & 1 − u_{1}^{2} \\ −u_{2} & −u_{1} \end{matrix}] .

(c) The first step of Newton’s method:

evaluate  A = Df(1, 1) = [\begin{matrix}−2 & 0 \\ −1 & −1 \end{matrix}] and b = −f(1, 1) = [\begin{matrix}0 \\ −1 \end{matrix}]
solve  Ah = b, [\begin{matrix}−2 & 0 \\ −1 & −1 \end{matrix}] [\begin{matrix}h_{1} \\ h_{2} \end{matrix}] = [\begin{matrix}0 \\ −1 \end{matrix}]

\{ \begin{matrix}− 2h_{1} = 0, \\ − h_{1} − h_{2} = −1, \end{matrix} h = [\begin{matrix}0 \\ 1 \end{matrix}]

update  u^{(1)} = u^{(0)} + h = [\begin{matrix}1 \\ 1 \end{matrix}] + [\begin{matrix}0 \\ 1 \end{matrix}] = [\begin{matrix}1 \\ 2 \end{matrix}] = ū

Bingo! We found one of the solutions. 1.7. (a) The solutions are given by

f(u) = [\begin{matrix}u_{1}(1 − u_{2}) \\ u_{2}(1 − u_{1}) \end{matrix}] = [\begin{matrix}0 \\ 0 \end{matrix}] .

We find two solutions ū = [\begin{matrix}0 \\ 0 \end{matrix}] and ū = [\begin{matrix}1 \\ 1 \end{matrix}].

(b) The Jacobian is

Df(u) = [\begin{matrix}1 − u_{2} & −u_{1} \\ −u_{2} & 1 − u_{1} \end{matrix}] .

(c) The first step of Newton’s method:

evaluate  A = Df(2, 2) = [\begin{matrix}−1 & −2 \\ −2 & −1 \end{matrix}] and b = −f(2, 2) = [\begin{matrix}2 \\ 2 \end{matrix}]
solve  Ah = b, [\begin{matrix}−1 & −2 \\ −2 & −1 \end{matrix}] [\begin{matrix}h_{1} \\ h_{2} \end{matrix}] = [\begin{matrix}2 \\ 2 \end{matrix}] ,

\{ \begin{matrix}− h_{1} − 2h_{2} = 2, \\ − 2h_{1} − h_{2} = 2, \end{matrix} h = [\begin{matrix}−2/3 \\ −2/3 \end{matrix}]

update  u^{(1)} = u^{(0)} + h = [\begin{matrix}2 \\ 2 \end{matrix}] + [\begin{matrix}−2/3 \\ −2/3 \end{matrix}] = [\begin{matrix}4/3 \\ 4/3 \end{matrix}]

Getting closer to one of the solutions ū!

2006-11-02 /stig