Damped Oscillation Equation Fitting

19 views (last 30 days)
Haardhik
Haardhik on 29 Oct 2023
Commented: Haardhik on 9 Nov 2023
I am attempting to fit Data Tables in Excel to a Sin curve in order to determine its Omega for a physics experiment, where I need the equation to be in the form
x' = Ae^-bt sin(ωt' + φ)
Unfortunately I have no idea about Matlab but when I plot them as Column vectors, it is the same as I would recieve in excel and nowhere close to a sin curve. Any and all help is appreciated, as I have been through the previous questions that are similar and they pick up from a sin curve, unable to find instructions to that step, I am posting this here. Attached is the different trials I wish to do this for.
  2 Comments
Star Strider
Star Strider on 6 Nov 2023
Your data may be undersampled and could show aliasing. If you can, repeat the experiment with a much faster sampling frequency (100 times the current sampling frequency would work), then do the analysis.
Haardhik
Haardhik on 7 Nov 2023
This was the fastest possible sampling frequency on the equipment, it's just a high school physics lab with a basic Vernier Caliper Motion Sensor.

Sign in to comment.

Answers (2)

Sam Chak
Sam Chak on 29 Oct 2023
You can probably estimate the value of omega (angular frequency) by counting the number of zero-crossing events. If you want to fit the damped sinusoidal equation to the data, you can try this:
Tab = readtable('No Mass (final).xlsx')
Warning: Column headers from the file were modified to make them valid MATLAB identifiers before creating variable names for the table. The original column headers are saved in the VariableDescriptions property.
Set 'VariableNamingRule' to 'preserve' to use the original column headers as table variable names.
Tab = 112×6 table
Time_s_ Position_m_ Time_s__1 Position_m__1 Time_s__2 Position_m__2 _______ ___________ _________ _____________ _________ _____________ 0.04 -0.0193 0.06 -0.01769 0.06 -0.013135 0.1 0.022134 0.1 0.020726 0.1 0.016225 0.14 -0.022867 0.14 -0.01769 0.14 -0.014782 0.18 0.022683 0.2 0.01661 0.18 0.01403 0.24 -0.017379 0.24 -0.01769 0.24 -0.013135 0.28 0.02186 0.28 0.019902 0.28 0.015677 0.32 -0.021495 0.32 -0.015495 0.32 -0.013135 0.36 0.020488 0.38 0.017158 0.38 0.012658 0.42 -0.017654 0.42 -0.017416 0.42 -0.01341 0.46 0.021311 0.46 0.01853 0.46 0.014853 0.5 -0.020398 0.5 -0.013849 0.5 -0.011763 0.54 0.018293 0.56 0.017433 0.56 0.013207 0.6 -0.017928 0.6 -0.016593 0.6 -0.013135 0.64 0.021037 0.64 0.017158 0.64 0.01403 0.68 -0.017928 0.7 -0.012751 0.68 -0.010117 0.74 0.016372 0.74 0.017433 0.74 0.013481
%% Trial 1
t = Tab.Time_s_(1:end-2); % Time
y = Tab.Position_m_(1:end-2); % Position
[~, count] = zerocrossrate(y, Method="comparison"); % count zero-crossing
T = t(end)/(count/2); % Period of a sine wave
omega = 2*pi/T % Estimate the angular frequency
omega = 68.4867
% Curve-Fitting
fo = fitoptions('Method','NonlinearLeastSquares',...
'Lower',[ 0, 0, omega-5],...
'Upper',[0.03, 1, omega+5],...
'StartPoint',[max(y) 0.5, omega]);
ft = fittype('a*exp(-b*x)*cos(c*x)', 'options',fo);
[curve, gof] = fit(t, y, ft) % check if omega ≈ c
curve =
General model: curve(x) = a*exp(-b*x)*cos(c*x) Coefficients (with 95% confidence bounds): a = 0.02388 (0.02315, 0.0246) b = 0.2507 (0.2357, 0.2657) c = 69.05 (69.03, 69.07)
gof = struct with fields:
sse: 1.7646e-04 rsquare: 0.9907 dfe: 107 adjrsquare: 0.9905 rmse: 0.0013
% Making plots
plot(t, y, '.'), hold on
plot(curve), grid on
xlabel('Time'), ylabel('Position'), title('Trial 1')
legend('Data', 'Fitted Sine')
  7 Comments
Sam Chak
Sam Chak on 6 Nov 2023
The reason I asked for the initial velocity is to verify the values of a, b, c, and d found by Alex. Additionally, your data is undersampled because there are fewer than 2 data points per cycle (period). Please follow @Star Strider's advice to obtain more data.
% Using Alex Sha's finding
a = 0.023540378407752; % amplitude
b = 0.247177601901692; % decay rate constant
c = 69.1369475428933; % omega
d = 1.32644349058437; % phi
% period
T = 2*pi/c
T = 0.0909
% number of points per second
numPts = 110/5 % 110 data points over 5 sec
numPts = 22
% How many point in a period?
T*numPts % less than two points
ans = 1.9994
% Proposed fitting model
f = @(t) a*exp(-b*t).*sin(c*t + d);
t = linspace(0, 5, 19811); % minimum 19810 points over 5 sec to have 360 points per period
plot(t, f(t)), xlabel('t'), title('f(t) = a*exp(-b*t).*sin(c*t + d)')
% Velocity by formula (df/dt)
df = a*c*exp(-b*t).*cos(c*t + d) - a*b*exp(-b*t).*sin(c*t + d);
plot(t, df), xlabel('t'), title('df/dt')
% initial velocity based on df formula
df0 = df(1)
df0 = 0.3881
% initial velocity computed directly from {a, b, c, d}
df0 = a*c*cos(d) - a*b*sin(d)
df0 = 0.3881
% Test if initial velocity is -0.939057777778
True_df0 = -0.939057777778;
ToF = logical(abs(df0 - True_df0) < 1e-4) % True (1) / False (0)
ToF = logical
0
Haardhik
Haardhik on 9 Nov 2023
I realized that this data is sifted only for the peaks, If I were to apply the same to the actual data that sampled every 0.02 seconds, would it be possible?
I also noticed that you used a Cos wave in ur Fittype and hence will attach the original data herewith.

Sign in to comment.


Alex Sha
Alex Sha on 29 Oct 2023
if taking fitting function as: x=a*exp(-b*t)*sin(w*t+phi)
for trial-1
Sum Squared Error (SSE): 9.1948847955911E-5
Root of Mean Square Error (RMSE): 0.000914274913678052
Correlation Coef. (R): 0.999583948620075
R-Square: 0.999168070338902
Parameter Best Estimate
--------- -------------
a 0.023540378407752
b 0.247177601901692
w 69.1369475428933
phi 1.32644349058437
for trial-2
Sum Squared Error (SSE): 0.000115663498814737
Root of Mean Square Error (RMSE): 0.0010162233075687
Correlation Coef. (R): 0.999720103291284
R-Square: 0.999440284924735
Parameter Best Estimate
--------- -------------
a 0.0200937306323941
b 0.236486196043985
w 69.1712756320905
phi 0.896068977752904
for trial-3
Sum Squared Error (SSE): 9.64366701017117E-5
Root of Mean Square Error (RMSE): 0.000936320992461801
Correlation Coef. (R): 0.999614594177523
R-Square: 0.999229336892693
Parameter Best Estimate
--------- -------------
a 0.0158395749665591
b 0.215537219553407
w 69.1871168416771
phi 7.30505431226722
Note: phi is not unique
  4 Comments
Haardhik
Haardhik on 6 Nov 2023
I assume you used the CFTool function. I too attempted the same but in the custom equation part, I always got it wrong, could you please tell me how you reached here. I attempted the fitting fucntion you gave above and still failed.
Haardhik
Haardhik on 9 Nov 2023
If you could provide me the base code you used for this, it would be very helpful

Sign in to comment.

Products


Release

R2023b

Community Treasure Hunt

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

Start Hunting!