# Driving force of driven damped harmonic oscillator missing

7 views (last 30 days)
Jerry Yeung on 23 Nov 2022
Edited: Sam Chak on 24 Nov 2022
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

Sam Chak on 24 Nov 2022
Edited: Sam Chak on 24 Nov 2022
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

### More Answers (1)

VBBV on 23 Nov 2022
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 on 24 Nov 2022
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

### Categories

Find more on Coordinate Transformations and Trajectories in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by