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. ---- [[ [Category Mathematics] | [Category Numerical Analysis] ]]