Plotting harmonic oscillators for different values of Beta

조회 수: 4 (최근 30일)
Emily McLain
Emily McLain 2022년 9월 12일
답변: Mathieu NOE 2022년 9월 12일
I need to plot harmonic oscillators depending on their beta values. I'm having trouble plotting the data that's in the if/elseif statements. My code is down below.
w0= 4*pi;
x0= 100;
beta= [0, 2, 5, 10, 20,50];
time=linspace(0, 120, 6);
%%Undamped Oscillator
if beta==0
C = 50;
pos_t = C*exp((w0*time)*1i) + C*exp(-(w0*time)*1i);
const = C*w0*1i;
vel_t = const*exp((w0*time)*1i) - const*exp(-(w0*time)*1i);
%%Under-damped
elseif beta< w0
delta = arctan(beta/w0);
A = x0/cos(delta);
pos_t = A*exp(-beta*time)*cos(w0*time - delta);
vel_t = -A*exp(-beta*time)*(beta*cos(w0*time- delta) + w0*cos(omega_0*time));
%%Over-damped
elseif (beta> w0)
gamma_minus = beta - sqrt(beta*beta - w0*w0);
gamma_plus = beta + sqrt(beta*beta - w0*w0);
C2 = x0/(1-(gamma_minus/gamma_plus));
C1 = X0 - C2;
pos_t = C1*exp(-gamma_minus*time) + C2*exp(-gamma_plus*time);
vel_t = -gamma_minus*C1*exp(-gamma_minus*time) - gamma_plus*C2*exp(-gamma_plus*time)
else
pos_t = 0;
vel_t = 0;
end
plot(time, pos_t)
  댓글 수: 1
Chris
Chris 2022년 9월 12일
편집: Chris 2022년 9월 12일
First, try setting beta to a single value.
Otherwise, you would need to for loop through each value and test them individually (you can do this once the plots are all working).
Once beta is no longer a vector, run the script and read the error messages.
For starters, arctan() isn't a standard matlab function. You probably want atan().
Then, you'll need to do some element-wise operations. Replace * and / and ^ with .* and ./ and .^
Otherwise, Matlab thinks you want to do linear algebra, when all you want is to apply an operation to each element of a vector/array.
Good luck! If you have a specific error you can't get past, let us know.

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

답변 (1개)

Mathieu NOE
Mathieu NOE 2022년 9월 12일
hello
as mentionned by @Chris there was some missing dots in element wise operations . Others minor bugs included wrong time vector definition (swap end value and number of points) and other minor things - see comments in the code below
then added the for loop with legend showing how your results look like vs beta
hope it helps !
w0= 4*pi;
x0= 100;
beta_list= [0, 2, 5, 10, 20,50]; % list of beta values
% time=linspace(0, 120, 6);
time=linspace(0, 3, 100); % LINSPACE(X1, X2, N) generates N points between X1 and X2.
% plot figure
figure(1),
hold on
for ci =1:numel(beta_list)
beta = beta_list(ci);
%%Undamped Oscillator
if beta==0
C = 50;
pos_t = C*exp((w0*time)*1i) + C*exp(-(w0*time)*1i);
const = C*w0*1i;
vel_t = const*exp((w0*time)*1i) - const*exp(-(w0*time)*1i);
%%Under-damped
elseif beta<= w0 % replaced < with <=
delta = atan(beta/w0);
A = x0/cos(delta);
pos_t = A*exp(-beta*time).*cos(w0*time - delta);
% vel_t = -A*exp(-beta*time).*(beta*cos(w0*time- delta) + w0*cos(omega_0*time));
vel_t = -A*exp(-beta*time).*(beta*cos(w0*time- delta) + w0*cos(w0*time)); % replaced omega_0 with w0
%%Over-damped
elseif (beta> w0)
gamma_minus = beta - sqrt(beta*beta - w0*w0);
gamma_plus = beta + sqrt(beta*beta - w0*w0);
C2 = x0/(1-(gamma_minus/gamma_plus));
% C1 = X0 - C2;
C1 = x0 - C2; % replaced X0 with x0
pos_t = C1*exp(-gamma_minus*time) + C2*exp(-gamma_plus*time);
vel_t = -gamma_minus*C1*exp(-gamma_minus*time) - gamma_plus*C2*exp(-gamma_plus*time);
else
pos_t = 0;
vel_t = 0;
end
plot(time, pos_t)
leg{ci} = (['Beta = ' num2str(beta)]); % legend (char) stored in cell array
end
legend(leg); % display legend once all data is plotted
hold off

카테고리

Help CenterFile Exchange에서 Digital Filter Analysis에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by