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

2 次查看（过去 30 天）
Sameer
Sameer 2017-11-20
回答： David Goodmanson 2017-11-21
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 Center 和 File 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