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