MATLAB Examples

Test C: 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

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

Define the ODE as Latex Equations

Eq = '$$2\,x^2\,\ddot{y} - x\,\dot{y} - 2 y = 0$$, with, $$y(1) = 5$$ and $$\dot{y}(1) = 0$$' ;
Sol = '$$y = x^2 + \frac{4}{\sqrt{x}}$$' ;
%

Set up the Nodes

In this test a set of evenly spaced nodes is generated a-priori. The ode45 solver in MATLAB generates a set of very dense nodes at the start of the solution of the this equation. The new method can solve this problem without the necessity for a particularly dense set of nodes.

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

Compute a Solution using the MATLAB IVP Solver for Comparison

The MATLAB standard ode45 method is used to compute a solution for camprative purposes.

[xML,Y] = ode45(@eulerODE02,[xMin xMax],[5 0]);
%
yML = Y(:,1);
%

Define an Inline Solution 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

ls = 13;
D = dopDiffLocal( x, ls, ls );
%
% Compute the linear differential operator.
%
L = 2*diag(x.^2)*D*D - diag(x)*D - 2*eye(noPts) ;
%

Compute the Constraining Matrix

Define the two initial values as constraints

c1 = zeros( noPts, 1 );
c1(1) = 1; % polace the first constraint at the start of the support
%
% Select the differential for the first node as the second constraint
%
c2 = D(1,:)';
%
% Concatinate the constraining vectors to form the constraining
% matrix.
%
C = [c1, c2];
%
% Define the vector of constraint values.
%
d = [5;0];

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 programmen
% to simplify the understanding of the code in light of the text in
% the paper.
%
M = Nc * pinv( L * Nc );
yp = yc - M * L * yc;
%

Solve the differential Equation

This is the code required for the repeated solution of the differential equation

g = zeros( size( x ));  % This is a homogeneous differential equation
%
% solve the equation
%
y = M * g + yp;
% %
% % Compute the error with respect to the analytical solution
% %
% e = f(x) - y;
% nE = norm( e );
%

Present the results

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', 'NorthWest' );

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.

eML = f(xML) - yML;
NeML = norm( eML );
%
eNew = f(x) - y;
NeNew = norm( eNew );
%
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,-10,0]);
%

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, 'IVP2Sol' );
    figure2eps( fig2, 'IVP2ErrorLog' );
end;
%