Critical points of a polynomial function

Hi everybody,

I would like to know the Matlab code to obtain the value of x at the maximum "y" value for a 3rd order polynomial. Basically, I have to obtain the area under the curve of a plot by integrating from 0 to the "x" value at the maximum "y" point.

I will appreciate any suggestion,

Thanks

Sergio

Answer 1

[Dr. Seis]

Does something like this work:

x_range = [0,10]; 

p = [1,4,-2,1]; % the 3rd order polynomial coefficients for (x^3+4*x^2-2*x+1)
k = polyder(p); % find derivative of above polynomial
r = roots(k); % find roots of polynomial created from derivative above
r = r(r>=x_range(1)); % get rid of roots less than x_min
if ~isempty(r)
  r = r(r<=x_range(2)); % get rid of roots greater than x_max
end

if ~isempty(r)
  try_points = zeros(1,2+length(r)); % pre-allocate
  try_points(1:2) = x_range; % add x_min and x_max to list to try
  for ii = 1 : length(r)
      try_points(2+ii) = r(ii);
  end
else
  try_points = x_range;
end

y = polyval(p,try_points);
[max_y,xi] = max(y);

peak_at_xy = [try_points(xi),max_y]; % x & y coordinates of peak

This program basically takes the derivative of your polynomial (coefficients defined in "p") to find the peaks a valleys, then it compiles a list of "x" values inside your range in "x" (i.e., "x_range"). It then computes y for that list of "x" values to try and finds the "x" that gives the maximum "y" on the interval defined by "x_range". Of course, you would define the maximum "x" in "x_range" to be whatever you want.

[Sergio]

Thanks for your soon answer.

I also need to integrate the function from 0 until the "x" value corresponding to the maximum "y" value. Do you know the matlab code for the integration with limits ?. Finally, I would like to know the coefficient of correlation of the trend line that fits better the x and y database. I sometimes increase the order of the polynomial to obtain a best fit but I also want to know the coefficient of correlation of the equation.

I will appreciate any help you can provide me.

Thanks,

Sergio

[Dr. Seis]

As far as the integration function, you should check out the help pages for "quad" - just above the section "Algorithm" describes a way of passing an anonymous function handle [i.e., F = @(x)x.^3+4*x.^2-2*x+1 in the example above] to the "quad" function over some integration limits.

Then is would be as simple as Q = quad(F,0,try_points(xi)) after defining your function "F" and finding the appropriate "x" value.

As far as finding the best fit polynomial coefficients (and some statistics on how well fit those coefficients are) you should check out "polyfit" in the help pages.

In the future, though, you should ask separate questions on separate pages to increase the number of people who might see it and be able to answer it better.

[Mak Dukan]

Hi all,

I am trying to understand better this code.

My questions are both mathematical and technical. First for the technical:

1) "r = r(r>=x_range(1)); % get rid of roots less than x_min": How does this code get rid of these roots? Why do we "write x_range(1)" and not some other number?

2) "r = r(r<=x_range(2)); % get rid of roots greater than x_max": Same question as above

3) Why do we first try for point 1, then point 2 and then for the entire lenght? Why not immediatrly try for entire lenght?

Then the mathematical questions (in danger of sounding naive):

1) Would this code look any different if we already know that our function does not have any other max besides one single point? In my case, I know how the function looks like in Excel, but want to derive its max point in matlab because I think it is much easier for multiple calculations.

2) I suppose that we by all means need to make the derivative of the polynomial first? I mean the derivative gives us the area under the curve, but it also gives us the max point? Right? This is why we looks for the roots or multiple possible solutions?

Thanks (am trying to understand what I am doing, as opposed to just copy pasting code)!

Mak

[Mak Dukan]

In addition:

in my case I know the exact x and y values of a line and want to use matlab to give me a max point on the curve. I have modified the above code for my purpose, and it works:

clear all

%read merchant and CfD Shareholder value (2) from an Excel file
filename = 'M:\Mak Folder\PAPERS\Paper C\#### EXCEL MAIN ####\Mod_(44)_R002_Base_FitF.xlsm';
CfD_Sv = xlsread(filename,'Overview','C14:M14');
Mer_Sv = xlsread(filename,'Overview','C15:M15');

%make values for debt levels
DL = [0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9];

x_range = [0.4,0.9]; 

CfD_p = polyfit(DL,CfD_Sv,3);
Mer_p = polyfit(DL,Mer_Sv,3);

%find max for CfD

k = polyder(CfD_p); % find derivative of above polynomial
r = roots(k); % find roots of polynomHHial created from derivative above
r = r(r>=x_range(1)); % get rid of roots less than x_min
if ~isempty(r)
  r = r(r<=x_range(2)); % get rid of roots greater than x_max
end
if ~isempty(r)
  try_points = zeros(1,2+length(r)); % pre-allocate
  try_points(1:2) = x_range; % add x_min and x_max to list to try
  for ii = 1 : length(r)
      try_points(2+ii) = r(ii);
  end
else
  try_points = x_range;
end
y = polyval(CfD_p,try_points);
[max_y,xi] = max(y);
peak_at_xy = [try_points(xi),max_y]; % x & y coordinates of peak

Would be great to hear some comments to improve this? Maybe I am doing something redundant?