Code covered by the BSD License  

Highlights from
FiniteDifferenceExplorer

image thumbnail
from FiniteDifferenceExplorer by JM Modisette
Launches an interactive GUI for comparing forward, backward and centered finite difference.

FiniteDifferenceExplorer()
function FiniteDifferenceExplorer()
% GUI for displaying convergence rates of three finite difference routines:
% forward difference, backward difference and centered difference.  A set
% of test functions are provided to illustrate convergence rates and the
% impact of singularities.  
%   Copyright 2013 The MathWorks, Inc.


%% Initialize the parameters
fun{1} = @(x) exp(x);
fun{2} = @(x) sin(2*pi*x);
fun{3} = @(x) sin(4*pi*x);
fun{4} = @(x) exp(x) + 0.05*sin(8*pi*x);
fun{5} = @(x) (x <= 1/3).*sin(pi*x) + (x > 1/3).*0.5.*sin(pi*x);
handles.fun = fun;

dfun{1} = @(x) exp(x);
dfun{2} = @(x) 2*pi*cos(2*pi*x);
dfun{3} = @(x) 4*pi*cos(4*pi*x);
dfun{4} = @(x) exp(x) + 0.2*(2*pi*cos(8*pi*x));
dfun{5} = @(x) (x <= 1/3).*pi*cos(pi*x) + (x > 1/3).*0.5.*pi*cos(pi*x);
handles.dfun = dfun;

handles.FunctionNames = {
    'f(x) = exp(x)'
    'f(x) = sin(2 pi x)'
    'f(x) = sin(4 pi x)'
    'f(x) = exp(x) + 0.05*sin(8 pi x)'
    'f(x) = (x <= 1/3).*sin(pi*x) + (x > 1/3).*0.5.*sin(pi x)'
    };

handles.xdom = [-0.6 0.6];

ylim = [
    1e-5 1e0
    1e-2 1e1
    1e-1 1e2
    1e-2 1e1
    1e-4 1e1
    ];
handles.ylim = ylim;

RunType = {
    'forward'
    'backward'
    'centered'
    };
handles.RunType = RunType;

Methods = {
    'Forward Difference'
    'Backward Difference'
    'Centered Difference'
    };
handles.Methods = Methods;

handles.VisibleFlag = ones(length(RunType),1);

handles.nnint = [1:5];

handles.Color = [
    %    0.00 0.00 1.00
    0.00 0.50 0.00
    1.00 0.00 0.00
    0.00 0.75 0.75
    0.75 0.00 0.75];

handles.LineWidth            = 1.5;
handles.MarkerSizeError      = 20;
handles.MarkerSizeCurrentn   = 40;
handles.MarkerSizeMethodPnts = 10;
handles.FontSize             = 12;

%% build up the GUI
fh = figure('Visible','on','Name','FiniteDifferenceExplorer');%,...
%    'Position',[2360,500,1000,550],...
%    'units','normalized');
%set(fh,'MenuBar','none');
%    'units','normalized');
%    'Position',[2360,500,1000,550],...
set(fh,'MenuBar','none');
pos_orig = get(fh,'Position');
Width = 1000; Height = 550;
set(fh,'Position',[pos_orig(1), pos_orig(2)+pos_orig(4)-Height, Width, Height]);
set(fh,'units','normalized');

% add invisible axis to hold text in case we want latex
handles.textaxis = axes('parent',fh,...
    'units','normalized',...
    'position',[0 0 1 1],...
    'visible','off');

bottom = 0.2;
width  = 0.3;
horiz_space = 0.05;
height = 0.7;
% add axis for plotting solution
handles.PlotAxes = axes('Parent',fh,...
    'Position',[horiz_space, bottom, width, height]);
handles = AddLinesToPlotAxes(handles);

% add axis for plotting error convergence
handles.ErrorAxes = axes('Parent',fh,...
    'Position',[2*horiz_space+width, bottom, width, height]);
handles = AddLinesToErrorAxes(handles);

% add the slider for N
nmin = min(handles.nnint);
nmax = max(handles.nnint);
handles.NControlSlider = uicontrol(fh,'Style','slider',...
    'Max',nmax,'Min',nmin,...
    'Value',1,...
    'SliderStep',1/(nmax-nmin)*[1 1],...
    'FontSize',handles.FontSize,...
    'units','normalized',...
    'Position',[horiz_space, 0.05, width, 0.05],...
    'CallBack',@UpdatePlotsFromNControlSlider);
axes(handles.textaxis)
handles.NControlText = text(horiz_space+width/2,.03,...
    'N = 2^1 = 4',...
    'unit','normalized',...
    'HorizontalAlignment','center',...
    'FontSize',handles.FontSize);

% add function selector
handles.FunctionSelector = uicontrol(fh,'Style','popupmenu',...
    'String',handles.FunctionNames,...
    'Value',1,...
    'FontSize',handles.FontSize,...
    'units','normalized',...
    'Position',[2*horiz_space+2*width+0.025,0.8, .25, .1],...
    'CallBack',@GetErrorData);

% add MethodPanel
handles.MethodPanel = uibuttongroup('Parent',fh,...
    'Title','Finite Differencing Methods',...
    'Position',[2*horiz_space+2*width+0.025,0.4, .25, .4]);

controlleft = 0.1;
controlwidth = 0.9;
controlheight = 0.1;

for i=1:length(handles.Methods)
    handles.MethodBoxes(i) = uicontrol(handles.MethodPanel,'Style','checkbox',...
        'String',handles.Methods{i},...
        'Value',1,...
        'FontSize',handles.FontSize,...
        'units','normalized',...
        'ForegroundColor', handles.Color(i,:),...
        'Position',[controlleft, 1-.2*i, controlwidth, controlheight],...
        'CallBack',@UpdatePlotsFromMethods);
%        'BackgroundColor', handles.Color(i,:),...
end

% add location of x-derivative
handles.xLocationEdit = uicontrol(fh,'Style','edit',...
    'String','0.0',...
    'FontSize',handles.FontSize,...
    'units','normalized',...
    'Position',[2*horiz_space+1.625*width, 0.05, 0.075, 0.05],...
    'CallBack',@GetErrorData);
axes(handles.textaxis)
handles.xLocationText = text(2*horiz_space+1.6*width,.075,'Derivative x-location:',...
    'unit','normalized',...
    'HorizontalAlignment','right',...
    'FontSize',handles.FontSize);

% add checkbox for rate triangle
%handles.RateTriangleFlag = 0;
handles.RateTriangle = uicontrol(fh,'Style','checkbox',...
        'String','Display Convergence Triangles',...
        'Value',0,...
        'FontSize',handles.FontSize,...
        'units','normalized',...
        'Position',[2*horiz_space+2*width+0.025,0.2, .25, .05],...
        'CallBack',@UpdatePlots);
    
    

%% store handles as appdata in fh
setappdata(fh,'handles',handles);

%% make initial plot
GetErrorData(handles.FunctionSelector)

end

function handles = AddLinesToPlotAxes(handles)
PlotAxes = handles.PlotAxes;
% set up the x limit
set(PlotAxes,'XLim',handles.xdom);
xlabel('x');
ylabel('f(x)');

% retrieve the standard colors
Color        = handles.Color;
LW           = handles.LineWidth;
MSmethodpnts = handles.MarkerSizeMethodPnts;


% build necessary line handles
for i=1:length(handles.RunType)
    hMethodLine(i) = line;
    set(hMethodLine(i),'xdata',[],'ydata',[]);
    set(hMethodLine(i),'LineStyle','-')
    set(hMethodLine(i),'Color',Color(i,:))
    set(hMethodLine(i),'LineWidth',LW)
    
    hMethodPnts(i) = line;
    set(hMethodPnts(i),'xdata',[],'ydata',[]);
    set(hMethodPnts(i),'LineStyle','none')
    set(hMethodPnts(i),'Color',Color(i,:))
    set(hMethodPnts(i),'Marker','o')
    set(hMethodPnts(i),'MarkerSize',MSmethodpnts)
    
end

% create function line handle
hFunction  = line;
set(hFunction ,'xdata',[],'ydata',[]);
set(hFunction ,'LineStyle','--');
set(hFunction ,'Color',[0 0 1]);
set(hFunction ,'LineWidth',LW);

% store handles for later use
handles.hMethodLines      = hMethodLine;
handles.hMethodPnts       = hMethodPnts;
handles.hFunction         = hFunction;
end

function handles = AddLinesToErrorAxes(handles)
ErrorAxes = handles.ErrorAxes;
% change yscale to log
set(ErrorAxes,'YScale','log')
% set x-limits
set(ErrorAxes,'XLim',[handles.nnint(1)-1 handles.nnint(end)+1]);
% set axis labels
xlabel('number of refinements');
ylabel('error');

% retrive default colors
Color   = handles.Color;
LW      = handles.LineWidth;
MSerror = handles.MarkerSizeError;
MScurn  = handles.MarkerSizeCurrentn;

% build necessary line handles
for i=1:length(handles.RunType)
    % set up handles for the Error Lines
    hErrorLine(i) = line;
    set(hErrorLine(i),'xdata',[],'ydata',[]);
    set(hErrorLine(i),'Marker','.','LineStyle','-','Visible','off')
    set(hErrorLine(i),'Color',Color(i,:));
    set(hErrorLine(i),'LineWidth',LW);
    set(hErrorLine(i),'MarkerSize',MSerror);
    
    % set up handles for the error convergence rate triangles
    hErrorRateTriangle(i) = line;
    set(hErrorRateTriangle(i),'xdata',[],'ydata',[]);
    set(hErrorRateTriangle(i),'LineStyle','--','Visible','off')
    set(hErrorRateTriangle(i),'Color',Color(i,:));
    set(hErrorRateTriangle(i),'LineWidth',LW);
    
    % set up handles for the text of the error convergence rate
    hErrorRateTriangleText(i) = text;
    set(hErrorRateTriangleText(i),'horizontalalignment','left');
    set(hErrorRateTriangleText(i),'Color',Color(i,:));
    
    % set up handles for the current n data point
    hErrorCurentn(i) = line;
    set(hErrorCurentn(i),'xdata',[],'ydata',[]);
    set(hErrorCurentn(i),'Marker','.','Visible','off')
    set(hErrorCurentn(i),'Color',Color(i,:));
    set(hErrorCurentn(i),'MarkerSize',MScurn);
end

% store the handles for later use
handles.hErrorLines = hErrorLine;
handles.hErrorRateTriangles = hErrorRateTriangle;
handles.hErrorRateTrianglesText = hErrorRateTriangleText;
handles.hErrorCurrentn = hErrorCurentn;

end

function UpdatePlotsFromNControlSlider(hObject, ~)
% get handles
fh = get(hObject,'Parent');
handles = getappdata(fh,'handles');
n = get(handles.NControlSlider,'Value');
% make sure n stays an integer
n = round(n);
set(handles.NControlSlider,'Value',n);
str = sprintf('N = 2^%d = %d',n,2^n);
set(handles.NControlText,'String',str);
UpdatePlots(hObject)
end

function UpdatePlotsFromMethods(hObject,~)
MethodBox = get(hObject,'Parent');
UpdatePlots(MethodBox);
end

function GetErrorData(hObject, ~)
% get handles
fh = get(hObject,'Parent');
handles = getappdata(fh,'handles');
% compute the error
handles = ComputeError(handles);
setappdata(fh,'handles',handles) ;
UpdatePlots(hObject);
end

function UpdatePlots(hObject, ~)
% get handles
fh = get(hObject,'Parent');
handles = getappdata(fh,'handles');

% update the visible flag
for i=1:length(handles.Methods)
    %get(handles.MethodBoxes(i),'Value')
    handles.VisibleFlag(i) = get(handles.MethodBoxes(i),'Value');
end

% update the plots
UpdateErrorPlot(handles)
UpdateFunctionPlot(handles)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% based on new error data update the plot
function UpdateErrorPlot(handles)
onf = {'off','on'};

% get relavent handles
hErrorLine             = handles.hErrorLines;
hErrorRateTriangle     = handles.hErrorRateTriangles;
hErrorRateTriangleText = handles.hErrorRateTrianglesText;
hErrorCurentn          = handles.hErrorCurrentn;
n                      = get(handles.NControlSlider,'Value');
RTF                    = get(handles.RateTriangle,'Value');
%handles.RateTriangleFlag;
nnint                  = handles.nnint;

% loop over the methods
for t = 1:length(handles.RunType)
    % grab the error for this method
    error = handles.Error(t,:);
    % grab the visible flag
    VF = handles.VisibleFlag(t);

    
    % add the error line data
    set(hErrorLine(t),'xdata',nnint,'ydata',error);
    set(hErrorLine(t),'Visible',onf{VF+1})
    
    % add dot for current number of elements
    set(hErrorCurentn(t),'xdata',nnint(n),'ydata',error(n));
    set(hErrorCurentn(t),'Visible',onf{VF+1})
    
    % add the triangles if they are needed
    xx = nnint(end-1:end);
    yy = error(end-1:end);
    xtri = [xx(1) xx(1) xx(2) xx(1)];
    ytri = 0.9*[yy(1) yy(2) yy(2) yy(1)];
    set(hErrorRateTriangle(t),'xdata',xtri,'ydata',ytri);
    set(hErrorRateTriangle(t),'Visible',onf{(VF&&RTF)+1})
    
    % add the text label for the triangles if needed
    rate = log(yy(2)/yy(1))/log((2^xx(1))/(2^xx(2)));
    set(hErrorRateTriangleText(t),'Position',[1/0.98*xx(1),0.7*sqrt(yy(2)*yy(1))]);
    set(hErrorRateTriangleText(t),'String',sprintf('%.2f',rate));
    set(hErrorRateTriangleText(t),'Visible',onf{(VF&&RTF)+1})

end

set(handles.ErrorAxes,'YLim',handles.ylim(get(handles.FunctionSelector,'Value'),:));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% based on request update the function plot
function UpdateFunctionPlot(handles)
onf = {'off','on'};

% get relavent handles
hMethodLine = handles.hMethodLines;
hMethodPnts = handles.hMethodPnts;
hFunction   = handles.hFunction;
n           = get(handles.NControlSlider,'Value');
funindex    = get(handles.FunctionSelector,'Value');
fun         = handles.fun;
xdom        = handles.xdom;
xloc        = str2double(get(handles.xLocationEdit,'String'));

% update exact function data
xref = linspace(xdom(1),xdom(2),1000);
yref = fun{funindex}(xref);
set(hFunction ,'xdata',xref,'ydata',yref);

% loop over the methods and update the method lines and points
for t = 1:length(handles.RunType)
    type = handles.RunType{t};
    y0 = fun{funindex}(xloc);
    h = (1/2)^n;
    dfdxh = FiniteDifferenceDerivative(fun{funindex},xloc,h,type);
    ydom = GetDerivativeLineYDOM(xdom,xloc,y0,dfdxh);

    % fill in the method line data
    set(hMethodLine(t),'xdata',xdom,'ydata',ydom);
    set(hMethodLine(t),'Visible',onf{handles.VisibleFlag(t)+1})
    
    % get the finite difference points
    xp = GetxPnts(xloc,h,type);
    yp = fun{funindex}(xp);
    set(hMethodPnts(t),'xdata',xp,'ydata',yp);
    set(hMethodPnts(t),'Visible',onf{handles.VisibleFlag(t)+1})
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the error associated with current funindex for all methods
function handles = ComputeError(handles)
% grab the location we want to evaluate the derivative about
xloc     = str2double(get(handles.xLocationEdit,'String'));
funindex = get(handles.FunctionSelector,'Value');
fun      = handles.fun{funindex};
dfun     = handles.dfun{funindex};

% loop over the Derivative Methods
for t = 1:length(handles.RunType)
    type = handles.RunType{t};
    
    % loop over our refinement levels and compute the error
    for i = 1:length(handles.nnint)
        n = handles.nnint(i);
        h = (1/2)^n;
        dfdx = dfun(xloc);
        dfdxh = FiniteDifferenceDerivative(fun,xloc,h,type);
        error(t,i) = abs(dfdx-dfdxh);
    end
end
handles.Error   = error;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the finite diffference derivative
function dfdx = FiniteDifferenceDerivative(fun,x,h,type)
    switch lower(type)
        case 'forward'
            dfdx = 1/h*(fun(x+h)-fun(x));
        case 'backward'
            dfdx = 1/h*(fun(x)-fun(x-h));
        case 'centered'
            dfdx = 1/(2*h)*(fun(x+h)-fun(x-h));
        otherwise
            error('unknown differentiation type');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the finite diffference derivative
function xp = GetxPnts(x,h,type)

    switch lower(type)
        case 'forward'
            xp = [x,x+h];
        case 'backward'
            xp = [x-h,x];
        case 'centered'
            xp = [x-h,x,x+h];
        otherwise
            error('unknown differentiation type');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% based on request update the function plot
function ydom = GetDerivativeLineYDOM(xdom,x0,y0,dfdx)
b = y0-dfdx*x0;
ydom = dfdx.*xdom + b;
end


Contact us