File Exchange Pick of the Week

August 26th, 2008

This puzzler is very straightforward, so I hope to hear from some of the newer MATLAB users. I think the solution is very linear, most people would come up with the same solution. As a hint, be sure to look at our help for POLYDER and ROOTS.

The challenge is this: Given the coefficients of two second order polynomials, find the local extrema of the first one. Knowing the X where this point happens, put a marker on both curves at that X value. Plot both curves for three units to the left and right of the inflection point of the first curve. [Edited Aug 27th in response to J. Paul R’s observation that I misused terminology, switching inflection point for local extrema. MATLAB t-shirt for the catch. Thank you!]

Here is the code to create the coefficients:

coef1 = rand(1,3)-0.5;
coef2 = rand(1,3)-0.5;

<pre> <code>

all the code so someone can just copy and paste it from the comments.

</code> </pre>

08:36 UTC | Posted in MATLAB Basics, Puzzler | Permalink |

You can follow any responses to this entry through the RSS 2.0 feed. You can skip to the end and leave a response. Pinging is currently not allowed.

11 Responses to “Puzzler: Working with polynomials”

1. Dan Christensen replied on August 26th, 2008 at 09:59 UTC :
   

   Here’s my approach using matrices for the calculation

    
   
   der1 = polyder(coef1);
   infl = roots(der1);
   
   x = linspace(infl-3,infl+3);
   
   Y =  [coef1 ; coef2] * [x.^2 ; x ; ones(1,length(x))];
   yi = [coef1 ; coef2] * [infl.^2 ; infl ; 1];
   
   figure
   plot(x,Y(1,:),'r')
   hold on
   plot(x,Y(2,:),'k')
   plot(infl,yi(1),'ks')
   plot(infl,yi(2),'k^')
   
    

2. J.Paul R replied on August 26th, 2008 at 11:16 UTC :
   

   The real puzzler is: “what is this puzzler?” Your text ask for an “inflection point” while your example shows an extrema. Additional, 3 coefficients specify a 2nd degree polynomial which has no inflection points. So are you looking for extrema or inflection points? Or perhaps the extrema of the first derivative, which is, generally, an inflection point of the original polynomial?

3. J.Paul R replied on August 26th, 2008 at 11:23 UTC :
   

   BTW, here’s my code for inflection points. I don’t think there are any edge cases to worry about with 3rd order polynomials. (Like x=0 for x^4)

    
   coef1 = rand(1,4)-0.5;
   coef2 = rand(1,4)-0.5;
   
   der2 = polyder(polyder(coef1));
   infpt = roots(der2);
   
   x = linspace(infpt-3,infpt+3);
   
   Y1 = polyval(coef1,x);
   Y2 = polyval(coef2,x);
   yi1 = polyval(coef1,infpt);
   yi2 = polyval(coef2,infpt);
   
   figure(1)
   clf
   plot(x,Y1,'r')
   hold on
   plot(x,Y2,'k')
   plot(infpt,yi1,'ks')
   plot(infpt,yi2,'k^')
   
    

4. Danilo replied on August 26th, 2008 at 15:07 UTC :
   

   coef1 = rand(1,3)-0.5;
   coef2 = rand(1,3)-0.5;
   
   infl = roots(polyder(coef1));
   
   xx = linspace(infl-3, infl+3);
   xx(50) = infl;
   y1 = polyval(coef1, xx);
   y2 = polyval(coef2, xx);
   
   plot(xx, y1, 'r', xx, y2, 'k',...
       infl, y1(50), 'sk', infl, y2(50), '^k');

5. Mike Briggs replied on August 27th, 2008 at 12:15 UTC :
   

   OK, here is my totally AR* comment: extrema is plural; for a quadratic, there can be only one local extremum (the singular version of extrema).

   Yours in good grammar, Mike Briggs (-:

   * AR: anal-retentive

6. Luca Balbi replied on August 28th, 2008 at 02:31 UTC :
   

   A 2nr order polynomial is just the old parabola from high school. It extremum is called the vertex and has analytically known coordinates. Its x coordinate is just

   
   (-coef1(2)/(2*coef1(1)))
   

   My solution is the following:

   
   coef1 = rand(1,3)-0.5;
   coef2 = rand(1,3)-0.5;
   
   extremum = (-coef1(2)/(2*coef1(1)));
   
   x = linspace(extremum-3,extremum+3);
   
   Y =  [coef1 ; coef2] * [x.^2 ; x ; ones(1,length(x))];
   yi = [coef1 ; coef2] * [extremum.^2 ; extremum ; 1];
   
   figure
   plot(x,Y(1,:),'r')
   hold on
   plot(x,Y(2,:),'k')
   plot(extremum,yi(1),'ks')
   plot(extremum,yi(2),'r^')
   

7. david replied on August 28th, 2008 at 03:19 UTC :
   

   perhaps some error checking is in order. After all, it is possible that our randomly generated quadratic happens to be linear ( coef1(1) == 0 ). In this case, the first three codes above codes throw a cryptic error “Inner matrix dimensions must agree.” from which it is not obvious what has gone wrong (roots(polyder) returns an empty matrix). The last code runs without error, but produces nothing as extremum has value +/-Inf

   Since this error is fairly unlikely, simply running the code a bunch of times with random input probably won’t notice it. What caught my attention was the last code, which happily sticks a random number from an interval containing zero at the bottom of a fraction.

   A quick check to make we have a quadratic before diving in to extreme point finding may be in order:

    
   if (coef1(1) == 0)
      error('Error: Generated polynomial is not quadratic!');
   end
    

   Regards,
   David

8. Luca Balbi replied on August 28th, 2008 at 07:41 UTC :
   

   While we’re at it…
   Checking for a number to be zero is tricky in itself.

   We’re better off using

   
   if (abs(coef1(1)) < eps)
      error('Error: Generated polynomial is not quadratic!');
   end
   

   revised code:

   
   coef1 = rand(1,3)-0.5;
   coef2 = rand(1,3)-0.5;
   
   if (abs(coef1(1)) < eps)
      error('Error: Generated polynomial is not quadratic!');
   end
   
   extremum = (-coef1(2)/(2*coef1(1)));
   
   x = linspace(extremum-3,extremum+3);
   
   Y =  [coef1 ; coef2] * [x.^2 ; x ; ones(1,length(x))];
   yi = [coef1 ; coef2] * [extremum.^2 ; extremum ; 1];
   
   figure
   plot(x,Y(1,:),'r')
   hold on
   plot(x,Y(2,:),'k')
   plot(extremum,yi(1),'ks')
   plot(extremum,yi(2),'r^')
   

   by the way, infinity is the correct answer in the limit of the 2nd degree coefficient going to zero, as long as the first degree coefficient is not zero: in that case infinity is not the solution either.

9. Tareq replied on August 28th, 2008 at 08:06 UTC :
   

   coef1 = rand(1,3)-0.5;
   coef2 = rand(1,3)-0.5;

   lex1=roots(polyder(coef1))
   lex2=roots(polyder(coef2))

   hold all
   plot(lex1-3:0.1:lex1+3,polyval(coef1,lex1-3:0.1:lex1+3),’-');
   plot(lex1-3:0.1:lex1+3,polyval(coef2,lex1-3:0.1:lex1+3),’-');

   plot(lex1,polyval(coef1,lex1),’o');
   plot(lex1,polyval(coef2,lex1),’d');

   xlim([lex1-3 lex1+3])

10. Frank replied on August 30th, 2008 at 09:44 UTC :
    

    This is my solution:

     
    
    coef1 = rand(1,3)-0.5;
    coef2 = rand(1,3)-0.5;
    
    der1 = polyder(coef1);
    ex1 = roots(der1);
    
    x = linspace(ex1-3,ex1+3);
    
    curve1 = polyval(coef1,x);
    curve2 = polyval(coef2,x);
    
    marker1 = polyval(coef1,ex1);
    marker2 = polyval(coef2,ex1);
    
    plot(x,curve1,x,curve2,ex1,marker1,'ks',ex1,marker2,'k^')
    
     

11. Doug replied on September 2nd, 2008 at 09:36 UTC :
    

    Wow!

    I have been out for a week, so I have not been able to respond to all the great comments here. I find it interesting how the challenge evolved.

    First many people solve the problem and then people start to see edge cases where the square term is zero. Then people find that the test for this edge case can be found more effectively in a different way.

    This kind of hardening of the code is exactly why we have frequent “code reviews” at The MathWorks to make sure many eyes see code before it enters the testing stage and finally the product itself.

    Thanks everyone!
    Doug

Leave a Reply

<pre class="code">
a = magic(3);
sum(a)
</pre>

These postings are the author's and don't necessarily represent the opinions of The MathWorks.