Driving force of driven damped harmonic oscillator missing

7 visualizaciones (últimos 30 días)
Jerry Yeung
Jerry Yeung el 23 de Nov. de 2022
Editada: Sam Chak el 24 de Nov. de 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

Respuesta aceptada

Sam Chak
Sam Chak el 24 de Nov. de 2022
Editada: Sam Chak el 24 de Nov. de 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

Más respuestas (1)

VBBV
VBBV el 23 de Nov. de 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
  5 comentarios
Jerry Yeung
Jerry Yeung el 23 de Nov. de 2022
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
VBBV el 24 de Nov. de 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

Iniciar sesión para comentar.

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by