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.

Contents

The example Equation being Considered

This is a linear differential equation with variable coefficients.

$$2\,x^2\,\ddot{y} - x\,\dot{y} - 2 y = 0$$

given

$$y(1) = 5$$

and

$$\dot{y}(1) = 0$$

The analytical solution to this problem is:

$$ y = x^2 + \frac{4}{\sqrt{x}}$$

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;
%