☰

数值稳定性问题

tag-icon • tag-icon macros calculations latex3 xparse l3fp
数值稳定性问题

在应该计算三阶导数值的宏中，我使用了以下公式：

(f (x + 3 * delta) -3 * f (x + 2 * delta) + 3 * f (x + delta) -f (x)) / (delta) ^ 3

由于分母中含有delta的3次方（delta应该尽量小），导致表达式变得很不稳定，只有在delta的一个很小的范围内才能近似的得到合理的解。

有没有什么办法可以改善这种情况？我不需要非常高的准确度，但如果结果不太依赖于所选的增量，那就更好了。

非常感谢您的帮助。

在此处输入图片描述

\documentclass{article}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[ngerman]{babel}
\usepackage{siunitx}

\begin{document}

\ExplSyntaxOn

\NewDocumentCommand{\FuncValue}{ t' t' t' O{} m m }
{% #1 = option list, #2 = value, #3 = function
\group_begin:
 \keys_set:nn { th/FV } { #4 }
  \tl_set:Nn \l_th_funk_in_tl { #5 }
    \regex_replace_all:nnN { pi } { \c{pi} } \l_th_funk_in_tl
    \regex_replace_all:nnN { ([+|\-]*)(\c[^BE].*)(/)(\c[^BE].*) } { \1\c{frac} \cB\{ \2\cE\} \cB\{ \4\cE\} } \l_th_funk_in_tl
   \IfBooleanTF{#3}% derivative of 3rd, 2nd, first order and function
    { 
     f^{\prime\prime\prime}\left(\l_th_funk_in_tl\right) =
     \th_funcDDD_value:nn { #5 } { #6 }
    } 
    {  }
\group_end:
}

\cs_generate_variant:Nn \cs_set:Nn { NV }

\keys_define:nn { th/FV }
{
 round .int_set:N  = \l__th_FV_round_int,
 round .initial:n  = 3,
 delta .tl_set:N   = \l__th_FV_delta_tl,
 delta .initial:n  = 1e-4,
}

\cs_new_protected:Nn \th_funcDDD_value:nn
{
  \tl_set:Nn \l_th_funkDDD_tl { #2 }
    \regex_replace_all:nnN { x } { \cP\#1 } \l_th_funkDDD_tl
  \cs_set:NV \__th_FVDDD_function:n  \l_th_funkDDD_tl
   \fp_eval:n
    {
      round( 
             ( 
                    \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 3*\l__th_FV_delta_tl ) } } }   %    f(x+3*delta)
             - ( 3* \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 2*\l__th_FV_delta_tl ) } } } ) % -3*f(x+2*delta)
             + ( 3* \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { \fp_eval:n { #1 } + \l__th_FV_delta_tl } } } )       % +3*f(x+delta)
                  - \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { #1 } } }                                             % -  f(x)
              )                                         % ( f(x+3*delta)-3*f(x+2*delta)+3*f(x+delta)-f(x) ) / ( delta )^3
              / 
              ( \fp_eval:n {\l__th_FV_delta_tl } )**3 
             , \l__th_FV_round_int )
  }
}

\ExplSyntaxOff

$ f(x) = 4\cdot x^2-2\cdot  x^3+x^5 ; \qquad \FuncValue'''{pi/2}{4*x^2-2*x^3+x^5} $

$ f(x) = 4\cdot x^2-2\cdot  x^3+x^5 ; \qquad \FuncValue'''[delta=1e-5]{pi/2}{4*x^2-2*x^3+x^5} $

$ f(x) = 4\cdot x^2-2\cdot  x^3+x^5 ; \qquad \FuncValue'''[delta=1e-6]{pi/2}{4*x^2-2*x^3+x^5} $

\end{document}

在 Christian Hupfer 关于中心差分和 Andrew Swann 关于有限差分系数的建议下，我找到了前三个导数的合理稳定（当然不是在所有情况下）解。在示例中，计算值与使用 Maple 计算的值没有太大差别。

在此处输入图片描述

\documentclass{article}
\usepackage[margin=1.3cm]{geometry}
\usepackage{siunitx}
\usepackage{xintexpr}

\def\dprime{\prime\prime}
\def\trprime{\prime\prime\prime}

\ExplSyntaxOn

\NewDocumentCommand{\FuncValue}{ t' t' t' O{} m m }
{% #1 für 1.Abl, #2 für 2.Abl, #3 für 3.Abl  #4 = option list, #5 = value, #6 = function
\group_begin:
\keys_set:nn { thomas/FuncValue } { #4 }
\tl_set:Nn \l_thomas_funk_in_tl { #5 }
\regex_replace_all:nnN { \. } { , }        \l_thomas_funk_in_tl
\regex_replace_all:nnN { pi } { \c{pi} }   \l_thomas_funk_in_tl
\regex_replace_all:nnN { ([+|\-]*)(\c[^BE].*)(/)(\c[^BE].*) } { \1\c{frac} \cB\{ \2\cE\} \cB\{ \4\cE\} }   \l_thomas_funk_in_tl
\tl_set:Nn \l_thomas_funk_out_tl { #6 }
\regex_replace_all:nnN { \. } { , }             \l_thomas_funk_out_tl
\regex_replace_all:nnN { \* } { \c{cdot} }      \l_thomas_funk_out_tl
\regex_replace_all:nnN { (\c{cdot})(x) } { x }  \l_thomas_funk_out_tl
\regex_replace_all:nnN { (\c{cdot})(sqrt) } { sqrt } \l_thomas_funk_out_tl
\regex_replace_all:nnN { pi } { \c{pi} }        \l_thomas_funk_out_tl
\regex_replace_all:nnN { sin|sind } { \c{sin} } \l_thomas_funk_out_tl
\regex_replace_all:nnN { cos|cosd } { \c{cos} } \l_thomas_funk_out_tl
\regex_replace_all:nnN { tan|tand } { \c{tan} } \l_thomas_funk_out_tl
\regex_replace_all:nnN { ln } { \c{ln} }        \l_thomas_funk_out_tl
\regex_replace_all:nnN { (sqrt)(\{\c[^BE].*\}) } { \c{sqrt} \2 }               \l_thomas_funk_out_tl
\regex_replace_all:nnN { (sqrt\( )(\c[^BE].*)(\)) } { \c{sqrt} \cB\{ \2\cE\} } \l_thomas_funk_out_tl
\regex_replace_all:nnN { (\^)(\()(\c[^BE].*)(\)) } { \1\cB\{\3\cE\} }          \l_thomas_funk_out_tl
\regex_replace_all:nnN { (\() (\c[^BE].*) (\)/) (\c[^BE].*) } { \c{frac} \cB\{ \2\cE\} \cB\{ \4\cE\} }    \l_thomas_funk_out_tl
\regex_replace_all:nnN { (\c[^BE].*) (/\() (\c[^BE].*) (\)) } { \c{frac} \cB\{ \1\cE\} \cB\{ \3\cE\} }    \l_thomas_funk_out_tl
\regex_replace_all:nnN { (\c{frac}\{\c[^BE].*\}) (\{\( ) (\c[^BE].*) (\)\}) } { \1 \cB\{ \3\cE\} }        \l_thomas_funk_out_tl
\regex_replace_all:nnN { ( (\d+|\d*\.\d+)\~* ) (/) ( (\d+|\d*\.\d+)\~*)  } { \c{frac} \cB\{ \2\cE\} \cB\{ \4\cE\} } \l_thomas_funk_out_tl
%\regex_replace_all:nnN { (\c[^BE].*) (/) (\c[^BE].*) } { \c{frac} \cB\{ \1\cE\} \cB\{ \3\cE\} }           \l_thomas_funk_out_tl
\regex_replace_all:nnN { [+-]?\d+e[+-]?\d+ } { \c{num} \cB\{ \0\cE\} }                                    \l_thomas_funk_out_tl

 \ensuremath % <-- is it really needed? Egreg! No, it is not really needed. Is there a problem with ensuremath in this case?
  {
   \bool_if:NT \l__thomas_FuncValue_fkt_bool
     {
       \l__th_FV_name_tl(x) = \l_thomas_funk_out_tl; \qquad
     }
   \IfBooleanTF{#3}%
   { \bool_if:NT \l__thomas_FuncValue_abl_bool
     { \l__th_FV_name_tl^{\trprime}\negthinspace\left(\textstyle\l_thomas_funk_in_tl\right) = }
     \thomas_funcDDD_value:nn { #5 } { #6 } }
     {
   \IfBooleanTF{#2}%
   {  \bool_if:NT \l__thomas_FuncValue_abl_bool
     { \l__th_FV_name_tl^{\dprime}\negthinspace\left(\textstyle\l_thomas_funk_in_tl\right) = }
     \thomas_funcDD_value:nn { #5 } { #6 } }
   {
   \IfBooleanTF{#1}%
   {  \bool_if:NT \l__thomas_FuncValue_abl_bool
     { \l__th_FV_name_tl^{\prime}\negthinspace\left(\textstyle\l_thomas_funk_in_tl\right) = }
     \thomas_funcD_value:nn { #5 } { #6 } }
   {  \bool_if:NT \l__thomas_FuncValue_abl_bool
     { \l__th_FV_name_tl\negthinspace\left(\textstyle\l_thomas_funk_in_tl\right) = }
     \thomas_function_value:nn { #5 } { #6 } }
   }
   }
   }
\group_end:
}

\cs_new_protected:Nn \thomas_function_value:nn
{
 \tl_set:Nn \l_thomas_funk_tl { #2 }
 \regex_replace_all:nnN { x } { (x) } \l_thomas_funk_tl
 \regex_replace_all:nnN { x } { \cP\#1 } \l_thomas_funk_tl
 \cs_set:NV \__thomas_functionValue_function:n  \l_thomas_funk_tl
 \bool_if:NTF \l__thomas_FuncValue_frac_bool
 { \xintSignedFrac { \xintIrr {
     \fp_eval:n
  {
   round( \__thomas_functionValue_function:n { \fp_eval:n { #1 } } , \l__th_FV_round_int )
  } } }
  }
 { \num{
 \fp_eval:n
  {
   round( \__thomas_functionValue_function:n { \fp_eval:n { #1 } } , \l__th_FV_round_int )
  }
  } }
}

\cs_generate_variant:Nn \cs_set:Nn { NV }

\keys_define:nn { thomas/FuncValue }
{
 fkt   .bool_set:N = \l__thomas_FuncValue_fkt_bool,
 fkt   .initial:n  = false,
 fkt   .default:n  = true,
 frac  .bool_set:N = \l__thomas_FuncValue_frac_bool,
 frac  .initial:n  = false,
 fkt   .default:n  = true,
 abl   .bool_set:N = \l__thomas_FuncValue_abl_bool,
 abl   .initial:n  = true,
 abl   .default:n  = true,
 round .int_set:N  = \l__th_FV_round_int,
 round .initial:n  = 3,
 eps   .tl_set:N   = \l__th_FV_epsilon_tl,
 eps   .initial:n  = 1e-3,
 name  .tl_set:N   = \l__th_FV_name_tl,
 name  .initial:n  = f,
}

\cs_new_protected:Nn \thomas_funcD_value:nn
{
 \tl_set:Nn \l_th_funkD_tl { #2 }
 \regex_replace_all:nnN { x } { (x) } \l_th_funkD_tl
 \regex_replace_all:nnN { x } { \cP\#1 } \l_th_funkD_tl
 \cs_set:NV \__th_FVD_function:n  \l_th_funkD_tl
 \tl_set:Nn \l_th_funkD_Wert_tl {
   \fp_eval:n
    {
     round(
            (
               ( \fp_eval:n { 1/280  } * \fp_eval:n { \__th_FVD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 4*\l__th_FV_epsilon_tl ) } } } )  %
             - ( \fp_eval:n { 4/105 }  * \fp_eval:n { \__th_FVD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 3*\l__th_FV_epsilon_tl ) } } } )  %
             + ( \fp_eval:n { 1/5 }    * \fp_eval:n { \__th_FVD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 2*\l__th_FV_epsilon_tl ) } } } )  %
             - ( \fp_eval:n { 4/5 }    * \fp_eval:n { \__th_FVD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 1*\l__th_FV_epsilon_tl ) } } } )  %
             + ( \fp_eval:n { 4/5 }    * \fp_eval:n { \__th_FVD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 1*\l__th_FV_epsilon_tl ) } } } )  %
             - ( \fp_eval:n { 1/5 }    * \fp_eval:n { \__th_FVD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 2*\l__th_FV_epsilon_tl ) } } } )  %
             + ( \fp_eval:n { 4/105 }  * \fp_eval:n { \__th_FVD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 3*\l__th_FV_epsilon_tl ) } } } )  %
             - ( \fp_eval:n { 1/280 }  * \fp_eval:n { \__th_FVD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 4*\l__th_FV_epsilon_tl ) } } } )  %
            )                                         % ( f(x+3*epsilon)-3*f(x+2*epsilon)+3*f(x+epsilon)-f(x) ) / ( epsilon )^3
             /
               ( \fp_eval:n {\l__th_FV_epsilon_tl } )
             , \l__th_FV_round_int )
  }
  }
  \bool_if:NTF \l__thomas_FuncValue_frac_bool
   { \xintSignedFrac { \xintIrr { \l_th_funkD_Wert_tl } } }
   { \num{ \l_th_funkD_Wert_tl } }
}

\cs_new_protected:Nn \thomas_funcDD_value:nn
{
 \tl_set:Nn \l_th_funkDD_tl { #2 }
 \regex_replace_all:nnN { x } { (x) } \l_th_funkDD_tl
 \regex_replace_all:nnN { x } { \cP\#1 } \l_th_funkDD_tl
 \cs_set:NV \__th_FVDD_function:n  \l_th_funkDD_tl
  \tl_set:Nn \l_th_funkDD_Wert_tl {
   \fp_eval:n
    {
     round(
            (
               ( \fp_eval:n { -1/560 } * \fp_eval:n { \__th_FVDD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 4*\l__th_FV_epsilon_tl ) } } } )   %
             + ( \fp_eval:n { 8/315 }  * \fp_eval:n { \__th_FVDD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 3*\l__th_FV_epsilon_tl ) } } } )   %
             - ( \fp_eval:n { 1/5 }    * \fp_eval:n { \__th_FVDD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 2*\l__th_FV_epsilon_tl ) } } } )   %
             + ( \fp_eval:n { 8/5 }    * \fp_eval:n { \__th_FVDD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 1*\l__th_FV_epsilon_tl ) } } } )   %
             - ( \fp_eval:n { 205/72 } * \fp_eval:n { \__th_FVDD_function:n { \fp_eval:n { #1 } } } )                                               %
             + ( \fp_eval:n { 8/5 }    * \fp_eval:n { \__th_FVDD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 1*\l__th_FV_epsilon_tl ) } } } )   %
             - ( \fp_eval:n { 1/5 }    * \fp_eval:n { \__th_FVDD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 2*\l__th_FV_epsilon_tl ) } } } )   %
             + ( \fp_eval:n { 8/315 }  * \fp_eval:n { \__th_FVDD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 3*\l__th_FV_epsilon_tl ) } } } )   %
             - ( \fp_eval:n { 1/560 }  * \fp_eval:n { \__th_FVDD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 4*\l__th_FV_epsilon_tl ) } } } )   %
            )                                         % ( f(x+3*epsilon)-3*f(x+2*epsilon)+3*f(x+epsilon)-f(x) ) / ( epsilon )^3
             /
               ( \fp_eval:n {\l__th_FV_epsilon_tl } )**2
             , \l__th_FV_round_int )
  }
  }
   \bool_if:NTF \l__thomas_FuncValue_frac_bool
   { \xintSignedFrac { \xintIrr { \l_th_funkDD_Wert_tl } } }
   { \num{ \l_th_funkDD_Wert_tl } }
}

\cs_new_protected:Nn \thomas_funcDDD_value:nn
{
 \tl_set:Nn \l_th_funkDDD_tl { #2 }
 \regex_replace_all:nnN { x } { (x) } \l_th_funkDDD_tl
 \regex_replace_all:nnN { x } { \cP\#1 } \l_th_funkDDD_tl
 \cs_set:NV \__th_FVDDD_function:n  \l_th_funkDDD_tl
 \tl_set:Nn \l_th_funkDDD_Wert_tl {
   \fp_eval:n
    {
     round(
            (
               ( \fp_eval:n { -7/240  } * \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 4*\l__th_FV_epsilon_tl ) } } } )  %
             + ( \fp_eval:n { 3/10 }    * \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 3*\l__th_FV_epsilon_tl ) } } } )  %
             - ( \fp_eval:n { 169/120 } * \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 2*\l__th_FV_epsilon_tl ) } } } )  %
             + ( \fp_eval:n { 61/30 }   * \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { \fp_eval:n { #1 } - ( 1*\l__th_FV_epsilon_tl ) } } } )  %
             - ( \fp_eval:n { 61/30 }   * \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 1*\l__th_FV_epsilon_tl ) } } } )  %
             + ( \fp_eval:n { 169/120 } * \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 2*\l__th_FV_epsilon_tl ) } } } )  %
             - ( \fp_eval:n { 3/10 }    * \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 3*\l__th_FV_epsilon_tl ) } } } )  %
             + ( \fp_eval:n { 7/240 }   * \fp_eval:n { \__th_FVDDD_function:n { \fp_eval:n { \fp_eval:n { #1 } + ( 4*\l__th_FV_epsilon_tl ) } } } )  %
            )                                         %
             /
               ( \fp_eval:n {\l__th_FV_epsilon_tl } )**3
             , \l__th_FV_round_int )
  }
  }
    \bool_if:NTF \l__thomas_FuncValue_frac_bool
   { \xintSignedFrac { \xintIrr { \l_th_funkDDD_Wert_tl } } }
   { \num{ \l_th_funkDDD_Wert_tl } }
}

\NewDocumentCommand{\FuncValueSet}{ m }{ \keys_set:nn { thomas/FuncValue } { #1 } }

\ExplSyntaxOff

\begin{document}

\FuncValueSet{round=10}

\makebox[8.5cm][l]{\FuncValue[fkt]{21/3}{(4*x^2-2)/(3*x^3+x^5)}}
\makebox[4.5cm][l]{\FuncValue'{21/3}{(4*x^2-2)/(3*x^3+x^5)} }
                   \FuncValue''{21/3}{(4*x^2-2)/(3*x^3+x^5)}

\makebox[8.5cm][l]{\FuncValue[fkt]{-7}{(4*x^2-2)/(3*x^3+x^5)} }
\makebox[4.5cm][l]{\FuncValue'{-7}{(4*x^2-2)/(3*x^3+x^5)} }
                   \FuncValue''{-7}{(4*x^2-2)/(3*x^3+x^5)}

\makebox[8.5cm][l]{\FuncValue[fkt]{-2/3}{1/(2*x^2-2*x^3+x^5)} }
\makebox[4.5cm][l]{\FuncValue'{-2/3}{1/(2*x^2-2*x^3+x^5)} }
                   \FuncValue''{-2/3}{1/(2*x^2-2*x^3+x^5)}

\makebox[8.5cm][l]{\FuncValue[fkt]{2/3}{1/(2*x^2-2*x^3+x^5)} }
\makebox[4.5cm][l]{\FuncValue'{2/3}{1/(2*x^2-2*x^3+x^5)} }
                   \FuncValue''{2/3}{1/(2*x^2-2*x^3+x^5)}

\FuncValue'''[fkt,eps=2e-2]{pi/2}{-2*x^3+x^5}\quad mit Maple18 berechnet: \num{136.04406601634037929}

\end{document}

答案1

如果您只使用多项式，那么您可以分析地实现多项式的导数。这可以使用 TeX 语言来完成。例如：

\input opmac
\input apnum

\def\derivative#1to#2{\expandafter\derA#1\relax#2}

\def\derA#1\relax{\def\tmpb{#1&=}%
   \replacestrings+{&+}\replacestrings-{&-}\replacestrings{ }{}%
   \edef\tmpb{\expandafter}\expandafter\derB\tmpb
}
\def\derB#1&#2{%
   \ifx=#2% the end
      \ifx&#1&\else\derC#1x!&\fi
      \ifx\tmpb\empty\def\tmpb{0}\fi
      \expandafter\derF
   \else
      \ifx&#1&\else\derC#1x!&\fi
      \expandafter\derB\expandafter#2%
   \fi
}
\def\derC#1x#2&{\ifx!#2\else % derivative of constant is zero
   \derD#1x#2x^{}&\fi
}
\def\derD#1x^#2#3&{\ifx&#2&\derE#1&% x^1
   \else
      \evaldef\c{#1#2}\ifnum\apSIGN>0 \edef\c{+\c}\fi
      \evaldef\e{#2-1}%
      \edef\tmpb{\tmpb\c \ifnum\e=0 \else *x\ifnum\e>9 ^{\e}\else\ifnum\e>1 ^\e\fi\fi\fi}
   \fi
}
\def\derE#1x#2&{\evaldef\c{#11}\ifnum\apSIGN>0 \edef\c{+\c}\fi \edef\tmpb{\tmpb\c}}
\def\derF#1{\let#1=\tmpb \message{>>>> \tmpb}}
\def\evalfunc#1=#2(#3){\let\tmpb=#2\replacestrings{x}{(#3)}%
   \expandafter\evaldef\expandafter#1\expandafter{\tmpb}\apROUND#1{10}}

\mathcode`*="2201

\def\p{4*x^2-2*x^3+x^5+7}

\derivative\p to\q
\derivative\q to\q
\derivative\q to\q
\evalfunc\v=\q(\PIhalf)

$ f(x) = \p ; \qquad f'''(x) = \q; \quad f'''(\pi/2) = \v $

\bye

打印结果：

衍生物

答案2

LaTeX3 FPU 旨在执行“便捷”的算法，以支持 LaTeX 运行中的一系列“合理”任务。特别是，它着手实现 IEEE754 的要求。“便捷”的算法包括计算图像的旋转、将表格中的数字列相加，ETC。：当然不是数学分析。从这个意义上说，公平的比较是典型的电子表格的用途。如果您在著名的电子表格应用程序中尝试问题中的演示，您会发现稳定性的变化相同：这完全超出了 FPU 的范围。因此，这是专业应用程序的任务。

相关内容