Code covered by the BSD License  

Highlights from
FiniteDifferenceExplorer

image thumbnail

FiniteDifferenceExplorer

by

 

24 Jul 2013 (Updated )

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