Damped Oscillation Equation Fitting

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.

[Star Strider]

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]

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.

Answer 1

[Sam Chak]

Hi @Haardhik

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')

[Haardhik]

Thank you so much for your time and answering and i apologize for the late response,

However, the Omega and c in the curve fitting and the calculation differ, I don't understand why that would be the case.

[Sam Chak]

Hi @Haardhik

The omega is the initial estimation of c or ω based on counting the number zero-crossing events. The value of omega is used in the StartPoint of the fitoptions() function. The value c is the one you should take. However, my proposed model 'a*exp(-b*x)*cos(c*x)' for the position assumes that the initial velocity is zero. If it is not, then follow @Alex Sha's approach.

[Haardhik]

The Initiial Velocity is not zero in this case. However, Thank you so much for the help, it helped me understand MATLAB quite a bit.

[Sam Chak]

Then show the initial velocity, .

[Haardhik]

It's initial velocity according to Vernier Pro Motion Sensor is -0.939057777778 m/s. There is a whole table of data for the same, would that prove useful in anyway. I assumed it would be irrelevant.

[Sam Chak]

Hi @Haardhik,

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]

@Sam Chak

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.

Answer 2

[Alex Sha]

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

[Haardhik]

Sorry for the late response, I only get limited net access

However, Attempting your method returns an error stating

Error using sum
Dimension argument must be a positive integer scalar, a vector of unique positive
integers, or 'all'.
 

[Haardhik]

Sorry for not mentioning earlier, Thank you so much for your time

[Haardhik]

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]

@Alex Sha

If you could provide me the base code you used for this, it would be very helpful