필터 지우기
필터 지우기

Why does the following simple Harmonic oscillator (without damping) behave as a damped harmonic oscillator?

조회 수: 1 (최근 30일)
function resonance
omega = 100; %
b = 0.0; %
A = 0.0; % driving amplitude per unit mass
omega0 = 0; % driving frequency
tBegin = 0; % time begin
tEnd = 80; % time end
x0 = 0.2; % initial position
v0 = 0.8; % initial velocitie
a = omega^2; % calculate a coeficient from resonant frequency
% Use Runge-Kutta 45 integrator to solve the ODE
[t,w] = ode45(@derivatives, [tBegin tEnd], [x0 v0]);
x = w(:,1); % extract positions from first column of w matrix
v = w(:,2); % extract velocities from second column of w matrix
plot(t,x);
title('Damped, Driven Harmonic Oscillator');
ylabel('position (m)');
xlabel('time (s)');
% Function defining derivatives dx/dt and dv/dt
% uses the parameters a, b, A, omega0 in main program but changeth them not
function derivs = derivatives(tf,wf)
xf = wf(1); % wf(1) stores x
vf = wf(2); % wf(2) stores v
dxdt = vf; % set dx/dt = velocity
dvdt = -a * xf - b * vf + A * sin(omega0*tf); % set dv/dt = acceleration
derivs = [dxdt; dvdt]; % return the derivatives
end
end
%

채택된 답변

David Goodmanson
David Goodmanson 2017년 11월 21일
Hi Sameer,
I believe this is a standard numerical accuracy issue. Here is some shortened code that is basically the same as yours. If you run it as is, you get similar behavior to yours.
[t x] = ode45(@(t,x) fun(t,x), [0 2000], [1 0]);
plot(t,x(:,1))
function dxdt = fun(t,x)
dxdt = [0 1;-5 0]*x;
end
If you replace the first line with the two lines
options = odeset('reltol',1e-5);
[t x] = ode45(@(t,x) fun(t,x), [0 2000], [1 0],options);
then with the tighter tolerance the amplitude is basically constant all the way across. Of course it is going to take a bit longer. There is also an abstol you can vary, and it appears that the default values for reltol and abstol are 1e-3 and 1e-6. The odeset function without any arguments gives you a list of what all the optional parameters are.

추가 답변 (0개)

카테고리

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