MATLAB Examples

Test D: Demonstrate the New Methor of Solving ODEs with an IVP

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 Discrete orthogonal polynomial toolbox, it is available at:

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

In this script we investigate the effect of the support length on the residual of the solution.

The example Equation being Considered

This is a linear differential equation with variable coefficients.

given

and

The analytical solution to this problem is:

Prepare to Execure Code

%\section{Clear up the Workspace} 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');%% % 

Set up the Nodes

Define the number of points (nodes) used in the solution

noPts = 69 ; % The same number of nodes are used as with the ode45 % % Define the interval and use linearly spaced nodes. % xMin = 1 ; xMax = 10 ; x = linspace( xMin, xMax, noPts )'; % 

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('x.^2 + (4./(sqrt(x)))') ; % ya = f(x); % 

Compute the Discrete Linear Differential Operator

Setup the differentiating matrix with support length ls = 13

lss = 3:2:25; for k=1:length(lss) ls = lss(k); D = dopDiffLocal( x, ls, ls ); L = 2*diag(x.^2)*D*D - diag(x)*D - 2*eye(noPts) ; c1 = zeros( noPts, 1 ); c1(1) = 1; % polace the first constraint at the start of the support c2 = D(1,:)'; C = [c1, c2]; d = [5;0]; yc = pinv( C' ) * d; Nc = null( C' ); M = Nc * pinv( L * Nc ); yp = yc - M * L * yc; g = zeros( size( x )); % This is a homogeneous differential equation y = M * g + yp; e = f(x) - y; nE(k) = norm( e ) / norm( ya ); end; % 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.2 0.83 0.77]; set(0,'DefaultaxesPosition',MyAxesPosition); % fig1 = figure; plot( lss, log10(nE), 'k' ); xlabel('Support length $$l_s$$'); ylabel('$$log_{10} (\epsilon)$$'); grid on; % 
Warning: With a support length greater than 13 there may be problems with the Runge phenomena. Warning: With a support length greater than 13 there may be problems with the Runge phenomena. Warning: With a support length greater than 13 there may be problems with the Runge phenomena. Warning: With a support length greater than 13 there may be problems with the Runge phenomena. Warning: With a support length greater than 13 there may be problems with the Runge phenomena. Warning: With a support length greater than 13 there may be problems with the Runge phenomena. 

Save the figures as .eps files for possible publication.

h = @figure2eps; hInfo = functions( h ); if isempty( hInfo.file ) disp( 'You will need to download the function figure2eps if you'); disp( 'wish to generate the eps files for the figures' ); disp( 'The function is available at:'); disp( 'http://www.mathworks.at/matlabcentral/fileexchange/42388'); else figure2eps( fig1, 'IVP2ErrorLs'); end; %