# The Newton idea for systems¶

Let us use Newton’s method on the system defined by the function

```
using LinearAlgebra
```

```
function nlvalue(x)
return [ exp(x[2]-x[1]) - 2,
x[1]*x[2] + x[3],
x[2]*x[3] + x[1]^2 - x[2]
]
end
```

```
nlvalue (generic function with 1 method)
```

Here is a function that computes the Jacobian matrix.

```
function nljac(x)
J = zeros(3,3)
J[1,:] = [-exp(x[2]-x[1]), exp(x[2]-x[1]), 0]
J[2,:] = [x[2], x[1], 1]
J[3,:] = [2*x[1], x[3]-1, x[2]]
return J
end
```

```
nljac (generic function with 1 method)
```

(These functions could be written as separate files, or embedded within another function as we have done here.) Our initial guess at a root is the origin.

```
x = [0.,0.,0.]
```

```
3-element Array{Float64,1}:
0.0
0.0
0.0
```

We need the value (residual) of the nonlinear function, and its Jacobian, at this value for \(\mathbf{x}\).

```
f = nlvalue(x)
```

```
3-element Array{Float64,1}:
-1.0
0.0
0.0
```

```
J = nljac(x)
```

```
3×3 Array{Float64,2}:
-1.0 1.0 0.0
0.0 0.0 1.0
0.0 -1.0 0.0
```

We solve for the Newton step and compute the new estimate.

```
s = -(J\f)
newx = x + s
```

```
3-element Array{Float64,1}:
-1.0
0.0
0.0
```

Here is the new residual.

```
nlvalue(newx)
```

```
3-element Array{Float64,1}:
0.7182818284590451
0.0
1.0
```

We don’t seem to be especially close to a root. Let’s iterate a few more times. To do that, we start keeping the root estimates as elements of a vector (i.e., a vector of vectors). Same for the residuals.

```
x = [x,newx]
```

```
2-element Array{Array{Float64,1},1}:
[0.0, 0.0, 0.0]
[-1.0, 0.0, 0.0]
```

```
y = [f,nlvalue(newx)]
```

```
2-element Array{Array{Float64,1},1}:
[-1.0, 0.0, 0.0]
[0.7182818284590451, 0.0, 1.0]
```

```
for n = 2:7
f = nlvalue(x[n])
push!(y,f)
s = -nljac(x[n])\f
push!(x,x[n]+s)
end
```

We find the infinity norms of the residuals.

```
residual_norm = [ norm(f,Inf) for f in y ] # max in dimension 1
```

```
8-element Array{Float64,1}:
1.0
1.0
1.0
0.20229273399879752
0.010252098544050403
2.1556384274173945e-5
1.9899884518004285e-10
1.3877787807814457e-17
```

We don’t know an exact answer, so we will take the last computed value as its surrogate.

```
r = x[end]
```

```
3-element Array{Float64,1}:
-0.45803328064126886
0.23511389991867648
0.10768999090411434
```

The following will subtract r from every column of x.

```
error = [ norm(xk-r) for xk in x[1:end-1] ]
```

```
7-element Array{Float64,1}:
0.525994454581891
0.6005027934725727
0.15184335534811869
0.010200335585503865
2.4987473439459857e-5
1.5326877679492286e-10
1.3877787807814457e-17
```

Now we take infinity norms and check for quadratic convergence.

```
ratios = [ log(error[k+1]) / log(error[k]) for k in 1:length(error)-1 ]
```

```
6-element Array{Float64,1}:
0.7937993447128695
3.695980885448303
2.4326597889977144
2.3110932374368462
2.1325411123294264
1.7176219080468489
```

For a brief time, we see ratios around 2, but then the limitation of double precision makes it impossible for the doubling to continue.