Solve a 2nd Order ODE: Damped, Driven Simple Harmonic Oscillator
This example builds on the first-order codes to show how to handle a second-order equation. We use the damped, driven simple harmonic oscillator as an example:
In a second order system, we must specify two initial conditions. In the plot shown below we choose
We assume MKS units, but this is unimportant for our discussion. The system is resonant when the driving frequency matches the natural frequency. We choose our parameters near resonance:
The following figure shows the result. The dynamics show initial transient behavior which gives way to resonant oscillations.
% resonance.m % Numerically integrate second-order ODE: Damped, driven harmonic oscillator function resonance omega = 1; % resonant frequency = sqrt(k/m) b = 0.2; % drag coeficient per unit mass A = 0.1; % driving amplitude per unit mass omega0 = 1.2; % 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
We compare the differences between the first order code (logisticV1.m) this code. Because the derivatives routine must reference both position and velocity, we 'package' these variables into a single array we call wf. Inside the derivatives function, we then 'unpack' the variables of interest (position xf and velocity xf in this case):
xf = wf(1); % wf(1) stores x vf = wf(2); % wf(2) stores v
Notice that because variables x and v are used in the main program, we choose different variable names in the function: xf and vf. Otherwise, the main program variables x and v will be overwritten and complete monkey chaos will ensue. We stick an 'f' on each of the variables to indicate they reside in our function...
Next, we calculate the derivatives. Because we are solving a second order system, we need to return two derivatives:
dxdt = vf; % set dx/dt = velocity dvdt = -a * xf - b * vf + A * sin(omega0*tf); % set dv/dt = acceleration
In a second-order system the first derivative dxdt is just be set to the velocity variable. The velocity derivative dvdt contains the interesting information defining the acceleration.
We package up our derivatives and stuff them in a new array we call derivs:
derivs = [dxdt; dvdt]; % return the derivatives
This variable derivs must be the same variable name used in defining the function:
function derivs = derivatives(tf,w)
Now back to the main program. The ode45 function returns our packaged dependent variables and the independent variable, which we call w and t. As before, we need to unpack our variable w to retreive the position x and velocity v data:
[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
That's it! The variables x and v contain the solutions we desired. We can now plot them or use them in in other calculations.