function varargout = newton_stability(varargin)
% NEWTON_STABILITY MATLAB code for newton_stability.fig
%      NEWTON_STABILITY, by itself, creates a new NEWTON_STABILITY or raises the existing
%      singleton*.
%
%      H = NEWTON_STABILITY returns the handle to a new NEWTON_STABILITY or
%      the handle to
%      the existing singleton*.
%
%      NEWTON_STABILITY('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in NEWTON_STABILITY.M with the given input arguments.
%
%      NEWTON_STABILITY('Property','Value',...) creates a new NEWTON_STABILITY or raises the
%      existing singleton*.  Starting from the left, property value pairs
%      are
%      applied to the GUI before newton_stability_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to newton_stability_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help newton_stability

% Last Modified by GUIDE v2.5 08-Aug-2012 09:31:57

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @newton_stability_OpeningFcn, ...
                   'gui_OutputFcn',  @newton_stability_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before newton_stability is made visible.
function newton_stability_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to newton_stability (see VARARGIN)

% Choose default command line output for newton_stability
handles.output = hObject;

% if we pass a command line argument, this is the initial number of nodes
% to use

if( length(varargin) == 1 )
    nodecount=varargin{1};
    hSlider=handles.slider1;
    minv=get(hSlider,'Min');
    maxv=get(hSlider,'Max');

    scamin=1; % scaled minimum
    scamax=50; % scaled maximum

    raw=(nodecount-scamin)*(maxv-minv)/(scamax-scamin)+minv;
    
    set(hSlider,'Value',raw);
end

    

% Get slider position and update the label accordingly
getSliderPosition(handles.slider1);
refreshSliderLabel(handles.slider1,handles.text2);

updatePlot(handles);

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes newton_stability wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = newton_stability_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on slider movement.
function slider1_Callback(hObject, eventdata, handles)
% hObject    handle to slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

getSliderPosition(hObject);
refreshSliderLabel(hObject,handles.text2);

updatePlot(handles);






% --- Executes during object creation, after setting all properties.
function slider1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end


function updatePlot(handles)
ud=get(handles.slider1,'UserData');

nodecount=ud.nodecount;

f=@(x)(sin(pi*x));

x=linspace(0,0.2,nodecount);
y=f(x);
xx=linspace(0,0.2,1000);
yy=f(xx);
yyn=newtonpoly(x,y,xx);
yyl=lagrangepoly(x,y,xx);
plot(handles.axes1,xx,log10(max(eps,abs(yyn-yy))),'-r',xx,log10(max(eps,abs(yyl-yy))),'-b')%,xx,log10(max(eps,abs(yyl2-yy))),'c+');
legend({'Newton','Lagrange'},'Location','SouthOutside'); %,'Lagrange (lagrange)'});
xlabel('x');
ylabel('log_{10}(error)','Interpreter','tex');
title('2 Ways of Evaluating the Interpolating Polynomial for sin(\pi x) on [0,0.2]','Interpreter','tex');
%ylim([-16 0]);
axis image
xlim ([0,0.2])
axis square

function getSliderPosition(hSlider)
raw=get(hSlider,'Value');
minv=get(hSlider,'Min');
maxv=get(hSlider,'Max');

scamin=1; % scaled minimum
scamax=50; % scaled maximum

sca=(raw-minv)*(scamax-scamin+1)/(maxv-minv)+scamin;
sca=floor(sca);

ud=get(hSlider,'UserData');
ud.nodecount=sca;
set(hSlider,'UserData',ud);

function refreshSliderLabel(hSlider,hLabel)
ud=get(hSlider,'UserData');

set(hLabel,'String',sprintf('%d',ud.nodecount));

% Copyright (c) 2009, Daniel Ellis
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of the Columbia University nor the names 
%       of its contributors may be used to endorse or promote products derived 
%       from this software without specific prior written permission.
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
% POSSIBILITY OF SUCH DAMAGE.

function [P,R,S] = lagrangepoly(X,Y,XX)
%LAGRANGEPOLY  Lagrange interpolation polynomial fitting a set of points
%   [P,R,S] = LAGRANGEPOLY(X,Y)  where X and Y are row vectors
%   defining a set of N points uses Lagrange's method to find 
%   the N-1th order polynomial in X that passes through these 
%   points.  P returns the N coefficients defining the polynomial, 
%   in the same order as used by POLY and POLYVAL (highest order first).
%   Then, polyval(P,X) = Y.  R returns the x-coordinates of the N-1
%   extrema of the resulting polynomial (roots of its derivative),
%   and S returns the y-values  at those extrema.
%
%   YY = LAGRANGEPOLY(X,Y,XX) returns the values of the polynomial
%   sampled at the points specified in XX -- the same as
%   YY = POLYVAL(LAGRANGEPOLY(X,Y)). 
%
%   Example:
%   To find the 4th-degree polynomial that oscillates between 
%   1 and 0 across 5 points around zero, then plot the interpolation
%   on a denser grid inbetween:
%     X = -2:2;  Y = [1 0 1 0 1];
%     P = lagrangepoly(X,Y);
%     xx = -2.5:.01:2.5;
%     plot(xx,polyval(P,xx),X,Y,'or');
%     grid;
%   Or simply:
%     plot(xx,lagrangepoly(X,Y,xx));
%
%   Note: if you are just looking for a smooth curve passing through 
%   a set of points, you can get a better fit with SPLINE, which 
%   fits piecewise polynomials rather than a single polynomial.
%
%   See also: POLY, POLYVAL, SPLINE

% 2006-11-20 Dan Ellis dpwe@ee.columbia.edu
% $Header: $

%  For more info on Lagrange interpolation, see Mathworld: 
%  http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html

% Make sure that X and Y are row vectors
if size(X,1) > 1;  X = X'; end
if size(Y,1) > 1;  Y = Y'; end
if size(X,1) > 1 || size(Y,1) > 1 || size(X,2) ~= size(Y,2)
  error('both inputs must be equal-length vectors')
end

N = length(X);

pvals = zeros(N,N);

% Calculate the polynomial weights for each order
for i = 1:N
  % the polynomial whose roots are all the values of X except this one
  pp = poly(X( (1:N) ~= i));
  % scale so its value is exactly 1 at this X point (and zero
  % at others, of course)
  pvals(i,:) = pp ./ polyval(pp, X(i));
end

% Each row gives the polynomial that is 1 at the corresponding X 
% point and zero everywhere else, so weighting each row by the 
% desired row and summing (in this case the polycoeffs) gives 
% the final polynomial
P = Y*pvals;

if nargin==3
  % output is YY corresponding to input XX
  YY = polyval(P,XX);
  % assign to output
  P = YY;
end

if nargout > 1
  % Extra return arguments are values where dy/dx is zero
  % Solve for x s.t. dy/dx is zero i.e. roots of derivative polynomial
  % derivative of polynomial P scales each power by its power, downshifts
  R = roots( ((N-1):-1:1) .* P(1:(N-1)) );
  if nargout > 2
    % calculate the actual values at the points of zero derivative
    S = polyval(P,R);
  end
end

function [ val ] = newtonpoly( x, f, xx )
%newton_poly evaluates the Newton interpolant
%   INPUT ARGUMENTS
%       x = the interpolation nodes
%       f = the values at these nodes
%   OUTPUT ARGUMENTS
%       val = the value

n=length(x);

dd=zeros(n,n);
dd(:,1)=f;

for l=2:n
    for k=1:(n-l+1)
        dd(k,l)=(dd(k+1,l-1)-dd(k,l-1))/(x(k+l-1)-x(k));
    end
end

p=fliplr(dd(1,:));

% Horner scheme for evaluation
ap = p(1)*ones(size(xx));

for l=1:(n-1)
    ap=p(l+1)+ap.*(xx-x(n-l));
end

val=ap;
