fitting critically damped simple harmonic oscillator

by

 

Four-parameter fit to the equation A*(1 + t/T2).*exp(-t/T1) + offset

fit_critical_damping( varargin )
function varargout = fit_critical_damping( varargin )
%[A, T1, T2, offset, Ssq] = fit_critical_damping( t, y, options )
%fitting critically damped oscillations to the equation
%y(t) = A (1 + t/T2) exp( -t/T1 ) + offset
%where A is the amplitude, T1 is a time constant, T2 is the decay time, along with
%the offest. Ssq is the sum of residuals (i.e. fitting error).
%
%options are 'notify' to show progress and 'plot' to give a plot
%you can also supply a guess or series of guesses via a row vector or
%matrix, m, where size( m, 2 ) = 4 (i.e. there are four columns);
%
%if you just run program without arguments it will do an example
%
%examples are
%fit_critical_damping( t, y, 'notify', 'plot')
%[A, T1, T2, offset, Ssq] = fit_critical_damping( t, y, [1 2 3 4] )
%
%Note: this script uses LMFnlsq from MATLAB central file exchange
%sak@wpi.edu 2/8/12

fnc = @(A, T1, T2, offset, t ) A*(1 + t/T2).*exp(-t/T1) + offset;

if nargin == 0
    disp( 'doing sample fitting' );
    t = linspace( 3, 100, 1000 );
    [A, T1, T2, offset] = deal( 1, 1, -.1, 0.1 );
    y = fnc( A, T1, T2, offset, t )+randn(size( t ) )*.01;
    varargin = { t, y, 'plot', 'Xnotify' };
end;
[t, y] = deal( shiftdim( varargin{1} ), shiftdim( varargin{2} ) );
varargin(1:2) = [];

plot_flag = strcmpi( varargin, 'plot' );
varargin( plot_flag ) = [];

if any( strcmpi( varargin, 'notify' ) );
    notify = @(str) fprintf( str );
else
    notify = @(str) '';
end
varargin( strcmpi( varargin, 'notify' ) ) = [];

if any( plot_flag )
    figure( sum( mfilename ) );
    clf;
    plot( t, y, '.r' );
end

if any( cellfun( 'isclass', varargin, 'double' ) )
    guess = varargin{cellfun( 'isclass', varargin, 'double' )};
else
    [T1, T2] = ndgrid( diff( t([1 end]) )*10.^[-3:2], diff( t([1 end]) )*[-10.^[-4:3] 10.^[-4:3]] );
    guess = [repmat( y(1), numel( T1 ), 1), reshape( T1, [], 1), reshape( T2, [], 1), repmat( mean(y(ceil(0.8*end):end)), numel( T1 ), 1 ) ];
end

param = guess(1,:);
J = @(param) [ [ [ones(size(t)), t.*param(1)/param(2)^2].*repmat( (1+t/param(3)), 1, 2), -param(1)/param(3)^2*t].*repmat( exp(-t/param(2) ), 1, 3), ...
    ones(size(t))];
warning off;
param = nan( size( guess ) );
[Ssq, CNT] = deal( nan( size( guess, 1), 1) );
tic;
for i=1:size(guess, 1);
    notify( sprintf( 'trying param A = %g, T1 = %g, T2 = %g, phase = %g, offset = %g ', guess(i,:) ) );
    [param(i,:), Ssq(i), CNT(i)] = LMFnlsq( @(param) fnc( param(1), param(2), param(3), param(4), t ) - y, guess(i,:), 'Jacobian', J, 'MaxIter', 20 );
    notify( sprintf( 'error %g \n', Ssq(i) ) );
end
notify(sprintf( '%.2f sec \n', toc ) );
warning on;
[best, i] = min( Ssq );

%refining solution
[param(i,:), Ssq(i), CNT(i), Res, XY] = LMFnlsq( @(param) fnc( param(1), param(2), param(3), param(4), t ) - y, param(i, :), 'Jacobian', J );
A = param(i,1);
T1 = param(i,2);
T2 = param(i,3);
offset = param(i, 4);
if any( plot_flag)
    hold on;
    plot( t, fnc(A, T1, T2, offset, t ), '-k' );
    title( sprintf( 'A = %g, T1 = %g, T2 = %g, offset = %g ', A, T1, T2, offset ) );
end

varargout = {A, T1, T2, offset, Ssq};
varargout = varargout(1:nargout);

Contact us