Driving force of driven damped harmonic oscillator missing

조회 수: 7 (최근 30일)
Jerry Yeung
Jerry Yeung 2022년 11월 23일
편집: Sam Chak 2022년 11월 24일
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
Sam Chak 2022년 11월 24일
편집: Sam Chak 2022년 11월 24일
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

추가 답변 (1개)

VBBV
VBBV 2022년 11월 23일
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
Jerry Yeung
Jerry Yeung 2022년 11월 23일
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 2022년 11월 24일
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

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by