26.3 Newton’s Method

26.3.1 Derivation

In general, Newton’s method uses the same vector notation as already introduced for linear systems. In this case we find a set of n nonlinear equations for n independent unknowns according to Eq. 26.1. We now consider our unknowns as vector $\overrightarrow{x}$ and our functions as a vector function according to

$\overrightarrow{x}=\left(\begin{array}{l}{x}_0\hfill \\ x_1\hfill \\ x_{{}_2}\hfill \\ \cdots \hfill \\ x_n\hfill \end{array}\right),\overrightarrow{f}=\left(\begin{array}{l}{f}_0\hfill \\ f_1\hfill \\ f_{{}_2}\hfill \\ \cdots \hfill \\ f_n\hfill \end{array}\right)$

Please note that in contrast to the matrix notation that we used in section 25.2.1.1, we cannot simply create a matrix here because $\overrightarrow{f}$ is a function vector, i.e., a vector with a function for each entry. This vector is not (as we are used to working with) a vector with fixed numbers as entries. We can now rewrite our original system of equations as

$\overrightarrow{f}\left(\overrightarrow{x}\right)=\overrightarrow{0}$

Because $\overrightarrow{f}$ is given, we need to find a vector $\overrightarrow{x}$ to which we apply the vector function and obtain in the zero vector $\overrightarrow{0}$, which is a vector with a 0 in all vector entries.

26.3.2 Finite Difference

One-Dimensional Case. In vector notation the finite difference is given as

$\Delta {\overrightarrow{x}}_n={\overrightarrow{x}}_{n+1}-\overrightarrow{x}{}_n$

We also know that for the linear case we can derive the value of the function f at x_n+1 by linearization. For this we take the derivative of the function at x_n, i.e., ${\frac{df}{dx}|}_{x_n}$ and multiply it by ∆x_n given by x_n+1− x_n

$f\left(x_{n+1}\right)=f\left(x_n\right)+{\frac{df}{dx}|}_{x_n}\left(x_{n+1}-x_n\right)$

Please note that this equation stems directly from the definition of the derivative (see Fig. 3.2a). For ∆x_n → 0 we will obtain the definition of the derivative (see Eq. 3.9).

n-Dimensional Case. Because we have an n-dimensional problem here, we need to take into account the changes of each individual equation f_i of our function vector with respect to each variable x_j. We therefore require the Jacobian matrix of the function vector $\overrightarrow{f}$. For consistency we will denote the matrix as D, finding

$\text{D}=\left(\begin{array}{cccc}\frac{\partial f_1}{\partial x_1}& \frac{\partial f_1}{\partial x_2}& \dots & \frac{\partial f_1}{\partial x_n}\\ \frac{\partial f_2}{\partial x_1}& \frac{\partial f_2}{\partial x_2}& \dots & \frac{\partial f_2}{\partial x_n}\\ \dots & \dots & \ddots & \dots \\ \frac{\partial f_n}{\partial x_1}& \frac{\partial f_n}{\partial x_2}& \dots & \frac{\partial f_n}{\partial x_n}\end{array}\right)$

We already encountered the Jacobian matrix during coordinate system transformations (see section 7.5). This matrix effectively represents the changes of our function vector $\overrightarrow{f}$ with respect to all independent variables. We will then apply the same linearization as for the one-dimensional case (see Eq. 26.9) in vector notation, finding

$\begin{array}{ll}\overrightarrow{f}\left({\overrightarrow{x}}_{n+1}\right)\hfill & =\overrightarrow{f}\left({\overrightarrow{x}}_n\right)+\overrightarrow{\nabla }{\overrightarrow{f}|}_{{\overrightarrow{x}}_n}\left({\overrightarrow{x}}_{n+1}-{\overrightarrow{x}}_n\right)\hfill \\ \hfill & =\overrightarrow{f}\left({\overrightarrow{x}}_n\right)+\text{D}\left({\overrightarrow{x}}_n\right)\left({\overrightarrow{x}}_{n+1}-{\overrightarrow{x}}_n\right)\hfill \end{array}$

Eq. 26.11 gives the analytical correlation between the solution ${\overrightarrow{x}}_n$ and the (improved) subsequent solution ${\overrightarrow{x}}_{n+1}$.

Solving for${\overrightarrow{x}}_{n+1}$. In Newton’s method we are seeking the solution $\overrightarrow{x}$ that makes $\overrightarrow{f}\left(\overrightarrow{x}\right)=\overrightarrow{0}$ We can therefore rewrite Eq. 26.11 as

$\overrightarrow{f}\left({\overrightarrow{x}}_{n+1}\right)=\overrightarrow{f}\left({\overrightarrow{x}}_n\right)+\text{D}\left({\overrightarrow{x}}_n\right)\left({\overrightarrow{x}}_{n+1}-{\overrightarrow{x}}_n\right)=\overrightarrow{0}$

which we solve for ${\overrightarrow{x}}_{n+1}$ finding

$\overrightarrow{f}\left({\overrightarrow{x}}_n\right)+\text{D}\left({\overrightarrow{x}}_n\right)\left({\overrightarrow{x}}_{n+1}-{\overrightarrow{x}}_n\right)=\overrightarrow{0}$

$\begin{array}{ll}\text{D}\left({\overrightarrow{x}}_n\right)\left({\overrightarrow{x}}_{n+1}-{\overrightarrow{x}}_n\right)\hfill & =-\overrightarrow{f}\left({\overrightarrow{x}}_n\right)\hfill \\ {\text{D}}^{-1}\left({\overrightarrow{x}}_n\right)\text{D}\left({\overrightarrow{x}}_n\right)\left({\overrightarrow{x}}_{n+1}-{\overrightarrow{x}}_n\right)\hfill & =-{\text{D}}^{-1}\left({\overrightarrow{x}}_n\right)\overrightarrow{f}\left({\overrightarrow{x}}_n\right)\hfill \\ {\overrightarrow{x}}_{n+1}\hfill & ={\overrightarrow{x}}_n-{\text{D}}^{-1}\left({\overrightarrow{x}}_n\right)\overrightarrow{f}\left({\overrightarrow{x}}_n\right)\hfill \end{array}$

Eq. 26.13 will give the exact solution to our problem. Unfortunately, it requires the inverse of the matrix D, which we are not even sure exists. We therefore need to circumvent this problem.

26.3.3 Simplification of the Newton’s Method: Quasi-Newton Method

A number of methods essentially simplify Eq. 26.13 such that the inverse matrix D⁻¹ that is required for evaluating D⁻¹ (${\overrightarrow{x}}_n$) need not be calculated. All of these methods are referred to as quasi-Newton methods because they use the underlying concept of the Newton method but simplify this step. Several noteworthy methods are commonly used, but we will focus on the most simple one. For this we use Eq. 26.12 instead of Eq. 26.13. We rewrite this equation to

$\text{D}\left({\overrightarrow{x}}_n\right)\left({\overrightarrow{x}}_{n+1}-{\overrightarrow{x}}_n\right)=\text{D}\left({\overrightarrow{x}}_n\right)\Delta {\overrightarrow{x}}_n=-\overrightarrow{f}\left({\overrightarrow{x}}_n\right)$

Since we can calculate D and we know ${\overrightarrow{x}}_n$, we obtain a system of n equations from Eq. 26.14 that we solve for n unknowns. The matrix D will contain all nonlinear terms of our problem, which will be converted to numbers once we apply the (known) vector ${\overrightarrow{x}}_n$ to it. We are then left with a set of linear equations that we can solve using one of the methods outlined in section 25.2. Solving these equations will give us the vector ${\overrightarrow{x}}_{n+1}$, which we can apply iteratively to Eq. 26.14. We terminate the process once the change in the solutions $\Delta {\overrightarrow{x}}_n={\overrightarrow{x}}_{n+1}-{\overrightarrow{x}}_n$ is below a certain threshold.

26.3.4 Applying the Newton’s Method to Our Example

We will now demonstrate the Newton’s method for the example of the position determination in LORAN. For this we rewrite Eq. 26.7 to a function vector

$\overrightarrow{f}=\left(\begin{array}{c}\sqrt{{\left(x+200\right)}^2+y^2}-\sqrt{{\left(x-200\right)}^2+y^2}-96.84\\ \sqrt{{\left(x+200\right)}^2+y^2}-\sqrt{{\left(x-220\right)}^2+{\left(y+100\right)}^2}-77.75\end{array}\right)$

Obviously the solution vector is

$\overrightarrow{x}=\left(\begin{array}{l}x\hfill \\ y\hfill \end{array}\right)$

Calculating the Matrix D. We also require the matrix D for Eq. 26.14. This we find according to Eq. 26.10

$\text{D}=\left(\begin{array}{ll}\frac{\partial f_x}{\partial x}\hfill & \frac{\partial f_x}{\partial y}\hfill \\ \frac{\partial f_y}{\partial x}\hfill & \frac{\partial f_y}{\partial y}\hfill \end{array}\right)$

which we carry out to find

$D=\left(\begin{array}{l}\begin{array}{c}\frac{x+200}{\sqrt{{\left(x-200\right)}^2+y^2}}-\frac{x-200}{\sqrt{{\left(x-200\right)}^2+y^2}}\frac{y}{\sqrt{{\left(x+200\right)}^2+y^2}}-\frac{y}{\sqrt{{\left(x-200\right)}^2+y^2}}\\ \end{array}\\ \frac{x+200}{\sqrt{{\left(x-200\right)}^2+{\left(y+200\right)}^2}}-\frac{x-200}{\sqrt{{\left(x-200\right)}^2+{\left(y+100\right)}^2}}\frac{y}{\sqrt{{\left(x+200\right)}^2+y^2}-\frac{y+200}{\sqrt{{\left(x+220\right)}^2+{\left(y+100\right)}^2}}}\end{array}\right)$

First Guess. We are now ready to start the numerical solving process. We first require an initial guess for $\overrightarrow{x}$. For the sake of simplicity, we will assume

${\overrightarrow{x}}_0=\left(\begin{array}{c}0\\ 0\end{array}\right)$

Next we find D (${\overrightarrow{x}}_0$) for the first estimate ${\overrightarrow{x}}_0$ which is

$\text{D}\left({\overrightarrow{x}}_0\right)=\left(\begin{array}{cc}2& 0\\ 1.91& -0.41\end{array}\right)$

Next we find $\overrightarrow{f}\left({\overrightarrow{x}}_0\right)$ to be

$\overrightarrow{f}\left({\overrightarrow{x}}_0\right)=\left(\begin{array}{c}-96.84\\ -119.41\end{array}\right)$

Using Eq. 26.14, we carry out the matrix/vector multiplication

$\left(\begin{array}{cc}2& 0\\ 1.91& -0.41\end{array}\right)\left({\overrightarrow{x}}_1-\overrightarrow{0}\right)=\left(\begin{array}{cc}-1& 0\\ -0.91& 0.41\end{array}\right)\left(\begin{array}{c}{x}_1\\ y_1\end{array}\right)=\overrightarrow{f}\left({\overrightarrow{x}}_n\right)=\left(\begin{array}{c}-96.84\\ -119.41\end{array}\right)$

from which we derive the two equations

$2x_1=96.84$

$1.91x_1-0.41y_1=119.41$

As you can see, we now have a linear system of equations that can be solved conveniently (see section 25.2). This system is so easy that we can actually just write the solution for ${\overrightarrow{x}}_1$ as

${\overrightarrow{x}}_1=\left(\begin{array}{c}48.42\\ -65.68\end{array}\right)$

Obviously, this result is not entirely correct. If we look at Fig. 26.2, we see that the position we are seeking is actually (50, −50). But we have definitely approached the solution starting from our initial guess (0, 0).

Iterative Scheme. Now that we have the improved estimate for the position of our vessel, we will repeat this process in order to obtain an even better approximation ${\overrightarrow{x}}_2$. However, it makes little sense to do this manually. The strength of numerical solutions is that they can be solved significantly faster by algorithms. Therefore we will implement a small solver in Maple that performs this task for us.

26.4.1 Helper Function: Matrix/Vector Multiplication

We will begin by implementing three helper functions. The first function will perform a matrix/vector multiplication, which is required for the Newton method. We carried out this multiplication manually in Eq. 26.17.

Just as a quick reminder to matrix/matrix calculus: A n×m matrix can be multiplied only by a m×o matrix. The first number indicates the number of lines; the second number indicates the number of columns. This means that a matrix A can be multiplied only by a second matrix B if the number of lines of the first matrix and the number of columns of the second matrix are identical. The result of this multiplication is a n × o matrix C given by

$\begin{array}{l}\left(\begin{array}{cccc}{A}_{1,1}& A_{1,2}& \dots & A_{1,n}\\ A_{2,1}& A_{2,2}& \dots & A_{2,m}\\ ⋮& ⋮& \ddots & ⋮\\ A_{n,1}& A_{n,2}& \dots & A_{n,m}\end{array}\right).\left(\begin{array}{cccc}{B}_{1,1}& B_{1,2}& \dots & B_{1,o}\\ B_{2,1}& B_{2,2}& \dots & B_{2,o}\\ ⋮& ⋮& \ddots & ⋮\\ B_{m,1}& B_{m,2}& \dots & B_{m,o}\end{array}\right)\\ =\left(\begin{array}{cccc}{A}_{1,1}{B}_{1,1}+\dots +A_{1,m}{B}_{m,1}& A_{1,1}{B}_{1,2}+\dots +A_{1,m}{B}_{m,2}& \dots & A_{1,1}{B}_{1,o}+\dots +A_{1,m}{B}_{m,o}\\ A_{2,1}{B}_{1,1}+\dots +A_{2,m}{B}_{m,1}& A_{2,1}{B}_{1,2}+\dots +A_{2,m}{B}_{m,2}& \dots & A_{2,1}{B}_{1,o}+\dots +A_{2,m}{B}_{m,o}\\ ⋮& ⋮& \ddots & ⋮\\ A_{n,1}{B}_{1,1}+\dots +A_{n,m}{B}_{m,1}& A_{n,1}{B}_{1,2}+\dots +A_{n,m}{B}_{m,2}& ⋮& A_{n,1}{B}_{1,o}+\dots +A_{n,m}{B}_{m,o}\end{array}\right)\end{array}$

In this multiplication the entry (i, k) of matrix C will be formed by summing the products of i^th line of matrix A with the k^th column of matrix B, respectively. In our example, we have two equations (Eq. 26.7) and two unknowns. Therefore the matrix we are dealing with is of type 2 × 2. The second matrix is the vector ${\overrightarrow{x}}_n$, i.e., a 2 × 1 matrix (Eq. 26.17). We implement the matrix/vector multiplication as shown in listing 26.1.

1  DoMatrixVectorMult:=proc(D::Matrix,v::Matrix,n) local i,j,ret:
2    ret:=Matrix(n,1):
3    for i to n do
4      for j to n do:
5        ret[i,1]:=D[i,j]*v[j,1]+ret[i,1]:
6      end do:
7     end do:
8     return ret;
9  end proc:

The function is termed DoMatrixVectorMult, and it expects a quadratic matrix with m = n number of lines and columns given by the parameter n. The input matrix D is given as a matrix as is the input vector v. The matrix is expected to be a n × n matrix and the vector to be a n × 1 matrix. The output ret therefore is a n × 1 matrix, i.e., a vector. The function then iterates over the lines i and sums the columns j. Please note that Maple has built-in functions for matrix multiplication, but the version used in listing 26.1 is so straightforward that it can be ported to any programming language.

26.4.2 Helper Function: Matrix Evaluation

The second helper function we will implement will evaluate a matrix function when supplied with a given value vector. We require this function for evaluating the matrix D using the value vector ${\overrightarrow{x}}_n$ as we did for finding D (${\overrightarrow{x}}_0$) in Eq. 26.16.

The function is shown in listing 26.2. It expects the n × n matrix function D in the form of a matrix of functions of n independent variables. The parameter values holds the numeric values to insert into the matrix function. The return value ret will of course be an n × n matrix again. The function then merely calls the individual entry (i, j) of D providing the values given for the respective independent variable. This is done by the respective Maple command seq, which yields a sequence, i.e., a comma-separated list of the values. As you can see in the way line 5 is written, the code simply calls the entry (i, j) of the matrix using the values provided. Because the entry of the matrix is expected to be a function, this notation is equivalent to calling a function in Maple with the values for the variables in bracket.

1  DoEvaluateMatrix:=proc(D::Matrix,values::Matrix,n) local i,j,ret:
2    ret:=Matrix(n,n):
3    for i to n do
4     for j to n do:
5      ret[i,j]:=D[i,j](seq(values[k,1],k=1..n)):
6     end do:
7    end do:
8  return ret;
9  d proc:

26.4.3 Helper Function: Vector Evaluation

The last helper function we require is very similar to the one we just implemented. It evaluates a function vector using a value vector. The function is shown in listing 26.3. As you can see, the only difference is that the helper function expects the input vector v to be a n × 1 matrix. Therefore we do not have to sum over the columns j, which makes the listing a bit more compact. Except for this detail the implementations are identical.

1  DoEvaluateVector:=proc(v::Matrix,values::Matrix,n) local i,ret:
2  ret:=Matrix(n,1):
3    for i to n do
4      ret[i,1]:=v[i,1](seq(values[k,1],k=1..n)):
5    end do:
6    return ret;
7  end proc:

26.4.4 Listing for Newton’s Method

Having implemented the preceding three helper functions, we are now ready to implement Newton’s method.

The code is shown in listing 26.4. It begins as usual by restarting the Maple server and including Core.txt (see line 1). We then implement the three helper functions DoMatrixVectorMult, DoEvaluateMatrix, and DoEvaluateVector whose listings we have already discussed. The implementation of these functions is not repeated in this listing; please refer to listing 26.1, listing 26.2, and listing 26.3 for details on these functions. The function DoNonlinearNewton implements Newton’s method. The function expects a list of nonlinear equations equations that needs to be solved for a list of independent variables given in variables using the initial values initial. The function also expects a maximum number of steps maxsteps to perform.

1  restart;read("Core.txt"):
2
3  #matrix/vector multiplication
4  DoMatrixVectorMult:=proc(D::Matrix,v::Matrix,n) local i,j,ret:
5
6  #function matrix evaluation using a value vector
7  DoEvaluateMatrix:=proc(D::Matrix,values::Matrix,n) local i,j,ret:
8
9  #function vector evaluation using a value vector
10  DoEvaluateVector:=proc(v::Matrix,values::Matrix,n) local i,ret:
11
12  #this function solves a non-linear system of equations using Newton's method
13  DoNonlinearNewton:=proc(equations::list,variables::list,initial::list,maxsteps) local n,i,j,k,s,xn, D,Dxn,f,fxn,dx,lhs,sol:
14    #check for consistency
15    n:=numelems(equations):
16    if numelems(variables)<>n or numelems(initial)<> n then
17      print ("Require same number of equations, variables and initial values"):
18      return:
19    end if:
20    #construct D
21    D:=Matrix(n,n):
22    for i to n do
23      for j to n do
24       D[i,j]:=unapply(diff(equations[i],variables[j]),seq(variables[k],k=1..n));
25      end do:
26    end do:
27    #construct the function vector f
28    f:=Matrix(n,1):
29    for i to n do
30      f[i,1]:=unapply(equations[i],seq(variables[k],k=1..n)):
31    end do:
32    #this holds the current value for xn - prime with the initial values
33    xn:=Matrix(n,1):
34    for i to n do
35      xn[i,1]:=initial[i]:
36    end do:
37    #differential vector
38    dx:=Matrix([seq([Delta[k]],k=1..n)]):
39    #now we iterate until maxsteps
40    for s to maxsteps do
41      #output to the user
42      print(sprintf(`step: %d, %s`,s,cat(seq(sprintf("%a=%a ",variables[k],xn[k,1]),k=1..n)))):
43      #calculate f(xn) and D(xn)
44      fxn:=DoEvaluateVector(f,xn,n):
45      Dxn:=DoEvaluateMatrix(D,xn,n):
46      #perform the matrix vector multiplication
47      lhs:=DoMatrixVectorMult(Dxn,dx,n):
48      #solve the new system of equations
49      try
50        sol:=solve({seq(lhs[k,1]=-fxn[k,1],k=1..n)},{seq(dx[k,1],k=1..n)}):
51        #adapt xn
52        for i to n do
53          xn[i,1]:=evalf(xn[i,1]+rhs(sol[i])):
54        end do:
55      catch:
56        print("no solution to linear systemml: ",seq(lhs[k,1]=-fxn[k,1],k=1..n));
57        return;
58     end try;
59   end do:
60  end proc:

In line 15 the function first checks for consistency. The number of equations n must match the number of unknown variables given as a list in variables and the number of initial values provided in initial. If this prerequisite is not met, the function quits by printing an error message. In line 21 the matrix D is created as an n × n matrix. The entry (i, j) of D is created by finding the derivative of the i^th equation given with respect to the j^th independent variable. The result of this diff operation is turned into a function of the independent variables (again given as a sequence) by performing an unapply operation.

In pretty much the same way, we then create the function vector f by performing an unapply operation to the i^th equation of equations line 30. In the next step, the current estimation ${\overrightarrow{x}}_n$ for $\overrightarrow{x}$ is created as a n × 1 vector (see line 33). The first assumption for ${\overrightarrow{x}}_n$ is the initial values given to the function as the parameter initial. The last step before starting the iteration is the definition of the differential vector dx according to Eq. 26.14 in line 38. Please note that we do not directly derive ${\overrightarrow{x}}_{n+1}$ here, but rather the difference vector $\Delta {\overrightarrow{x}}_n$ which is why we later have to add this delta value to the last estimate of the solution ${\overrightarrow{x}}_n$ (see line 53). dx is a n × 1 matrix that is filled with the values Delta[1], Delta[2], etc. They will be the independent variables for which we solve the system of linear equations we obtain.

The loop incorporates the iterative scheme of Newton’s method. We iterate s to a total number of maxsteps (see line 40). We first output the current step s and the values for the independent variables. This happens in line 42. This lines uses print to output a string created by sprintf. This function creates formatted string output. The function cat joins strings. Here it outputs the individual independent variables variables as a sequence using the form “name=current value”. This helps us keep track of the current value of our approximations.

The following lines implement the scheme we already performed manually. First, the function vector f and the matrix D are evaluated for the current value vector xn using the two helper functions DoEvaluateVector and DoEvaluateMatrix we previously implemented (line 44 and line 45). We then perform a matrix/vector multiplication using the helper function DoMatrixVectorMult we also implemented previously (see line 47). For this we use the evaluated matrix D (${\overrightarrow{x}}_n$) and the differential vector $\Delta {\overrightarrow{x}}_n$. In our manual example, we performed this multiplication in Eq. 26.17. We were then left with a system of linear equations in Eq. 26.18. The same thing happens at this point. We now need to solve this system. This happens in the code lines starting from line 49. In line 50 we call solve to make Maple solve the system of linear equations given by Eq. 26.14. Because we have a total of n equations, we provide these in form of a sequence that Maple should solve for the unknown variables dx[1], dx[2], etc. We store the solution in sol from which we extract the individual values for the variables using the function rhs that extracts the right-hand side of an equation. We then correct our original assumption for the k^th variable by the correction dx[k]. As stated this procedure is slightly different from the first manual step we performed in Eq. 26.17 where we used the vector ${\overrightarrow{x}}_{n+1}-{\overrightarrow{x}}_n$ instead of $\Delta {\overrightarrow{x}}_n$. This is why we need to correct ${\overrightarrow{x}}_n$ by $\Delta {\overrightarrow{x}}_n$ in line 53. Having corrected the assumption, this process is repeated for a total of maxsteps iteration.

It is also possible to check whether the value dx is smaller than a given threshold, i.e., if the solution can be considered sufficiently exact. At this point the function can be aborted.

You may have noted that the call to solve is embedded in a try/catch block (see line 49). This is due to the fact that the algorithm may run into equations that cannot actually be solved. We will look at an example of such a call in just a moment. If this happens, the algorithm obviously cannot continue and aborts by providing an error message to the user.

26.4.5 Using the Solver

We are now ready to apply the solver we just implemented to solve Eq. 26.7, which we obtained from the LORAN positioning system. Listing 26.5 shows how to apply the function DoNonlinearNewton to Eq. 26.7. We provide the two equations and the two unknowns x, y in list form to the function. We also need to provide an initial estimation.

1 DoNonlinearNewton([sqrt((x+200)^2+y^2)-sqrt((x-200)^2+y^2)-96.84,sqrt((x+200)^2+y^2)-sqrt((x-220) ^2+(y+100)^2)-77.75],[x,y],[0,0],20);
2 DoNonlinearNewton([sqrt((x+200)^2+y^2)-sqrt((x-200)^2+y^2)-96.83,sqrt((x+200)^2+y^2)-sqrt((x-220) ^2+(y+100)^2)-77.75],[x,y],[-100,600],20);

Here we will see a good example of why it is important to provide reasonable and often initial values which are close to the actual solution. We start with line 1 where we assume a first estimate of (0, 0). If this code is run, Maple will output the current step and the current values for x and y. The values are given in Tab. 26.1.

As you can see, the first run assumes our initial values (0, 0). At the end of the first iteration, we obtain the corrected values we obtained while performing the algorithm manually (see Eq. 26.19). The values are not exactly the same because during the manual derivation, we considered only 2 decimals, whereas our numerical solver operates at machine-precision. Interestingly, at the end of the second iteration, the results have almost converged to the correct solution. From Fig. 26.2 we know that the vessel is located at (50, −50). After the fifth iteration, the solver has reached a stable solution, and the output does not change. This is the point where the algorithm could be aborted. If higher precision is required, the general settings of Maple must be changed, i.e., the number of decimal places to which a numerical value is evaluated must be increased.

As can be seen, the implemented solvers find a solution to this set of nonlinear equations pretty fast. In the next step, we will try to solve the same equations using a different set of initial values.

A Common Case: A Solver Finds No Solution. For the second example, we use the initial values (−100, 600) (see line 2). As you can see from Fig. 26.1, this first estimate is not too far off the actual solution. If the ship was to determine its position for the first time, there would be no way of determining whether (0, 0) or (−100, 600) would be the better initial assumption. Running this line of code will output the values shown in Tab. 26.2.

As you can see, the solver will perform only 11 steps during which the solution values grow continually. At the end of step 11, the solver will abort with an error message saying

no solution to linear systemml: 0 = 96.83

Obviously, this system of equations cannot be solved as 0 6 = 96.83. This is a good example of what happens when a bad initial value is supplied to a solver and why we had to account for similar problems by employing the try/catch construction in line 49. Numerical solvers will eventually choke on some bad initial values throwing comparable error messages.

When employing a numerical solver, it is always a good idea to learn as much about the problem at hand as possible in order to provide good initial values. If this is not possible, the initial values must be estimated, and a trial-and-error approach must be used. For this initial values are varied until the solver converges. However, even then it is important to validate the results obtained in order to ensure that they are physically meaningful. In the case of the LORAN system, we may come across locations that give rise to more than one solution. The vessel then will need to estimate which of the solutions obtained are correct. This can be done by picking up more signals from other beacons and checking whether the solutions obtained have one common location. A second approach may be to track the location over time, in which case all solutions that are too far away from the last known location can be discarded.

Physically Nonmeaningful Results. We encountered an example where physically nonmeaningful results can be obtained from a numerical solver while studying the pendant drop equations (see section 24.5). We found that several of the results that we found numerically are physically impossible because the drop volumes do not increase and there is a discontinuity in the drop volumes (see Fig. 24.9b). In all cases the solutions found are mathematically correct but not physically correct. Therefore it is always important to validate a numerically obtained solution for its physical significance.