MATLAB: Damped harmonic oscillator fitting

fitting damped harmonic oscillator

Hello everyone,
I am trying to fit my data to a damped harmonic oscillator with functional form: mx=A*cos(omega*time-phi)*exp(-gamma*time)-B. I have attached how it looks my data once it is plotted.
Now, to perform the fitting, I have followed this very complete approach given by Star Strider: https://es.mathworks.com/matlabcentral/answers/484489-damped-harmonic-motion-curve-fit
My code following its scheme it is the following:
mx_breather=detrend(mx_inelastic_localized_breather); % Remove Linear Trend
max_mx_breather=max(mx_breather);
min_mx_breather=min(mx_breather);
difference_mx_breather=(max_mx_breather-min_mx_breather);                                                                   % Range of ‘y’
z_mx_breather=mx_breather-max_mx_breather+(difference_mx_breather/2);
zci_mx_breather=@(v)find(v(:).*circshift(v(:),[-1 0])<= 0); % Returns Approximate Zero-Crossing Indices Of Argument Vector
zt_mx_breather=physical_time_inelastic_breather(zci_mx_breather(mx_breather));
per_mx_breather=2*mean(diff(zt_mx_breather));                                                         % Estimate period
ym_mx_breather=mean(mx_breather);                                                                   % Estimate offset
fit_breather=@(b,x) b(1).*exp(-b(2).*x).*(cos(x.*b(3)-b(4)))+b(5);   % Objective Function to fit
fcn_breather=@(b) norm(fit_breather(b,physical_time_inelastic_breather)-mx_breather);                                                  % Least-Squares cost function
[coefficients_breather,nmrs_breather]=fminsearch(fcn_breather,[difference_mx_breather...
    ;-10;per_mx_breather;-1;ym_mx_breather]) % Minimise Least-Squares
xp_breather=linspace(min(physical_time_inelastic_breather),max(physical_time_inelastic_breather),500);
In principle, it does not give an error when it is compiled. It gives the following values for the fitting:
b(1)=2.69632459149622e-06
b(2)=-10.7656220403261
b(3)=2.02496387156646e-12
b(4)=-1.73792117849851
b(5)=-7.06802818282082e-18
However, when I plot it accompanying my simulated data (exposed in the previous plot), there is no trace of it. My environment for plotting it is the usual:
u7=figure(7)
plot(physical_time_inelastic_breather.*(10^(12)),mx_inelastic_localized_breather,...
    '-b','LineWidth',2)
hold on
plot(xp_breather,fit_breather(coefficients_breather,xp_breather),'--r','LineWidth',2)
hold off
xlabel('Time, $t \, \, \left( \mathrm{ps} \right)$','FontSize',14,'interpreter','latex')
ylabel('$m_x$','FontSize',14,'interpreter','latex')
xlim([physical_time_inelastic(1573).*(10^(12)) physical_time_inelastic(1707).*(10^(12))]);
t7=title('Spatially localized breather via annihilation','FontSize',14,'interpreter','latex') 
set(t7,'interpreter','latex','FontSize',14)
l7=legend('Simulated data','Fitted curve','FontSize',14,'Location','best');
set(l7,'interpreter','latex','FontSize',14)
% cs7='\begin{tabular}{lll}  $\mathrm{A} \cos \left( \omega t - \phi \right) e^{-\gamma t}+\mathrm{b}$ \\ $\mathrm{Nontopological \, \, solitons} \in \left[ \tilde{v}_-, \tilde{v}_+ \right]$ \\ $\mathrm{Nonlinear \, \, spin \, \, waves} \in \left[ \left. \tilde{v}_+, + \infty \right) \right.$ \end{tabular}';
% ha7=annotation('textbox',[0.25 0.75 0.43 0.15], 'Interpreter', 'latex','FontSize',14,'BackgroundColor','w','FaceAlpha',1);
% set(ha7,'String',cs7);
set(gca,'TickLabelInterpreter','latex','FontSize',14)
set(u7,'Units','Inches');
posu7=get(u7,'Position');
set(u7,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[posu7(3),posu7(4)])
saveas(gcf,fullfile(outputdir,['mx_Magnetization_Component_Evolution_Breather_Non_Zero_Region_Naked_Eye.pdf']))
Any idea?
I have attached my data in case someone wants to suggest other way to do it. The first column represents the time, or the variable physical_time_inelastic_breather in my code ( it is in seconds, so if you want to have picoseconds as me in the plot, you have to multiply each element of it by E12), and the second column is mx, that it is, the variable mx_inelastic_localized_breather in my code.

Best Answer

The exponential ddecay is likely multi-exponential, since a single expoinential provides a decent although inexact fit: 
D = load('Breather_Data.m');
x = D(:,1);
y = D(:,2);
y = detrend(y);                                                                 % Remove Linear Trend
yu = max(y);
yl = min(y);
yr = (yu-yl);                                                                   % Range of ‘y’
yz = y-yu+(yr/2);
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                            % Returns Approximate Zero-Crossing Indices Of Argument Vector
zt = x(zci(y));
per = 2*mean(diff(zt));                                                         % Estimate period
ym = mean(y);                                                                   % Estimate offset
fit = @(b,x)  b(1) .* exp(b(2).*x) .* (sin(2*pi*x.*b(3) + b(4))) + b(5);        % Objective Function to fit
fcn = @(b) norm(fit(b,x) - y);                                                  % Least-Squares cost function
[s,nmrs] = fminsearch(fcn, [yr; -1E+11;  1/per;  -1;  ym])                      % Minimise Least-Squares
xp = linspace(min(x),max(x), 500);
figure
plot(x,y,'b', 'LineWidth',1.5)
hold on
plot(xp,fit(s,xp), '--r')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
legend('Original Data',  'Fitted Curve')
text(min(xlim)+0.05*diff(xlim), min(ylim)+0.8*diff(ylim), sprintf('$y = %.3E\\cdot e^{%.3E\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.3E + %.3f) %+.3f$', s), 'Interpreter','latex')
This is slightly revised from the previous code, since your data required some special considerations because of the magnitude of the independent variable.
Related Question