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 Center および File ExchangeOrdinary Differential Equations についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by