www.pudn.com > starter.zip > polyinterp.m, change:2011-01-04,size:4217b


function [minPos,fmin] = polyinterp(points,doPlot,xminBound,xmaxBound) 
% function [minPos] = polyinterp(points,doPlot,xminBound,xmaxBound) 
% 
%   Minimum of interpolating polynomial based on function and derivative 
%   values 
% 
%   In can also be used for extrapolation if {xmin,xmax} are outside 
%   the domain of the points. 
% 
%   Input: 
%       points(pointNum,[x f g]) 
%       doPlot: set to 1 to plot, default: 0 
%       xmin: min value that brackets minimum (default: min of points) 
%       xmax: max value that brackets maximum (default: max of points) 
% 
%   set f or g to sqrt(-1) if they are not known 
%   the order of the polynomial is the number of known f and g values minus 1 
 
if nargin < 2 
    doPlot = 0; 
end 
 
nPoints = size(points,1); 
order = sum(sum((imag(points(:,2:3))==0)))-1; 
 
% Code for most common case: 
%   - cubic interpolation of 2 points 
%       w/ function and derivative values for both 
%   - no xminBound/xmaxBound 
 
if nPoints == 2 && order ==3 && nargin <= 2 && doPlot == 0 
    % Solution in this case (where x2 is the farthest point): 
    %    d1 = g1 + g2 - 3*(f1-f2)/(x1-x2); 
    %    d2 = sqrt(d1^2 - g1*g2); 
    %    minPos = x2 - (x2 - x1)*((g2 + d2 - d1)/(g2 - g1 + 2*d2)); 
    %    t_new = min(max(minPos,x1),x2); 
    [minVal minPos] = min(points(:,1)); 
    notMinPos = -minPos+3; 
    d1 = points(minPos,3) + points(notMinPos,3) - 3*(points(minPos,2)-points(notMinPos,2))/(points(minPos,1)-points(notMinPos,1)); 
    d2 = sqrt(d1^2 - points(minPos,3)*points(notMinPos,3)); 
    if isreal(d2) 
        t = points(notMinPos,1) - (points(notMinPos,1) - points(minPos,1))*((points(notMinPos,3) + d2 - d1)/(points(notMinPos,3) - points(minPos,3) + 2*d2)); 
        minPos = min(max(t,points(minPos,1)),points(notMinPos,1)); 
    else 
        minPos = mean(points(:,1)); 
    end 
    return; 
end 
 
xmin = min(points(:,1)); 
xmax = max(points(:,1)); 
 
% Compute Bounds of Interpolation Area 
if nargin < 3 
    xminBound = xmin; 
end 
if nargin < 4 
    xmaxBound = xmax; 
end 
 
% Constraints Based on available Function Values 
A = zeros(0,order+1); 
b = zeros(0,1); 
for i = 1:nPoints 
    if imag(points(i,2))==0 
        constraint = zeros(1,order+1); 
        for j = order:-1:0 
            constraint(order-j+1) = points(i,1)^j; 
        end 
        A = [A;constraint]; 
        b = [b;points(i,2)]; 
    end 
end 
 
% Constraints based on available Derivatives 
for i = 1:nPoints 
    if isreal(points(i,3)) 
        constraint = zeros(1,order+1); 
        for j = 1:order 
            constraint(j) = (order-j+1)*points(i,1)^(order-j); 
        end 
        A = [A;constraint]; 
        b = [b;points(i,3)]; 
    end 
end 
 
% Find interpolating polynomial 
params = A\b; 
 
% Compute Critical Points 
dParams = zeros(order,1); 
for i = 1:length(params)-1 
    dParams(i) = params(i)*(order-i+1); 
end 
 
if any(isinf(dParams)) 
    cp = [xminBound;xmaxBound;points(:,1)].'; 
else 
    cp = [xminBound;xmaxBound;points(:,1);roots(dParams)].'; 
end 
 
% Test Critical Points 
fmin = inf; 
minPos = (xminBound+xmaxBound)/2; % Default to Bisection if no critical points valid 
for xCP = cp 
    if imag(xCP)==0 && xCP >= xminBound && xCP <= xmaxBound 
        fCP = polyval(params,xCP); 
        if imag(fCP)==0 && fCP < fmin 
            minPos = real(xCP); 
            fmin = real(fCP); 
        end 
    end 
end 
% Plot Situation 
if doPlot 
    figure(1); clf; hold on; 
 
    % Plot Points 
    plot(points(:,1),points(:,2),'b*'); 
 
    % Plot Derivatives 
    for i = 1:nPoints 
        if isreal(points(i,3)) 
            m = points(i,3); 
            b = points(i,2) - m*points(i,1); 
            plot([points(i,1)-.05 points(i,1)+.05],... 
                [(points(i,1)-.05)*m+b (points(i,1)+.05)*m+b],'c.-'); 
        end 
    end 
 
    % Plot Function 
    x = min(xmin,xminBound)-.1:(max(xmax,xmaxBound)+.1-min(xmin,xminBound)-.1)/100:max(xmax,xmaxBound)+.1; 
    size(x) 
    for i = 1:length(x) 
        f(i) = polyval(params,x(i)); 
    end 
    plot(x,f,'y'); 
    axis([x(1)-.1 x(end)+.1 min(f)-.1 max(f)+.1]); 
 
    % Plot Minimum 
    plot(minPos,fmin,'g+'); 
    if doPlot == 1 
        pause(1); 
    end 
end