MATLAB Answers

Damped harmonic oscillator fitting

31 views (last 30 days)
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.

Accepted Answer

Star Strider
Star Strider on 10 Mar 2020
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.
  4 Comments
Star Strider
Star Strider on 11 Mar 2020
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by