Driving force of driven damped harmonic oscillator missing

Hello. Below is my code for a simple driven damped spring oscillator. For some reason the response is always that of a damped oscillator, regardless of the amplitude of the driving force. I've tried running the code with a much smaller oscillator frequency (and damping coefficient) and it works as expected, but with this set of parameters the effects of the driving force cannot be seen. My code borrows from the template code in the ode45 page. Any help is appreciated.

c0=299792458; %speed of light

lase_lambda = 1064e-9; %wavelength

lase_freq = c0/lase_lambda;

w0 = 2*pi*lase_freq; %Oscillation frequency

tend = 1e-13; %end time

ft = linspace(0,tend,1000); %time vector

f = 5*cos(ft*w0); %Driving force

options=odeset('RelTol',1e-6,'AbsTol',1e-6);

p0 = [1 0]; %initial position and velocity

[t,p] = ode45(@(t,y) SpringFunction(t,y,ft,f), [0 tend], p0, options);

plot(t,p(:,1)) %plot the first column (x) vs. time

function dydt = SpringFunction(t,y,ft,f)

f = interp1(ft,f,t);

c0=299792458; %speed of light

lase_lambda = 1064e-9; %wavelength

lase_freq = c0/lase_lambda;

w0 = 2*pi*lase_freq; %Oscillation frequency

tau2 = 500e-6;

c = 1/tau2*1e10; %Damping coefficient

dydt = zeros(size(y));

dydt(1) = y(2);

dydt(2) = f-c*y(2)-(w0^2)*y(1);

end

Answer 1

[Sam Chak]

Hi @Jerry Yeung

Basically, there is nothing wrong with your original code. The reason that the spring response is not affected by the Driving Force is because the magnitude of the Driving Force (5 N) to too small (almost negligible) relative to the Spring's initial Restoring Force , when is a very large value and.

So, either you select a smaller initial value for where 5 N is relatively significant, or you design a higher magnitude for the driving force. Two cases are included in the code. The Driving Force in Case #2 is state-dependent so that the initial values can be arbitrarily selected.

On practical implementation, you need to check if the initial position is realistically selected at 1 meter. If it is true, then you also need to check if you can find a machine that delivers that amount of large driving force.

tend    = 3e-13;                    % end time
rtol    = 3e-14;                    % relative tolerance
atol    = 1e-28;                    % absolute tolerance
options = odeset('RelTol', rtol, 'AbsTol', atol); 
% p0    = [1.6e-30 0];              % initial position and velocity     CASE #1
p0      = [1 0];                    % initial position and velocity     CASE #2
[t, p]  = ode45(@(t, y) SpringFunction(t, y), [0 tend], p0, options);
plot(t, p(:, 1)),                   % plot the first column (x) vs. time
grid on, xlabel('t')

function dydt = SpringFunction(t, y)
    c0          = 299792458;        % speed of light
    lase_lambda = 1064e-9;          % wavelength
    lase_freq   = c0/lase_lambda; 
    w0          = 2*pi*lase_freq;   % Oscillation frequency
    tau2        = 500e-6; 
    c           = 1/tau2*1e10;      % Damping coefficient
    % f         = 5*cos(t*w0);      % Driving force                 CASE #1
    wc          = 6e13;             % design parameter for Case #2
    f           = (c - 2*wc)*y(2) + (w0^2 - wc^2)*y(1);           % CASE #2
    dydt    = zeros(size(y));
    dydt(1) = y(2);
    dydt(2) = f - c*y(2) - (w0^2)*y(1);
end

Answer 2

[VBBV]

c0=299792458;   %speed of light
lase_lambda = 1064e-9; %wavelength
lase_freq = c0/lase_lambda; 
w0 = 2*pi*lase_freq;         %Oscillation frequency
tend = 1e-13;        %end time
ft = linspace(0,tend,1000);        %time vector
f = 10*cos(ft*w0)    %Driving force

f = 1×1000

10.0000 9.8434 9.3785 8.6198 7.5911 6.3247 4.8601 3.2434 1.5250 -0.2411 -1.9997 -3.6956 -5.2758 -6.6907 -7.8961 -8.8541 -9.5349 -9.9169 -9.9884 -9.7470 -9.2003 -8.3654 -7.2685 -5.9439 -4.4332 -2.7836 -1.0468 0.7228 2.4697 4.1393

options=odeset('RelTol',1e-6,'AbsTol',1e-6); 
p0 = [1 0]; %initial position and velocity
[t,p] = ode45(@(t,y) SpringFunction(t,y,ft,f), [0 tend], p0, options);
plot(t,p(:,1)) %plot the first column (x) vs. time

function dydt = SpringFunction(t,y,ft,f)
f = interp1(ft,f,t);
c0=299792458;   %speed of light
lase_lambda = 1064e-9; %wavelength
lase_freq = c0/lase_lambda;
w0 = 2*pi*lase_freq;    %Oscillation frequency
tau2 = 500e-6; 
c = 1/(tau2*1e10);    %Damping coefficient   check this value too hgh 
dydt = zeros(size(y));
dydt(1) = y(2);
dydt(2) = f-c*y(2)-(w0^2)*y(1);
end

[VBBV]

Your damping coefficient seems to be too unreasnobaly high value, May be you mean this

c = 1/(tau2*1e10); 

in place of

c = 1/tau2*1e10; 

[Jerry Yeung]

The system is still underdamped since my oscillation frequency is 1.7703e15 Hz, while my damping coefficient is 2e13.

[VBBV]

Do you mean the values for damping coefficient (2e13) and oscillator frequency(1.77e15) produce underdamped system ? or overdamped ? are these values fixed for your system ?

c0=299792458;   %speed of light
lase_lambda = 1064e-9; %wavelength
lase_freq = c0/lase_lambda; 
w0 = 2*pi*lase_freq    %Oscillation frequency

w0 = 1.7703e+15

tau2 = 500e-6; 
c = 1/(tau2*1e10)    %

c = 2.0000e-07

[Jerry Yeung]

The two values produce an underdamped system since in my code, c/w0 < 1. The oscillator frequency is fixed but the damping coefficient is not. I chose this value in hopes of observing how the driving force affects the system.

[VBBV]

do you mean the effec of driving force effect on vibration amplitude ?

c0=299792458;   %speed of light
lase_lambda = 1064e-9; %wavelength
lase_freq = c0/lase_lambda; 
w0 = 2*pi*lase_freq;         %Oscillation frequency
tend = 1e-13;        %end time
ft = 0.5;  

options=odeset('RelTol',1e-6,'AbsTol',1e-6); 
p0 = [0 0]; %initial position and velocity
F = [1 100 1000]  % different force amplitudes

F = 1×3

1 100 1000

for k = 1:length(F)
     f = F(k)*cos(ft*w0);    %Driving force
     [t,p] = ode45(@(t,y) SpringFunction(t,y,ft,f), [0 tend], p0, options);
     plot(t,p(:,1)) %plot the first column (x) vs. time  
     hold on
end
legend('F = 1','F = 100','F = 1000')

function dydt = SpringFunction(t,y,ft,f)
% f = interp1(ft,f,t);
c0=299792458;   %speed of light
lase_lambda = 1064e-9; %wavelength
lase_freq = c0/lase_lambda;
w0 = 2*pi*lase_freq;    %Oscillation frequency
tau2 = 500e-6; 
c = 1/(tau2*1e5);    %Damping coefficient   check this value too hgh 
dydt = zeros(size(y));
dydt(1) = y(2);
dydt(2) = f-c*y(2)-(w0^2)*y(1);
end