# Iterative method to solve nonlinear equations

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}

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

\end{document}


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

-
TeX is not a mathematical engine, it's a typesetting system. I'm not quite sure what question you are seeking to answer here, but wouldn't it be most sensible to do what seems to be a complex calculation in a suitable specialist tool and include the result in your typeset material? –  Joseph Wright Aug 17 '12 at 10:41
Thats for sure a possibility, but would't that be a reason not to progress on the work on l3fp? The goal in this case is to dynamically adopt the stepwide.... –  Peter Ebelsberger Aug 17 '12 at 10:46
The reason we wanted an FPU in LaTeX3 was primarily to support rotation, which does need floating point support (sine and cosine). There are often common cases for slightly more involved calculations 'within' TeX, for example as seen in pgfplots, TikZ or some of the spreadsheet-like packages. However, all of those cases involve balancing convenience (doing in TeX) with performance/accuracy, and some of these tasks using less accurate but faster dimen-based code. The task you seem to want to carry out is very complex, and I'm not sure TeX is the best tool. –  Joseph Wright Aug 17 '12 at 10:52
Quite apart from that, I'm not really sure what TeX question we are being asked to answer. –  Joseph Wright Aug 17 '12 at 10:55
Asking something like 'What is the most efficient way to do this calculation using a TeX-based system?' could be on topic, and could attract answers using pgfmath, l3fp, fp, Lua code, etc., with discussion of the underlying implementations. On the other hand, simply asking about the best algorithm is more a mathematical question. –  Joseph Wright Aug 17 '12 at 11:06

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 }
{
}
\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 }
{
$$\IfValueT {#1} { \label{equ:#1} } \frac{1}{\sqrt{\lambda}} = 2 * \log (\mathrm{Re} * \sqrt{\lambda}) - 0.8$$
$$\IfValueT {#2} { \label{equ:#2} } \frac{1}{\sqrt{\fpeval[9]{\xA}}} = 2 * \log (\fpeval[9]{\RE} * \sqrt{\fpeval[9]{\xA}}) - 0.8$$
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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}

-
Very impressive, Bruno! :) –  Paulo Cereda Aug 19 '12 at 0:37
@PauloCereda Thanks. I'm actually disappointed that I couldn't have a general algorithm that would work for arbitrary functions. One issue is to find the domain on which a function is defined. If that's known then bisection should work reasonably well. –  Bruno Le Floch Aug 19 '12 at 0:41
@BrunoLeFloch As I've commented on the question, that is more a maths question than a LaTeX one. –  Joseph Wright Aug 19 '12 at 9:21
Great, thank you very lot. –  Peter Ebelsberger Aug 20 '12 at 10:38

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}

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

\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'$$' + '\n' ans += f_tex + '\n' ans += r'$$' + '\n'
ans += r'$$' + '\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'$$' + '\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}

-
@GPoore Is there a reason why PythonTeX is not on CTAN? Having all packages at the same place is helpful to users. –  Bruno Le Floch Aug 20 '12 at 16:19
@BrunoLeFloch I will be putting it on CTAN soon, but I've been waiting till I feel comfortable declaring it out of beta. Parts of Python don't quite behave the same under Linux and Windows, and there are also the changes between Python 2 and 3, so it's been a little trickier than the typical package. –  G. Poore Aug 20 '12 at 20:17
I will try to install it manually and so I am ready I will give you a feedback, if ok? –  Peter Ebelsberger Aug 20 '12 at 21:30
@G.Poore That makes sense. You may also want to test the behaviour for non-Ascii strings, which are treated differently in Python 2 and 3, and are treated differently in pdfTeX (manipulates bytes, not characters) versus XeTeX and LuaTeX (Unicode-aware). –  Bruno Le Floch Aug 21 '12 at 12:44