Damped harmonic oscillator fitting

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.

Answer 1

[Star Strider]

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.

[Richard Wood]

Hey! Thank you very much for your comment. As you say, the fitting is not excellent. Do you think there it is another approach to deal with this fitting without dealing with a multiexponential decay model? Or do you think that not exist the possibility to fit to this profile using perhaps a more accurate approach? In the case that a multidimensional decay model it is implemented, how could I edit your code to incorporate that possibility?

[Star Strider]

My pleasure.

Fitting a multiexponential decay model requires knowledge of the process that created those data. I doubt that is straightforward to implement.

Feel free to experiment.

[Richard Wood]

Thank you for your comments. A question on which I would like to know your opinion. Imagine that your variable y (I am using your nomenclature) is zero for many values of x, before and after the region in which appreciable values of y other than zero appear that allow you to define a damped harmonic oscillator type profile. If that were the case, would your approach work without problem even if these useless values are present? At least in the accurate calculus of b(3).

[Star Strider]

As always, my pleasure.

I suspect my approach could have problems with the data having a large number of zero values or other constant values at either end of the signal, although I have never experimented with it in that context. My approach to that problem would simply be to remove them (likely detecting them with the diff function), although the exact method would depend on the data.