image thumbnail

Dassl Mex file compilation to Matlab 5.3 and 6.5!

by

 

16 Oct 2007 (Updated )

I have compiled a Dassl solver version (MATLAB MEX INTERFACE), in order to use with Matlab (version

run_Dassl_test_problem.m
%
%#############################################################################
%
%                            DASSL TEST PROBLEM 
%
%#############################################################################
%
% Testing using a DAE problem and an ODE problem (compare with ode45)
%
% requires files: dydt.m

% global debug level
%   0 : no output
%   1 : write output values only
%  10 : write all results from m-files
% 100 : confirm exit from each m-file with <RETURN>

%___________________________________________________________________________
%___________________________________________________________________________
%
% Edited by: Giovani Tonel, Taking Master Degree in Chemical Engineering
% Chemical Engineering Department - GIMSCOP/ALSOC
% Federal University of Rio Grande do Sul - UFRGS
% Porto Alegre, RS - Brazil
% Phone: +55 51 3308 4166 (GIMSCOP - Room 2)
% Email address: giotonel@enq.ufrgs.br  
% Website: http://www.enq.ufrgs.br/gimscop
% October 2007; Last revision: 15-October-2007
%___________________________________________________________________________
%___________________________________________________________________________

clear all
clc

disp(' ')
disp(' #############################################################################')
disp(' #############################################################################')
disp(' #########################                      ##############################')
disp(' #########################  DASSL TEST PROBLEM  ##############################')
disp(' #########################                      ##############################')
disp(' #############################################################################')
disp(' #############################################################################')
disp(' ')

global DEBUG;

DEBUG   = 1;


Nm      =   'dydt';         % derivative name needs EXACTLY 4 letters

t0      =   0.0;            % initial value for independent variable
tf      =   1.0;            % final value for independent variable
NEQ     =   3;              % number of equations in system
x0      =   [3 6 5]';       % initial state variables
RTOL    =   1.E-4;          % relative tolerance
ATOL    =   1.E-4;          % absolute tolerance

% RPAR  =   [3 4];          % optional arguments passed to file Nm

% testing for struct entering 
% rpar.a  =   [3 4];
% rpar.b  =   [1 5 3];
%
% Dassl doesn't accept a struct variable!!!
% 
% rpar must be a maximium vector with  100-by-1
rpar    =   [1 5 3 1 5 6];

% dassl_mat_6_5.dll         for 6.5 MATLAB version
% dassl_mat_5_3.dll         for 5.3 MATLAB version 

% Note: also plot_do needs to be changed on dydt.m .
plot_do = 0;

if plot_do
    
    % DAE problem
    %
    %####################### Calling dassl.dll ##############################
    %
    tic
    [IDID,NRS,YM]= dassl_mat_6_5(Nm, t0, tf, NEQ, x0, RTOL, ATOL, rpar);
    tdassl          = toc;
    %
    %
    %########################################################################
    %
else
    % ODE problem
    %
    %####################### Calling dassl.dll ##############################
    %
    tic
    [IDID,NRS,YM]   = dassl_mat_6_5(Nm, t0, tf, NEQ, x0, RTOL, ATOL, rpar);
    tdassl          = toc;
    %
    %
    %########################################################################
    %    
    % 
    %############################### ODE45 ##################################
    %
    opt = odeset('RelTol', RTOL,'AbsTol', ATOL);
    %
    tic
    [T, Y]		= 	ode45(@dxdt, [t0, tf], x0, opt, rpar);
    tode45      =   toc;
    %
    %########################################################################
    
    % Dassl's time faster factor than ode45
    %
    dtf     = tode45/tdassl;
    sprintf('Dassl''s time faster factor than ode45= %0.4f', dtf)
end



if ( IDID > 0 )
    disp(' ')
    disp(sprintf('Dassl: integration completed successfully, %d points',NRS));
    disp(' ')
   
    if ( DEBUG > 0 )
        
        % result vector contains t in row 1, y(k) in row k
        t = YM(1,:);
        y = YM(2:NEQ+1,:);
        for i = 1:size(YM,2)
            buf = sprintf('% .2e ',y(:,i));
            %disp(sprintf('pt %3d: t=% .2e y = %s',i,t(i),buf));
        end
        
        if plot_do
            % DAE problem solution  by Dassl
            %
            plot(t, y(1,:), t, y(2,:), t, y(3,:));
            title('DAE problem solution by Dassl')
            legend('y_1_,_D_a_s_s_l','y_2_,_D_a_s_s_l','y_3_,_D_a_s_s_l')
            text(0.5, 4, sprintf('Tsim,Dassl= %0.2f [s]', tdassl))
            
        else
            % ODE problem solution  by Dassl and ode45
            %
            plot(t, y(1,:), t, y(2,:), t, y(3,:), T, Y(:,1), T, Y(:,2), T, Y(:,3));
            legend('y_1_,_D_a_s_s_l','y_2_,_D_a_s_s_l','y_3_,_D_a_s_s_l',...
                'y_1_,_o_d_e_4_5','y_2_,_o_d_e_4_5','y_3_,_o_d_e_4_5')
            title('ODE problem solution comparison between Dassl and ode45')
            
            text(0.5, 4, sprintf('Tsim,Dassl= %0.4f [s]', tdassl))
            text(0.5, 3, sprintf('Tsim,ode45= %0.4f [s]', tode45))
            text(0.2, 2, sprintf('Dassl''s time faster factor than ode45= %0.1f', dtf)) 
        end
        
    end
    
else
    disp(' ')
    disp(sprintf('Dassl: integration failed: result flag = %d',IDID));
    
end


Contact us