3

In a macro that is supposed to calculate the value of the third derivative, I used the following formula:

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

Because of the 3rd power of delta in the denominator (delta should be as small as possible), the expression becomes very unstable and only in a very small range for delta can one obtain approximately a reasonable solution.

Is there a way to improve that? I do not need very high accuracy, but it would be good if the result did not depend so much on the chosen delta.

Many thanks for your help.

enter image description here

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

With the advice of Christian Hupfer on central differences and Andrew Swann on finite difference coefficient, I have found a reasonably stable (certainly not in all cases) solution for the first three derivatives. In the examples, the calculated values are not very different from those calculated using Maple.

enter image description here

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

closed as off-topic by egreg, CarLaTeX, Stefan Pinnow, TeXnician, Kurt Dec 23 '17 at 11:09

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question does not fall within the scope of TeX, LaTeX or related typesetting systems as defined in the help center." – egreg, CarLaTeX, Stefan Pinnow, TeXnician, Kurt
If this question can be reworded to fit the rules in the help center, please edit the question.

  • Do you have to determine the derivative numerically? If not, you can maybe use the analytical derivative instead. – Dr. Manuel Kuehner Dec 23 '17 at 6:04
  • @Dr.ManuelKuehner that would be the best, using the analytical derivative. Is there a way to calculate this automatically? – user139826 Dec 23 '17 at 6:06
  • 4
    This is a task for a proper software like Maple, Mathematica, Matlab (MuPad) or similar freeware solutions. In the given examples that you use, you can do it by hand. – Dr. Manuel Kuehner Dec 23 '17 at 6:34
  • 1
    @Thomas: Your Differences - Scheme is not suited for a 5th order polynomial. You need a higher order which 'proves' to be exact for polynomials to the same order – user31729 Dec 23 '17 at 6:41
  • @ChristianHupfer, I dont think so. The same problem is with all other functions, for example: f(x)=x^2 \cdot \sin(x)or f(x)=x^3 – user139826 Dec 23 '17 at 8:24
7

If you are working only with polynomials then you can implement the derivative of polynomials analytically. This can be done using TeX language. For example:

\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

This prints the result:

derivative

8

The LaTeX3 FPU is designed to carry out 'convenient' arithmetic to support a range of 'reasonable' tasks within a LaTeX run. In particular, it sets out to implement the requirements of IEEE754. 'Convenient' arithmetic covers things like working out rotations of images, adding up columns of numbers in tables, etc.: certainly not mathematical analysis. In that sense, a fair comparison is what a typical spreadsheet can be used for. If you try the demonstration in the question in a well-known spreadsheet application, you'll find that the same variations in stability are seen: this is simply out-of-scope for the FPU. As such, this is a task for a specialist application.

  • -@Joseph Wright thank you for your explanation. – user139826 Dec 23 '17 at 9:09