MATLAB Examples

Test E: Generate Data for a Reconstruction Problem

This script is supplemental material to the paper:

"An Algebraic Framework for the Real-Time Solution of Inverse Problems on Embedded Systems"

by: Matthew Harker, Christoph Gugg, Paul O'Leary and Gerhard Rath

This script require our Polynomial Toolbox available at:

http://www.mathworks.com/matlabcentral/fileexchange/41658

and the function to convert figures to eps:

http://www.mathworks.com/matlabcentral/fileexchange/42388

Contents

Defining the Problem

This function generates simulated data for the reconstruction of a structure being monitored by an array of inclinometers. A set of constraints are placed on the problem and a polynomial is generated which fulfils all these constraints. The derivative of the polynomial is explicitly known. Consequently, derivatives can be computed and the constraints defined. The task of the resconstruction is to reconstruct the curve.

Given the reconstruction and the function from which the derivatives were analytically derived enable the comparison of the reconstruction with the analytical solution.

This is a three point boundary value problem, however, with four constraints. The conditions it includes both Dirchlet and Neumann constraints, both homogeneous and non-homegeneous.

Prepare to execute code.

close all
clear all;
%
%
% Set some defaults
%
fontSize = 12;
set(0,'DefaultFigureColor',[1,1,1]);
set(0,'DefaultaxesFontName','Times');
set(0,'DefaultaxesFontSize',fontSize);
set(0,'DefaulttextFontName','Times');
set(0,'DefaulttextFontSize',fontSize);
set(0,'DefaultfigurePaperType','A4');
set(0,'DefaultTextInterpreter', 'latex');%%

Setup the Support

define a problem with two derivative constraints on on each end and an intermediate zero constraint.

minX = 0;
maxX = 1;
%
% This corresponds to a chain of 21 inclinometers
%
nrPts = 21;
x = linspace( minX, maxX, nrPts)';
at = 16;
%
% The final polynomial is of degree 8.
%
degree = 8;

Define the Constraints

The following four constraints are defined

$$ y(n_{16} = 0$$
$$ y(1) = -0.1$$
$$ D y(0) = 1 $$
$$ D y81) = 0 $$

t1 = [1,0,1];
t2 = [0,x(at),0];
t3 = [1,1,0];
t4 = [0,1,-0.1];
%
T = [t1;t2;t3;t4];

Generate the Particular and Homogeneous Solutions

Use the polyBox tools to generate the Vandermonde matrix corresponding to the constraints Vh and the particular solution yh.

[Vh, yp, pp, L, Sc] = polyconstrained( x, degree, T );
%

Generate the Test Function

[n,m] = size( Vh );
%
% There are the coefficients for the "free" porion of the polynomial.
% It is of degree 4. When the 4-constraints are added there will be a
% degree 8 polynomial.
%
cs = [1.1; 0.4; 0.5; -1.2; -0.3];
%
% Map the unconstrained problem to a constrained problem. The
% coefficients csv are for a polynomial which fulfills all the
% constraints.
%
csv = L * cs + pp;
%
% Compute the function y(x) and the analytical derivatives.
%
V = polyvander( x, degree );
M = polydiffcfs( degree );
%
y = V * csv;
%
% It is important to note that these are analytical derivatives and
% not numerical approximations to derivatives.
%
dy = V * M * csv;
%
% Save the results for use diring reconstruction
%
save BVPFitData x y dy T degree cs csv;
%

Present the Test Function

Output the polynomials as strings

disp('Unconstrained polynomial');
strCS = poly2str( cs, 'x')
disp('Constrained polynomial');
strCSV = poly2str( csv, 'x')
disp('Gradient of the Constrained polynomial');
strDCSV = poly2str( M * csv, 'x')
%
FigureSize=[1 1 10 6];
set(0,'DefaultFigureUnits','centimeters');
set(0,'DefaultFigurePosition',FigureSize);
set(0,'DefaultFigurePaperUnits','centimeters');
set(0,'DefaultFigurePaperPosition',FigureSize);
MyAxesPosition=[0.13 0.1 0.86 0.86];
set(0,'DefaultaxesPosition',MyAxesPosition);
%
fig1 = figure;
%
subplot(2,1,1);
plot( x, y, 'k');
grid on;
xlabel('$$x$$');
ylabel('$$y(x)$$');
axis([0, 1, -0.2, 0.3]);
hold on;
%
% Plot the constraints
%
plot( x(at), y(at), 'ko', 'MarkerFaceColor', 'k');
plot( x(end), y(end), 'ko', 'MarkerFaceColor', 'k');
%
subplot(2,1,2);
plot( x, dy, 'k');
grid on;
xlabel('$$x$$');
ylabel('$$D y(x)$$');
hold on;
plot( x(1), dy(1), 'ko', 'MarkerFaceColor', 'k');
plot( x(end), dy(end), 'ko', 'MarkerFaceColor', 'k');
axis([0,1,-1.5,1.5]);
%
% Save the Figure as EPS if the figure2eps function has been installed
%
h = functions( @figure2eps );
if ~ isempty( h.file )
    figure2eps( fig1, 'BVPdata');
end;
Unconstrained polynomial
strCS =
   1.1 x^4 + 0.4 x^3 + 0.5 x^2 - 1.2 x - 0.3
Constrained polynomial
strCSV =
   -0.46985 x^8 + 0.41127 x^7 + 0.34891 x^6 + 0.03827 x^5 + 1.0323 x^4
   - 1.5886 x^3 - 0.88426 x^2 + 1 x + 0.011895
Gradient of the Constrained polynomial
strDCSV =
   -3.7588 x^7 + 2.8789 x^6 + 2.0934 x^5 + 0.19135 x^4 + 4.1292 x^3
   - 4.7657 x^2 - 1.7685 x + 1