Damped harmonic oscillator fitting

9 次查看（过去 30 天）
Roderick
Roderick 2020-3-9
评论： Star Strider 2020-3-11
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.

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2020-3-10
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 个评论
显示 2更早的评论
Roderick
Roderick 2020-3-11
编辑：Roderick 2020-3-11
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
Star Strider 2020-3-11
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Interpolation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by