Iterative method to solve nonlinear equations

Question

I am trying to solve the following equation (coming from fluid dynamics) using TeX.

1/sqrt(lambda) = 2 * log(Re * sqrt(lambda)) - 0.8

where lambda is the unknown, sqrt means "square root", the log is decimal logarithm, and Re is a given number.

This can't be solved directly... It cannot be written as lambda = ....

How can I improve the code or the approach below to find the zero faster?

I have written code in the fp package and with Bruno Le Floch help a version with is programmed in the LaTeX3 equivalent l3fp.

The approach is to start from lambda small, and increase it until the difference between left- and right- side becomes (close to) zero. The problem is that there is a constant step, so I need to do many steps to find lambda. But I would like the amount of steps to be a minimum.

Here the MWE

\documentclass{article}
\usepackage{expl3}
% Start
\ExplSyntaxOn
\fp_new:N \lambdaa
\fp_new:N \epsa
\fp_new:N \schrittweite
\fp_new:N \zelem
\fp_new:N \hstep
\fp_new:N \RE
\int_new:N \stepsa
\int_new:N \stepsb
% This would end up in the definition of a command if you wish.
\fp_set:Nn \lambdaa { 1e-12 }
\fp_set:Nn \schrittweite { 1e-04 }
\fp_set:Nn \epsa { 0 }
\fp_set:Nn \zelem { 0 }
\fp_set:Nn \hstep { 0.00001 }
\fp_set:Nn \RE { 35800 }
\int_zero:N \stepsa
\int_zero:N \stepsb
%
\fp_set:Nn \differe
  {
    \lambdaa ** -0.5
    - ( 2 * ln(\lambdaa ** 0.5 * \RE) / ln(10) - 0.8 )
  }
\bool_while_do:nn % also \bool_do_while:nn may be useful
  { \fp_compare_p:n { \differe > \zelem } }
  {
    \fp_set:Nn \lambdab { \lambdaa }
    \fp_add:Nn \lambdaa { \hstep } % \lambdaa = \lambdaa + \hstep
    \int_incr:N \stepsb
    \fp_set:Nn \differe
      {
        \lambdaa ** -0.5
        - ( 2 * ln(\lambdaa ** 0.5 * \RE) / ln(10) - 0.8 )
      }
    \typeout{Lambda~zu~klein~[...]~\int_use:N\stepsb,~
      \fp_use:N\lambdab,~\fp_use:N\lambdaa}
    \typeout{diff~=~\fp_use:N\differe} % Note: "~" gives a space.
  }
\bool_while_do:nn
  { \fp_compare_p:n { \differe < \epsa }}
  {
    \int_incr:N \stepsa
    \fp_sub:Nn \lambdaa { \schrittweite }
    \fp_set:Nn \differe
      {
        \lambdaa ** -0.5
        - ( ln(\lambdaa ** 0.5 * \RE)/ln(10) * 2 - 0.8 )
      }
    \typeout{differe~ \fp_use:N\differe}
    \typeout{stepsa~ \int_use:N\stepsa}
    \typeout{lambdaa~ \fp_use:N\lambdaa}
  }
\cs_set_eq:NN \fpeval \fp_eval:n
\ExplSyntaxOff

\begin{document}

\begin{equation}
\label{equ:z6}
    \frac{1}{\sqrt{\lambda}} = 2 * log (Re * \sqrt{\lambda}) -0,8
\end{equation}
\begin{equation}
\label{equ:z7}
    \frac{1}{\sqrt{\fpeval{\lambdaa}}} = 2 * log (\fpeval{\RE} *
\sqrt{\fpeval{\lambdaa}}) -0,8
\end{equation}

\end{document}

It would be even nicer to write a function which solves equations like this in general...

Answer 1

Solving the equation f(lambda) = g(lambda) that you have is equivalent to finding zeros of the function f(lambda) - g(lambda). The method you used was to compute this function at points which are all equally spaced. This is rather wasteful, because when you are far from the true zero, there is no need to have such a small step size.

Below you will find the implementation of three different methods to find zeros. First yours. Then the bisection method, which seems to work really well for your situation. I've chosen the starting interval to be [0, 1e9], but you can reduce it a lot, to reduce the number of steps that the algorithm does. The third method is the secant method, which converges faster, hence lets you get much more precise results if you want: its main drawback is that your starting values must be rather close to the correct zero, otherwise, you will not get anything.

\documentclass{article}
\usepackage{expl3, xparse}
\ExplSyntaxOn

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The function \fp_until_do:nn only exists since 2012-08-16.   %
% If it does not exist, emulate it with slightly slower, but   %
% entirely equivalent code.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\cs_if_exist:NF \fp_until_do:nn
  {
    \cs_new:Npn \fp_until_do:nn #1
      { \bool_until_do:nn { \fp_compare_p:n {#1} } }
  }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We will find it useful to define a function to get the sign  %
% of a floating point number.  There, we do not use the \fp_   %
% prefix: only kernel code should do this.  The function we    %
% define is called \my_fp_sign:n.  It gives 1 for positive     %
% numbers and -1 for everything else.  This could be improved, %
% but at the cost of adding an auxiliary function.             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\cs_new:Npn \my_fp_sign:n #1
  { \fp_eval:n { (#1) > 0 ? 1 : -1 } }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Declare variables                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\fp_new:N \xA
\fp_new:N \xMid
\fp_new:N \xB
\fp_new:N \xC
\fp_new:N \fnA
\fp_new:N \fnMid
\fp_new:N \fnB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method number 1 (your method): start from lambda = 0, and    %
% increase lambda step by step, until "myfn(lambda)" changes   %
% sign.  I changed 1e-5 to 1e-3 in your code, to make the test %
% reasonably fast, of course, this is very bad for precision,  %
% but in any case, you should use one of the other methods:    %
% this one is too slow.                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\cs_new_protected:Npn \methodI
  {
    \fp_zero:N \xA
    \fp_until_do:nn { (\fn{\xA}) < 0 }
      {
        \fp_add:Nn \xA { 1e-3 }
      }
    \msg_term:n { Result~(1):~\fp_use:N \xA }
  }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method number 2 (bisection): start from an interval where    %
% you expect the solution to be.  Split the inteval in two at  %
% each step, and use the sign of fn to decide which half of    %
% the interval to keep.                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\cs_new_protected:Npn \methodII
  {
    % The interval starts at [0, 1e5] (change if needed)
    %
    \fp_set:Nn \xA { 0 }
    \fp_set:Nn \xB { 1e5 }
    %
    % Compute fn(A) and fn(B)
    %
    \fp_set:Nn \fnA { \fn { \xA } }
    \fp_set:Nn \fnB { \fn { \xB } }
    %
    % Until the interval's size is < 1e-9, set "Mid" to be
    % "(A+B)/2", then compute fn(Mid).  If fn(Mid) and fn(B)
    % have opposite signs, we wish to keep the interval
    % [Mid, B], so set A = Mid, fnA = fnMid.  Otherwise,
    % we keep the interval [A, Mid], so we set B = Mid and
    % fnB = fnMid.
    %
    \fp_until_do:nn
      { \xB - \xA < 1e-9 }
      {
        \fp_set:Nn \xMid { ( \xA + \xB ) / 2 }
        \fp_set:Nn \fnMid { \fn{\xMid} }
        \fp_compare:nTF
          { \my_fp_sign:n { \fnMid } = \my_fp_sign:n { \fnB } }
          {
            \fp_set_eq:NN \xB \xMid
            \fp_set_eq:NN \fnB \fnMid
          }
          {
            \fp_set_eq:NN \xA \xMid
            \fp_set_eq:NN \fnA \fnMid
          }
      }
    \msg_term:n { Result~(2):~\fp_use:N \xA }
  }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method number 3 (secant):  this is a discrete version of     %
% Newton's method, since Newton's method requires us to know   %
% how to compute a derivative.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\cs_new_protected:Npn \methodIII
  {
    % This method is more sensitive to initial conditions.
    % Small initial values seem to work well with your function.
    %
    \fp_set:Nn \xA { 1e-5 }
    \fp_set:Nn \xB { 2e-5 }
    %
    % Compute fn(A) and fn(B)
    %
    \fp_set:Nn \fnA { \fn{\xA} }
    \fp_set:Nn \fnB { \fn{\xB} }
    %
    % Until |A - B| < 1e-9, compute
    %
    %   C = B - fnB * (fnA - fnB) / (A - B)
    %
    % then store (B, C) into (A, B).
    %
    \fp_until_do:nn
      { abs(\xB - \xA) < 1e-9 }
      {
        \fp_set:Nn \xC
          { \xB - (\fn{\xB}) * (\xA - \xB) / ((\fn{\xA}) - (\fn{\xB})) }
        \fp_set_eq:NN \xA \xB
        \fp_set_eq:NN \xB \xC
      }
    \fp_set_eq:NN \xA \xC
    \msg_term:n { Result~(3):~\fp_use:N \xA }
  }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluating and displaying floating point numbers.            %
% This could be improved using the not-yet-on-CTAN module      %
% l3str-format.                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\NewDocumentCommand { \fpeval } { om }
  {
    \IfValueTF {#1}
      { \fp_to_tl:n { round(#2,#1) } }
      { \fp_to_tl:n {#2} }
  }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display the result                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\NewDocumentCommand { \displayresult } { oo }
  {
    \begin{equation}
      \IfValueT {#1} { \label{equ:#1} }
      \frac{1}{\sqrt{\lambda}}
      = 2 * \log (\mathrm{Re} * \sqrt{\lambda}) - 0.8
    \end{equation}
    \begin{equation}
      \IfValueT {#2} { \label{equ:#2} }
      \frac{1}{\sqrt{\fpeval[9]{\xA}}}
      = 2 * \log (\fpeval[9]{\RE} * \sqrt{\fpeval[9]{\xA}}) - 0.8
    \end{equation}
  }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Value of \RE, and definition of our function \fn.            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\fp_new:N \RE
\fp_set:Nn \RE { 35800 }
\cs_new:Npn \fn #1
  {#1 ** -0.5 - (2 * ln(#1 ** 0.5 * \RE) / ln(10) - 0.8)}

% \cs_new:Npn \fn #1 {#1 ** 2 - 1}

\ExplSyntaxOff

\begin{document}

\methodI
\displayresult[z6][z7]

\methodII
\displayresult[z8][z9]

\methodIII
\displayresult[z10][z11]

\end{document}

Answer 2

This isn't a pure TeX solution; it uses my PythonTeX package. But this sort of approach can be used to solve most equations, without requiring you to write your own solvers and without ever requiring you to leave your TeX editor (assuming you set up a shortcut key for the PythonTeX script).

First, solving the equation using a very simple approach:

\documentclass{article} 
\usepackage{amsmath}
\usepackage{pythontex}

\begin{document}

\begin{pycode}
from numpy import sqrt, log10
from scipy.optimize import brentq
#Can't use "lambda", cause it's a keyword
def f(Lambda, Re):
    return 1/sqrt(Lambda) - 2 * log10(Re * sqrt(Lambda)) + 0.8
def get_lambda(Re):
    return brentq(f, 10**-12, 10**9, args=Re)
\end{pycode}

\begin{equation}
\frac{1}{\sqrt{\lambda}}=2\log_{10}\left(\text{Re}\times\sqrt{\lambda}\right)-0.8
\end{equation}
\begin{equation}
\frac{1}{\sqrt{\py{get_lambda(35800)}}}=2\log_{10}\left(35800\times\sqrt{\py{get_lambda(35800)}}\right)-0.8
\end{equation}

\end{document}

Second, a more advanced approach, that automatically generates the TeX code for arbitrary Re:

\documentclass{article} 
\usepackage{amsmath}
\usepackage{pythontex}

\begin{document}

\begin{pycode}
from numpy import sqrt, log10
from scipy.optimize import brentq
#Can't use "lambda", cause it's a keyword
def f(Lambda, Re):
    return 1/sqrt(Lambda) - 2 * log10(Re * sqrt(Lambda)) + 0.8
#Create a TeX form of the equation
f_tex = r'\frac{1}{\sqrt{\lambda}}=2\log_{10}\left(\text{Re}\times\sqrt{\lambda}\right)-0.8'

#Create a function that will typeset the equation and a solved form
def typeset(Re):
    ans = ''
    ans += r'\begin{equation}' + '\n'
    ans += f_tex + '\n'
    ans += r'\end{equation}' + '\n'
    ans += r'\begin{equation}' + '\n'
    ans += f_tex.replace(r'\text{Re}', str(Re)).replace(r'\lambda', str(brentq(f, 10**-12, 10**9, args=Re))) + '\n'
    ans += r'\end{equation}' + '\n'
    return ans

#Could just do "typeset(35800)" here
#But will use "\py" command instead
#It could be used anywhere, with any value for Re
\end{pycode}

\py{typeset(35800)}

\end{document}