derivative of a polynomial function in C and ...

Papabravo

Joined Feb 24, 2006
21,157
I don't think he is wanting to do numeric estimation of the derivative. But even if he is, setting the delta to that value is inviting disaster. If x is greater than 1 then you will get zero for the numerator.
I may have misinterpreted the meaning of:
"...that compute the derivative of a polynomial function", and
"... I am also trying to look at the analytical method that compute the derivative of a function."
It is easy to conflate analytical with numerical on the first reading pass given the context.

It depends on how you perform the operations, and how many significant digits the algorithms use for computation versus displaying results. What makes you think this is a problem? Can you show an example? My spreadsheet is just fine with the method using values of x greater than 1. It is also true that real polynomials with positive integral coefficients diverge to ±∞, but that is the same concern regardless of the method you use to compute the value of the derivative.
 
Last edited:

WBahn

Joined Mar 31, 2012
29,976
For any floating point representation, there is a limited precision. This is usually designated by a value "epsilon" (often shortened to "eps") which is the smallest positive value such that (1+eps) is different than 1. For IEEE 754 single-precision floating point values, eps ~= 1.192e-7. So if you evaluate f(x) and f(x+dx) where x = 1.0 and dx = 1.0e-7 using floats, you will get the same result for both evaluations because, in both cases, your argument is 1.0. For IEEE-754 double-precision representation, epx ~= 2.220e-16. So if you use dx = 1.0e-15 with x = 1.0 you will not have zero in the numerator, but you will have significant roundoff error because the actual width of your slice and the assumed width of your slice differ significantly.

To see the effect, let's use f(x) = x. We know that df(x)/dx is exactly 1 independent of the choice of x or dx. So let's see what your recommendation yields.

C:
#include <stdio.h>

#define DX_FLT ((float)1e-7)
#define DX_DBL ((double)1e-15)

float f_flt(float x)
{
   return x;
}

double f_dbl(double x)
{
   return x;
}

int main(void)
{
   double x_flt, dfdx_flt;
   double x_dbl, dfdx_dbl;
 
  double x[] = {0.0, 0.001, 0.002, 0.005,
                     0.01,  0.02,  0.05, 
                     0.1,   0.2,   0.5,   
                     1.0,   2.0,   5.0,   
                    10.0,  20.0,  50.0, 
                   100.0};

   int i;

   for (i = 0; i < sizeof(x)/sizeof(*x); i++)
   {
      x_flt = x[i];
      x_dbl = x[i];
      dfdx_flt = ( f_flt( x_flt + DX_FLT) - f_flt(x_flt) ) / DX_FLT;
      dfdx_dbl = ( f_dbl(x_dbl + DX_DBL) - f_dbl(x_dbl) ) / DX_DBL;
      printf("%15.10f  %15.10f  %15.10f\n", x[i], dfdx_flt, dfdx_dbl);
   }
 
   getch();
   return 0;
}
Here is the output:

Code:
  0.0000000000  1.0000000000  1.0000000000
  0.0010000000  1.0000076028  1.0000680839
  0.0020000000  0.9988434496  1.0000680839
  0.0050000000  1.0011717560  1.0000680839
  0.0100000000  1.0058283688  0.9992007222
  0.0200000000  1.0058283688  0.9992007222
  0.0500000000  1.0058283688  0.9992007222
  0.1000000000  0.9685754663  0.9992007222
  0.2000000000  1.0430812714  0.9992007222
  0.5000000000  1.1920928816  0.9992007222
  1.0000000000  1.1920928816  1.1102230246
  2.0000000000  0.0000000000  0.8881784197
  5.0000000000  0.0000000000  0.8881784197
 10.0000000000  0.0000000000  1.7763568394
 20.0000000000  0.0000000000  0.0000000000
 50.0000000000  0.0000000000  0.0000000000
100.0000000000  0.0000000000  0.0000000000
As you can see, the value of eps is readily apparent in the float data (it's a bit masked in the double data, but still very evident). But even this is due to pronounced roundoff error. If you print out f(x+dx) - f(x) you should get 1.0e-7 but, instead, you get 1.192... e-7 due to the limited precision.

The ONLY value that produces correct results is at x=0 and that's because the difference yields full precision and doesn't suffer from the kind of extreme roundoff error that even "small" values of x do.

Even at x=1 there is nearly a 20% error in the calculation with floats and more than a 10% error in the calculation with doubles.

As expected, as soon as you get much above 1.0, things go to zero for the float and they go to zero for the double just above 10.0.
 

WBahn

Joined Mar 31, 2012
29,976
Since you said you were using a spreadsheet, here's an example using a spreadsheet (MS Excel):

x| 1.00000000000000000 |
dx| 0.00000000000001000 |(1e-14)
(x+dx) |1.00000000000001000|
(x+dx)-x| 0.00000000000000999 |
[(x+dx)-dx]/dx |0.99920072216264100 |
||
x |1.00000000000000000 |
dx |0.00000000000000100| (1e-15)
(x+dx) |1.00000000000000000 |
(x+dx)-x |0.00000000000000000|
[(x+dx)-dx]/dx |0.00000000000000000 |
 

Papabravo

Joined Feb 24, 2006
21,157
I hadn't considered this particular example so thanks for the illustration that even simple problems can present unexpected and counter intuitive behavior. I was aware of the problem with floating point addition and subtraction concerning numbers with different exponents. The process of de-normalization, to allow the operation to be performed, causes problems with retaining enough significant bits in both the mantissas prior to the addition or subtraction.
 

WBahn

Joined Mar 31, 2012
29,976
Yeah, it can be a big issue that sneaks up on you. It's amazing when you think about it. The interatomic spacing in most solids is on the order of 10^-10 m. So 15 sig figs is enough to resolve the distance between two atoms along a length that is more than twice the circumference of the earth. Yet there are numerous real world examples where that has not been sufficient resolution to prevent unacceptable roundoff errors from happening.

One of the big reasons that we need all those sig figs is precisely so that we can add/subtract values of different orders of magnitude and preserve enough sig figs in the result to still be useful.
 

Thread Starter

Eric007

Joined Aug 5, 2011
1,158
Be sure to think about the limitations of that approach.
Yes! I just wanted to test the functionality of something. In my code I have limited the exponent to 3 and I am not dealing with negative exponent and / or exponents with decimal point.

From a practical standpoint, what if the person wants to enter f(x) = x^27 - 1 ?

Then consider if you want to limit yourself to strict polynomial functions requiring that the exponents be a set of the whole numbers (non-negative integers). If you do (and if the order of the polynomial is going to be reasonably small) then this is a good way to go.
Ok. This could actually become a nice exercise where you involve all cases, adding complexity by, for instance, adding multiple variables x, y, z...:)
 

Thread Starter

Eric007

Joined Aug 5, 2011
1,158
Again, it comes down to how you want to represent the information.
Yep! I responded too quickly without actually reading your comment properly.

If you keep an array of the coefficients then you can compute the coefficients of the derivative very easily.

You have a couple of options. The first is to maintain an array in which the index is the exponent -- though this only works if the exponents are non-negative integers. But if the order of your polynomial is large and/or arbitrary, then it might be better to maintain parallel arrays with the coefficients in one and the exponents in the other.

If you want the person to enter the polynomial, then you have to decide how you want them to enter it. They could enter it as two lists containing coefficients in one and exponents in the other. Another would be to provide a list of coefficient/exponent pairs. Yet a third would be to enter the polynomial as a string and then parse it to generate your internal representation.
Like I said in the previous post, I think parsing would be the best option. If I find some time I will write the code, as an exercise, where the user enter the polynomial as a string then the program parse it so we won't have a problem for cases like f(x) = x^27 - 1 then I will add cases one at a time...

Also for the option where coefficients would be stored in an array, the array should be dynamic...

Thanks guys for all the comments... will read the rest soon...
 

Thread Starter

Eric007

Joined Aug 5, 2011
1,158
What I was actually doing is writing a simple code that finds the root of a polynomial using Newton method. The code worked perfect so no need to post a code but I do have a question though.

The attached equation will be iterated a number of times, let's say 0 <= i < n-1 . We either find the closest root or the sequence does not converge to a root of f(x).
Now my question is what can you conclude about the chosen x0 when f '(x0) in the attached equation yields 0 at iteration i = 0.

I just display a message saying the value cancel the first derivative of f(x).
For example f(x) = - x^2 + 2x with x0 = 1
 

Attachments

Last edited:

WBahn

Joined Mar 31, 2012
29,976
If f'(x)=0, then you can conclude that f(x) is constant and not a function of x.
He's not saying f'(x) = 0 for all x, he's saying that f'(xo) for a specific value of x0.

What you can say is that this is a local (possibly global) min or max. You can't say much about where any roots lie relative to this point.
 

vpoko

Joined Jan 5, 2012
267
He's not saying f'(x) = 0 for all x, he's saying that f'(xo) for a specific value of x0.

What you can say is that this is a local (possibly global) min or max. You can't say much about where any roots lie relative to this point.
Ahh, you are correct as always.

Eric007, the example polynomial you gave is quadratic, which means the corresponding graph would be parabolic. It has exactly one point where the slope of the tangent line will be zero (horizontal) and doesn't point towards the x axis.
 

WBahn

Joined Mar 31, 2012
29,976
Ahh, you are correct as always.
Ah, but were that true...


One thing that I said that bothered me at the time, but I needed to come back to it and convince myself that I was wrong (or at least not complete) is saying that if the derivative at a point is zero that it is a local (and possibly global) minima or maxima. This isn't true. While it is necessary that the derivative be zero at points that are local extrema, it is not sufficient. The function could be either increasing or decreasing at that point.

To prove that to myself I set about constructing such a function and it was trivially easy to do so. I picked a cubic function:

\(
y(x) \: = \: x^3 \: + \: ax^2 \: + \: bx \: + \: c
\)

And asked what would happen if the two points were the curve's derivative go to zero were set equal to each other.

That merely requires that we set

\(
b = \frac{a^2}{3}
\)

Doing so results in the derivative going to zero at

\(
\frac{dy}{dx} \: = \: 0 \: @ \: x \: = \: - \frac{a}{3}
\)

The derivative is positive on both sides of this point, so the function is continuously increasing and, hence, has no local maxima at all (it has a global min at -∞ and a global max at +∞ and that's it).
 

vpoko

Joined Jan 5, 2012
267
Yes, if I remember my differential calculus, it's known as a critical point. All extrema are at critical points, but not all critical points are extrema, they can also be inflection points.
 

WBahn

Joined Mar 31, 2012
29,976
I don't recall them being called "critical points", but that was many moons ago and those are the kinds of details that I forget pretty easily. Still, I'm curious if my calc book called them that or something else. I'll have to try to remember to look.
 

Papabravo

Joined Feb 24, 2006
21,157
I remember both terms, "critical points" and "inflection point" being used.

Thompson, S.P., Calculus Made Easy, 1910, uses "inflection point". The specific example is \(y=x^3 \ at \ x=0\)

The MIT online course: Strang, G., Calculus
http://ocw.mit.edu/ans7870/resources/Strang/Edited/Calculus/Calculus.pdf


Uses the term "critical point" on p. 99

Important The maximum always occurs at a stationay point (where dfldx = 0) or a rough point (no derivative) or an endpoint of the domain. These are the three types of critical points. All maxima and minima occur at critical points!

And uses the term "inflection point" on p.106 in the caption to Fig. 3-7

Fig. 3.7 Increasing slope =concave up (f" >0). Concave down is f" <0. Inflection point f" = 0.

The last point about the second derivative being 0 at an inflection point was the one I'd forgotten.
 

WBahn

Joined Mar 31, 2012
29,976
And a similar caveat for second derivatives and inflection points applies. All inflection points occur at places where the second derivative vanishes, but not all places where the second derivative vanishes are inflection points. The function y = x^4 is an example. The function's first and second derivatives are zero at x=0, yet the function lacks an inflection point altogether.

In general, the information provided by a given order of derivative at points where it is non-zero is absent at points where it IS zero without examining the next higher-order derivative at that point -- and if the next higher-order derivative is zero then you have to go up again and keep going until you get a non-zero derivative.
 

MrAl

Joined Jun 17, 2014
11,388
Yeah, it can be a big issue that sneaks up on you. It's amazing when you think about it. The interatomic spacing in most solids is on the order of 10^-10 m. So 15 sig figs is enough to resolve the distance between two atoms along a length that is more than twice the circumference of the earth. Yet there are numerous real world examples where that has not been sufficient resolution to prevent unacceptable roundoff errors from happening.

One of the big reasons that we need all those sig figs is precisely so that we can add/subtract values of different orders of magnitude and preserve enough sig figs in the result to still be useful.
Hi,

That's a very interesting way to look at it. It's almost ironic or something that we still have numerical problems in calculations these days.

For whomever is interested, the numerical derivative:
dy/dx=(f(x+h)-f(x))/h

is often called a "one sided derivative", and the other is:
dy/dx=(f(x)-f(x-h))/h

A slightly better one is the "central means derivative":
dy/dx=(f(x+h)-f(x-h))/(h+h)

which i think is second order whereas the first two are just first order. h+h is used to avoid multiplication.
There are higher order numerical derivatives, but they are harder to generate. I'd actually have to look them up again as to how to generate them because it's been a while since i needed them.
 

joeyd999

Joined Jun 6, 2011
5,234
A slightly better one is the "central means derivative":
dy/dx=(f(x+h)-f(x-h))/(h+h)
Funny. I've never heard of this, but it makes sense. Are there classes of equations where this is a better solution?

Please point me to an authoritative reference -- Google seems not to know of this either.
 
Top