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):

enter image description here

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,

  • $\begingroup$ @Amzoti thanks for the reply, I don't, I only have B values (water depth in m), and I need to get C values based on B. I tried a simple set of B values (500, 1000, 2000) to get C, but replacing those values in any of that equations didn't seem to work. To solve for Newtons method I would have an approximate value, in this case that would be 600 (because Tw equation and Tst_100 intersect around that value, for instance). $\endgroup$ – Gery May 19 '14 at 4:54
  • $\begingroup$ According to the paper, $C$ represents the thickness of the GHSZ, or Gas Hydrate Stability Zone, which is presumably some constant. $\endgroup$ – Mark McClure May 19 '14 at 11:47
  • $\begingroup$ It's possible to solve for $C$ explicitly using the Lambert W function. Would you be interested in such an answer? $\endgroup$ – Antonio Vargas May 20 '14 at 16:53
  • $\begingroup$ @AntonioVargas Im interested, it would be cool Antonio, at the beginning I was wanting to find an answer about how the 3 equations were solved with the Newtons method, but now I just want to get C based on B. If you can also tell me why Newtons method doesnt apply here, it would be much nicer. $\endgroup$ – Gery May 20 '14 at 17:01
  • 1
    $\begingroup$ As for the explicit solution, for the equation $$a_1 C - a_2 B^{-a_3} + a_4 \log(C+B) - a_5 = 0,$$ the solution you seek is given by $$C = -B + \frac{a_4}{a_1} W_{-1}\left[\frac{a_1}{a_4} \exp\left(\frac{a_1}{a_4} B + \frac{a_2}{a_4} B^{-a_3} + \frac{a_5}{a_4}\right)\right].$$ Here $W_{-1}$ is the lower real branch of the Lambert W function. This expression was obtained with Mathematica. $\endgroup$ – Antonio Vargas May 20 '14 at 18:20

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$:

enter image description here

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]];
   {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:

enter image description here

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)$.

  • $\begingroup$ Cool, this answer and the formula look great, thanks for that Antonio. I will try the approach in your answer to check if it works also for the other equations. If so, I will mark it as the answer. About the direct equation using the Lambert function, it looks like the solution I need (because I need to directly input B values going from ~100 to ~3000) but I am still confused by that W, do you have further details to get that W? I'd be great if you add that commentary to your answer. $\endgroup$ – Gery May 21 '14 at 4:55
  • $\begingroup$ Are you asking how to calculate $W_{-1}(x)$? If so, what are you doing your calculations in? (Mathematica, MATLAB, Python, etc.) In Mathematica it's implemented as LambertW[-1,x]. $\endgroup$ – Antonio Vargas May 21 '14 at 5:19
  • $\begingroup$ yes exactly, that was the question, so far I've been using awk among other command-line tools in linux, sometimes I use excel, but nothing more elaborated than that. So, Mathematica is new for me, and it seems that I need to buy one license for that (not an option so far), the same for Matlab. Python can be an option, but haven't tried it so far. In fact, I tend to solve my equations with a pen and a paper, quite archaic I know. $\endgroup$ – Gery May 21 '14 at 6:06
  • $\begingroup$ There are many software packages which will calculate it. Take a look at this list from the Lambert W's Wikipedia page. Unfortunately I don't think Excel is one that will do it automatically :). On that list, the free options are PARI/GP, octave, Maxima, and the GSL. $\endgroup$ – Antonio Vargas May 21 '14 at 15:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.