我正在尝试使用 TeX 求解以下方程（来自流体动力学）。
1/sqrt(lambda) = 2 * log(Re * sqrt(lambda)) - 0.8
其中lambda为未知数，sqrt表示“平方根”，log为十进制对数，Re为给定数。
这没法直接解决...不能写成lambda = ...。
我如何改进代码或下面的方法以更快地找到零？
我已经在包中编写了代码fp，并在 Bruno Le Floch 的帮助下编写了一个使用 LaTeX3 等效版本进行编程l3fp。
方法就是从小lambda处开始，然后逐渐增加，直到左侧和右侧之间的差异变为（接近）零。问题是有一个恒定的步骤，所以我需要做很多步骤才能找到lambda。但我希望步骤数量最少。
这里是 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}
编写一个可以解决此类方程的函数就更好了……
答案1
解出这个方程f(lambda) = g(lambda)相当于找到函数的零点f(lambda) - g(lambda)。你使用的方法是计算该函数在等距点处的值。这相当浪费，因为当你远离真实零点时，就没有必要使用这么小的步长。
下面你将看到三种不同的方法来找到零点。首先是你的。然后二分法，这似乎非常适合您的情况。我选择的起始间隔是[0, 1e9]，但您可以将其大大减少，以减少算法执行的步骤数。第三种方法是割线法，其收敛速度更快，因此如果您愿意的话，可以得到更精确的结果：其主要缺点是您的起始值必须非常接近正确的零，否则，您将什么也得不到。
\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}
答案2
这不是一个纯粹的 TeX 解决方案；它使用了我的PythonTeX包。但这种方法可用于解决大多数方程式，而不需要您编写自己的求解器，也不需要您离开 TeX 编辑器（假设您为 PythonTeX 脚本设置了快捷键）。
首先，用一个非常简单的方法解这个方程：
\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}
第二种方法更先进，可以自动生成任意 Re 的 TeX 代码：
\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}