MATLAB Examples

# Test A Speed: Solution of an IVP on Arbitrary Nodes Speed Test

This script is supplemental material to the paper:

"Model Based Design for the Real-Time Solution of Inverse Problems in Cyber Physical Systems"

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

This script requires our discrete orthogonal polynomial toolbox, it is available at:

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

and the figure2eps function to generate the eps figures:

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

## The example being Considered

This is a linear non-homogeneous third order differential equation.

given

and

The analytical solution to this problem is:

## Prepare to Execure 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');%% % 

## Define the number of Solutions to be Computed

noSols = 10000; 

## Define the range of the solution

Define the interval and use linearly spaced nodes.

xMin = 0 ; xMax = 8 ; % 

## Compute a Solution using the MATLAB IVP Solver

The MATLAB standard ode45 method is used to compute a solution for camprative purposes. This routine returns the vector of abscissae xML and the solution vector in y = Y(:,1). xML is then used as the vector of positions for which the new method should compute the solution. This demonstrates the capability of the system to compute solutions for arbitrary given nodes.

t = cputime; for k = 1 : noSols [xML,Y] = ode45(@eulerODE03,[xMin xMax],[3 -3 -47]); end; ode45Time = cputime - t; % yML = Y(:,1); % % Use the x coordinates from the MATLAB solution so the the comparison % of the methods is done for exactly the same node placements. % noPts = length( xML ); x = xML; 

## Define an Inline Folution for the Analytical Solution

The analytical solution is computed im comparison to the MATLAB solution and for the new solution presented in the paper.

f = inline('(3 - 25 * x.^2).* exp(-x) + 5 * (x.^3).*exp(-x)') ; % ya = f(x); % 

## Start Preparatory Computations

------------------------------------------------------------ This is where the preparatory computations start. These computations can be peformed advance and need not be computed at run-time.

## Compute the Discrete Linear Differential Operator

Setup the differentiating matrix with support length ls = 13

ls = 11; D = dopDiffLocal( x, ls, ls ); D2 = D * D; D3 = D * D2; % % Compute the linear differential operator. % L = D3 + 3*D2 + 3*D + eye( noPts ); % 

## Compute the Constraining Matrix

The details behind these computations can be found in:

@article{o2012framework, title={A Framework for the Evaluation of Inclinometer Data in the Measurement of Structures}, author={O'Leary, Paul and Harker, Matthew}, journal={Instrumentation and Measurement, IEEE Transactions on}, volume={61}, number={5}, pages={1237--1251}, year={2012}, publisher={IEEE} }

Define three vectors, one for each of the constraints: The value constraint

c1 = zeros( noPts, 1 ); c1(1) = 1; % polace the first constraint at the start of the support % % The constraints on the first and second derivatives % c2 = D(1,:)'; c3 = D2(1,:)'; % % Concatinate the constraining vectors to form the constraining % matrix. % C = [c1, c2, c3]; % % Define the vector of constraint values. % d = [ 3 ; -3; -47 ]; 

## Generate yc and M

yc = pinv( C' ) * d; % % generate an orthonormal set of basis functions for the null space of % C'. % Nc = null( C' ); % % The following code is not optimized for efficiency, it is programmed % to simplify the understanding of the code. It is a direct one-to-one % implementation of the equations as presented in the paper. % % The matrix M and the vector yh would now be saved and made available % to the embedded computation system. % M = Nc * pinv( L * Nc ); yh = yc - M * L * yc; % 

## Run-Time Computation

------------------------------------------------------------------- The matrix M and the vector yh would now be saved and made available to the embedded computation system. The conde in the next portion is that portion of code required to function in real time.

## Solve the differential Equation

The following line simulates acquiring data from the sensors.

g = 30*exp(-x); % % Solve the differential equation, this is the code required for the % repeated solution of the differential equation % t = cputime; for k = 1 : noSols y = M * g + yh; end; runTime = cputime - t; 

## Characterize and present the results

Compute the errors and the 2-norm of the errors for the Matlab solution and for the new metod

eML = f(x) - yML; NeML = norm( eML ); % eNew = f(x) - y; NeNew = norm( eNew ); % 

## Display the times

disp(['ode45 time = ',num2str(ode45Time)]) ; disp(['preparatory time = ',num2str(prepTime)]) ; disp(['run time = ',num2str(runTime)]) ; 
ode45 time = 31.4219 
Undefined function or variable 'prepTime'. Error in TestASpeed (line 201) disp(['preparatory time = ',num2str(prepTime)]) ;

## Generate the desired plots

The following plot shows the analytical solution, the ode45 solution adn the solution generated using the new method.

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.15 0.86 0.82]; set(0,'DefaultaxesPosition',MyAxesPosition); % fig1 = figure; plot( x, ya, 'k' ) ; hold on plot( x, y, 'ko' ) ; plot( xML, yML, 'kd'); % xlabel('$$x$$'); ylabel('$$y(x)$$'); grid on; range = axis; %plot( range(1:2), [0,0],'k'); %axis([0.5,10.5,0,110]); legend( 'Analytical', 'New', 'ode45', 'Location', 'SouthEast' ); 

## Plot the Error Functions

The following two plots show the residual error for the new method and for the ode45 solver. To ensure camparability the same number of nodes are used for the new solution as are generated by the ode45 solver.

fig2 = figure; plot( x, log10(abs(eNew)), 'k'); hold on; plot( xML, log10(abs(eML)), 'k-.'); xlabel('$$x$$') ylabel('$$\mathrm{log}_{10}(\epsilon)$$'); grid on; % newLeg = ['New $$| \epsilon |_2 = ',num2str(NeNew),'$$']; mlLeg = ['ode45 $$| \epsilon |_2 = ',num2str(NeML),'$$']; legend(newLeg, mlLeg, 'Location', 'SouthEast' ); axis([0,8,-15,0]);