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 )