Version 2 of Evaluating polynomial functions

Updated 2004-07-13 11:37:47 by TV

if { 0 } { Arjen Markus (13 july 2004) I have created a small package for dealing with polynomial functions: you can add, subtract, and multiply polynomial functions. You can divide a polynomial by another polynomial (with a remainder) and you can, of course, evaluate them at a particular coordinate.

The latter is not as simple as it may look at first. A naive computation like this:

`   f(x) = 1.0 + 2.0 * x + 3.0 * x**2 + 4.0 * x**3 + ... + 10.0 * x**9`

may lead to very inaccurate results plus it uses far too many multiplications. Horner's rule can be used to improve the results:

`   f(x) = 1.0 + x * (2.0 + x * (3.0 + x * (4.0 + ... x * 10.0)))))))))`

This is not however a guarantee to success: the various families of (classical) orthogonal polynomials, like Chebyshev and Legendre polynomials, cause significant trouble even with this rule, as kbk pointed out to me.

To appreciate how much trouble, let us do some experiments: }

``` #
# Load the scripts for general and orthogonal polynomials
#
source polynomials.tcl
source classic_polyns.tcl

#
# A simple polynomial: g(x) = (x-1.3)**20
# A second polynomial: h(x) = x*(x-1.1)*(x-2.1)*(x-3.1)*...*(x-19.1)
#
set f [::math::polynomials::polynomial {-1.3 1}]
set g 1
set h 1
for {set i 0 } {\$i < 20} {incr i} {
set g [::math::polynomials::multPolyn \$f \$g]
set h [::math::polynomials::multPolyn \$h \
[::math::polynomials::polynomial [list [expr {-\$i-0.1}] 1]]]
}

#
# Some checks:
#
puts "Degree:       [::math::polynomials::degreePolyn \$g]"
puts "Coefficients: [::math::polynomials::allCoeffsPolyn \$g]"

#
# Evaluate g at x=1.3 in two ways - naively and via Horner
#
puts "g:"
puts "x, \"exact\", naive, Horner -- relative errors"
foreach x {0.1 0.2 0.3 0.5 1.0 1.2 1.3 1.4 1.6 2.0 5.0} {
set result1 0.0
set power   0
foreach c [::math::polynomials::allCoeffsPolyn \$g] {
set result1 [expr {\$result1+\$c*pow(\$x,\$power)}]
incr power
}

set result2 [::math::polynomials::evalPolyn \$g \$x]
set exact   [expr {pow(\$x-1.3,20)}]
if { \$exact != 0.0 } {
set error1  [expr {(\$result1-\$exact)/\$exact}]
set error2  [expr {(\$result2-\$exact)/\$exact}]
puts [format "%4.2f %20.12e %20.12e %20.12e -- %10e %10e" \
\$x \$exact \$result1 \$result2 \$error1 \$error2]
} else {
set error1  -
set error2  -
puts [format "%4.2f %20.12e %20.12e %20.12e -- %10s %10s" \
\$x \$exact \$result1 \$result2 \$error1 \$error2]
}
}
puts "h:"
puts "x, \"exact\", naive, Horner -- relative errors"
foreach x {0.1 0.2 0.3 0.5 1.01 1.2 1.3 1.4 1.6 2.1 5.1} {
set result1 0.0
set power   0
foreach c [::math::polynomials::allCoeffsPolyn \$h] {
set result1 [expr {\$result1+\$c*pow(\$x,\$power)}]
incr power
}

set result2 [::math::polynomials::evalPolyn \$h \$x]
set exact   1.0
for { set i 0 } { \$i < 20 } { incr i } {
set exact   [expr {\$exact*(\$x-\$i-0.1)}]
}
if { \$exact != 0.0 } {
set error1  [expr {(\$result1-\$exact)/\$exact}]
set error2  [expr {(\$result2-\$exact)/\$exact}]
puts [format "%4.2f %20.12e %20.12e %20.12e -- %10e %10e" \
\$x \$exact \$result1 \$result2 \$error1 \$error2]
} else {
set error1  -
set error2  -
puts [format "%4.2f %20.12e %20.12e %20.12e -- %10s %10s" \
\$x \$exact \$result1 \$result2 \$error1 \$error2]
}
}

#
# Check Legendre 2, 4, ..., 50
#
puts "order (n), Pn(1), Pn(-1)"
puts "expected: 1 and 1"
for { set i 2 } { \$i <= 50 } { incr i 2 } {
set leg [::math::special::legendre \$i]
set p1 [::math::polynomials::evalPolyn \$leg 1.0]
set p2 [::math::polynomials::evalPolyn \$leg -1.0]

puts [format "%5d %20.12g %20.12g" \$i \$p1 \$p2]
}```

The result is:

``` Degree:       20
Coefficients: 190.049637749 -2923.84058075 21366.5273209 -98614.7414812 322394.347149 -793586.085291 1526127.0871 -2347887.82631 2934859.78288 -3010112.59783 2547018.352 -1781131.71469 1027575.98924 -486426.503784 187087.11684 -57565.26672 13837.8045 -2504.58 321.1 -26.0 1.0
g:
x, "exact", naive, Horner -- relative errors
0.10   3.833759992447e+01   3.833759992476e+01   3.833759992476e+01 -- 7.338657e-12 7.316973e-12
0.20   6.727499949326e+00   6.727499948326e+00   6.727499947888e+00 -- -1.486001e-10 -2.137647e-10
0.30   1.000000000000e+00   9.999999934125e-01   9.999999919215e-01 -- -6.587529e-09 -8.078473e-09
0.50   1.152921504607e-02   1.152913820723e-02   1.152913493922e-02 -- -6.664706e-06 -6.948161e-06
1.00   3.486784401000e-11  -1.930119685767e-05  -1.431931048046e-05 -- -5.535539e+05 -4.106748e+05
1.20   1.000000000000e-20  -1.147058508124e-04  -7.339982499843e-05 -- -1.147059e+16 -7.339982e+15
1.30   0.000000000000e+00  -2.580087965214e-04  -1.517559886679e-04 --          -          -
1.40   1.000000000000e-20  -5.542161987933e-04  -2.977163204889e-04 -- -5.542162e+16 -2.977163e+16
1.60   3.486784401000e-11  -2.267025187393e-03  -9.986696039732e-04 -- -6.501765e+07 -2.864157e+07
2.00   7.979226629761e-04  -2.498599886894e-02  -6.362375005551e-03 -- -3.231381e+01 -8.973674e+00
5.00   2.312248366666e+11   2.312248356551e+11   2.312248368273e+11 -- -4.374473e-09 6.947283e-10
h:
x, "exact", naive, Horner -- relative errors
0.10  -0.000000000000e+00   3.739352465198e+04   4.481800000000e+04 --          -          -
0.20  -8.460014212631e+15  -8.460014212595e+15  -8.460014212587e+15 -- -4.224343e-12 -5.158266e-12
0.30  -1.154825498147e+16  -1.154825498143e+16  -1.154825498142e+16 -- -3.410039e-12 -3.675014e-12
0.50  -9.999185058306e+15  -9.999185058245e+15  -9.999185058266e+15 -- -6.050493e-12 -4.009927e-12
1.01  -7.137691193017e+14  -7.137691188577e+14  -7.137691192674e+14 -- -6.220714e-10 -4.806876e-11
1.20   4.923817795711e+14   4.923817806832e+14   4.923817796032e+14 -- 2.258555e-09 6.511969e-11
1.30   7.371226583915e+14   7.371226601609e+14   7.371226584225e+14 -- 2.400446e-09 4.204643e-11
1.40   8.035389574478e+14   8.035389602065e+14   8.035389574776e+14 -- 3.433108e-09 3.705436e-11
1.60   6.341259826947e+14   6.341259889764e+14   6.341259827225e+14 -- 9.906112e-09 4.386116e-11
2.10  -5.923385583628e-02   3.500333996906e+07   2.383600000000e+04 -- -5.909347e+08 -4.024060e+05
5.10  -3.774706499371e-03   8.557366072250e+09   3.909520000000e+05 -- -2.267028e+12 -1.035715e+08
order (n), Pn(1), Pn(-1)
expected: 1 and 1
2                    1                    1
4        1.00000000002        1.00000000002
6        0.99999999994        0.99999999994
8        0.99999999992        0.99999999992
10        0.99999999907        0.99999999907
12        1.00000000349        1.00000000349
14        1.00000000299        1.00000000299
16       0.999999939378       0.999999939378
18         1.0000000112         1.0000000112
20        1.00000501373        1.00000501373
22        1.00001560093        1.00001560093
24       0.999974806438       0.999974806438
26        1.00033043498        1.00033043498
28        1.00687428185        1.00687428185
30        1.04755829826        1.04755829826
32        1.13065619275        1.13065619275
34        1.29242497941        1.29242497941
36        1.53122404724        1.53122404724
38        10.4985706725        10.4985706725
40        50.6674061331        50.6674061331
42       -504.307776538       -504.307776538
44       -971.900395524       -971.900395524
46        7574.82722195        7574.82722195
48        111470.061522        111470.061522
50          500314.6215          500314.6215```

The overall conclusions:

• Horner's rule is in many cases (slightly) more accurate, but can be way off the correct value (see for instance: h(2.1), which should have been zero - but even the "exact" value is definitely not zero).
• Evaluating Legendre polynomials should not be done with Horner's rule when the order gets over, say, 10 or 20.

TV (13 jul 04) Em, urh, well, the remarks starting of this page suggest things I don't really understand.

What's wrong with filling in a simple (for engineers)

```  N
----
\          i
>   a . x
/      i
----
i=0```

For any moderate N? That's a matter of usually floating point or maybe rationals or in some cases integer precision, for which such a formula has a normal condition, nothing special, unless maybe N is large or X is large or the sheer number of i's makes the mantissa not acurate enough. I don't see the problem.

With respect to approximations by polynoms, I don't see your point, but I do know what an important point can be: using polynomal or power basis to approximate functions, direct polynomal approximation like in a spline approach can start oscilating and becoming unstable when the number of terms starts to grow. Higher power point-evaluation based polynomials are not good idea, and often of course as always it will depend on your problem: most natural problems do not by nature *have* higher orders in them, so trying to approximate by using higher powers is doomed to be inacurate by nature, or leave you with an approximation where you're trying to find higher order terms which you'd basically find to be zero.

As engineer one would know of course that one of the by far preferable approximation methods for many practical problems, and of great theoretical value, which in a sense can be seen to have polynomial results, is the Taylor Expansion.

I've mentioned this expansion being available in Maxima on the page about that package.

With respect to numerical accuracy in general: according to (as I remember vividly from a 1st year physics practicum) normal error propagation rules, one can compute what the maximum error is for a certain evaluation of a given formula, and then figure out if the number of bits in the mantissa of the floating point operations is sufficient. Doing that analysis to some depth would seem intesting, but then one would probably start with integers, and include a decent rounding analyis.

Category Numerical Analysis

]