equations solved with Newton's method by finding the zeros of functions?

Question

I found this statement in one paper I read recently:

This problem can be solved by finding the zero of functions:

f(B) = -37.5 × C / 1000 - 373.41 × B^(-0.6269) + 9.1256 × Ln(C + B) - 50.029
f(B) = -37.5 × C / 1000 - 373.41 × B^(-0.6269) + 7.3637 × Ln(C + B) - 34.502
f(B) = -37.5 × C / 1000 - 373.41 × B^(-0.6269) + 7.0724 × Ln(C + B) - 29.970

Newton’s method was used to find the zero of these functions.

Here is the link to the PDF (page 6): http://tao.cgu.org.tw/pdf/v174p815.pdf. The other equations in that paper (the ones that form the above equations) are easy to solve, but not sure how the Newton's method was used to solve the three equations above (that paper said they were solved with that method).

I appreciate if someone could explain me if they indeed were solved with Newton's method. If so, an example using one of these equations would be great.

I would like to add a bit more of information that can be useful to make things a bit clearer. The image below shows how C (curves 1, 2 and 3 in the attached image) varies with B (water depth or seafloor in the image):

Basically, the intersection of the gas hydrate stability curves (in the paper I mentioned these are expressed by the logarithmic curves eg. Tst100) and the geothermal gradient profile (in the paper it is expressed as 37.5*B/1000) provides an estimate of the thickness of the GHSZ. The three equations above should solve that (3 different cases as the attached image). I need to find the values of C based on B.

Any help is grateful, thanks in advance,

Answer 1

I'll describe a rudimentary method for making the plots in Figure 4 of the paper. Since the three equations are so similar I'll just consider the first,

$$ -37.5 C/1000 - 373.41 B^{-0.6269} + 9.1256 \log(C+B) - 50.029 = 0. \tag{1} $$

We'll solve $(1)$ for $C$ numerically as $B$ ranges over the interval $[600,3500]$. For each new application of Newton's method we'll use linear extrapolation to determine the new initial guess.

We need two starting points. Let's take $B = 600$ and define the function

$$ \begin{align} g(C) &= -37.5 C/1000 - 373.41 B^{-0.6269} + 9.1256 \log(C+B) - 50.029 \\ &= -37.5 C/1000 + 9.1256 \log(C + 600) - 56.7986, \end{align} $$

which is the left hand side of $(1)$. Below is a plot of $g(C)$ for $0 \leq C \leq 100$:

It appears to have a root around $C = 70$, so taking $B = 600$ in $(1)$ and using Newton's method with an initial guess of $x_0 = 70$ we obtain the numerical root $C = 68.2908$. We'll call this $C_{600}$.

Similarly we find that, for $B = 601$, equation $(1)$ has a root at $C_{601} = 69.1582$.

Let's calculate our linear extrapolation formula. The line passing through the two points $(n-1,C_{n-1})$ and $(n,C_{n})$ is given by

$$ y = (C_{n}-C_{n-1})(x-n)+C_n, $$

so taking $x = n+1$ we find that our initial guess for $C_{n+1}$ will be

$$ C_{n+1} \approx 2C_n - C_{n-1}. \tag{2} $$

So, for example, our initial guess for using Newton's method with $B = 602$ will be

$$ C_{602} \approx 2C_{601} - C_{600} = 70.0256. $$

Now we calculate $C_n$ iteratively for $n = 602,603,\ldots,3500$. Here is some Mathematica code which will do just that:

pts = {{600, 68.29080925572698`}, {601, 69.15821186254433`}};
lastTwo =.;
guess =.;
For[b = 602, b <= 3500, b++,
  lastTwo = Take[pts, -2];
  guess = 2 lastTwo[[2, 2]] - lastTwo[[1, 2]];
  AppendTo[pts,
   {b, c /.
     FindRoot[-37.5 c/1000 - 373.41 b^(-0.6269) + 9.1256 Log[c + b] - 50.029,
       {c, guess}]}
   ]
  ];
ListPlot[pts, Joined -> True]

Mathematica complains a little because of the loss of precision from using decimals instead of fractions but it still produces the desired plot:

In the Mathematica code, the FindRoot function uses Newton's method with guess as the starting point for the iteration. The variable guess is calculated in the 6th line of the code using formula $(2)$.