FrequencySweep.jl

The FrequencySweep library is tailored towards the calculation of frequency response functions for systems of ordinary differential equations that are difficult or impractical to compute by typical analytical or numerical means (e.g. nonlinear models such as Duffing's equation).

As an illustration, consider Duffing's equation.

\[\ddot{x} + \delta \dot{x} + \alpha x + \beta x^{3} = \gamma \sin(\omega t)\]

To use this library the equation needs to be written in terms of two first order equations.

function duffing_eqn(du, u, p, f, t)
    alpha, beta, gamma, delta = p
    du[1] = u[2]
    du[2] = gamma * sin(f * t) - delta * u[2] - alpha * u[1] - beta * u[1]^3
end

The arguments to the function are as follows:

  • du: This is an output array where the output of each first order equation should be written.
  • u: This is an input array the same size as du providing the current estimates for each state variable.
  • p: A structure of user-defined type that can be used to pass model parameters.
  • f: The oscillation frequency the solver is analyzing. This value will be supplied in the same units as supplied by the user.
  • t: The independent variable, typically time.

Duffing's equation exhibits the jump phenomenon near resonance. The nature of the jump behavior is controlled by the parameters utilized. Setting $\alpha = \gamma = 1$, $\delta = 0.1$, and $\beta = 0.5$ produces the following frequency response.

The frf_sweep function computes an estimate to the frequency response function by solving the system of equations supplied by the user for a predetermined number of oscillations. The solver then computes the amplitude and phase of the last few cycles at end of the solution. These constraints can be defined by the user to better control the process. The solver computes the solution at each supplied frequency point sweeping in the direction defined by the user. The solver uses the final solution point of the previous solution as initial conditions for the next frequency point.

Here is how it all looks in Julia.

f1 = Array(0.5:0.05:2.0)    # sweeping forward in frequency
f2 = Array(2.0:-0.05:0.5)   # sweeping backward in frequency

sweep1 = frf_sweep(
    duffing_eqn,            # function containing the equations to solve
    [0.0, 0.0],             # initial condition vector
    f1                      # frequency points
)
sweep2 = frf_sweep(
    duffing_eqn,
    [0.0, 0.0],
    f2
)

The solution returned from frf_sweep is an M-by-N matrix where M is the number of frequency points and N is the number of equations defined. The solution is complex-valued such that each column supplies both magnitude and phase information for each equation. For the Duffing model utilized above, the frequency response of the first equation is as follows.

A word of caution however. Ensure that only non-zero frequency values are specified as the frequency value determines the length of solution as the solver uses the frequency values to determine the length of the ODE solution.

Function Documentation

FrequencySweep.frf_sweep — Function

rst = frf_sweep(fcn, init, freq)

Description

The frf_sweep function computes the frequency response functions for a system of ordinary differential equations (ODE's) by means of "sweeping" across a range of excitation frequencies. The algorithm solves the system of ODE's at each excitation frequency in the supplied range. It uses the solution of these equations to determine the amplitude and phase of the resulting wavefuorm once it has sufficiently damped out the transient component of the solution. The amplitude and phase are determined by a discrete Fourier transform of the last N cycles of the solution waveform for each ODE.

fcn: The solver requires the user to define a function supplying the ODE's to integrate. This function must be defined as follows.

function fcn(du, u, p, omega, t)
    # User code goes here
end

The arguments to the function are as follows:

  • du: This is an output array where the output of each first order

equation should be written.

  • u: This is an input array the same size as du providing the current

estimates for each state variable.

  • p: A structure of user-defined type that can be used to pass model

parameters.

  • f: The oscillation frequency the solver is analyzing. This value will

be supplied in the same units as supplied by the user.

  • t: The independent variable, typically time.

init: The user must also supply an initial set of initial conditions for the solution of the ODE's. This set of initial conditions is only used for the solution to the first frequency value in the range as the solver will update the initial condition vector for subsequent frequency value solutions.

freq: The user must also supply a range of frequency values at which to compute the frequency responses. Ensure that there are no zero-valued frequencies in the range as the solver does not accept such a value. The solver utilizes the frequency value as a method of determining the length of the ODE solution; therefore, smaller frequency values result in longer solutions. A note about units: It is recommended that the frequency values be supplied with units of radians per unit time; however, whatever units are provided will be the same units used for the value passed to fcn.

The function outputs an M-by-N matrix where M is the number of frequency points and N is the number of ODE's being solved. The result is complex-valued such that each column of the matrix describes both the amplitude and phase of each equation.

The user may supply a set of optional parameters to control both the ODE solver and the calculation of amplitude and phase for each equation. The optional parameters are as follows.

ODE Solver Optional Arguments

    • alg_hints: * Any hints to pass to the ODE solver (e.g. :stiff).
    • reltol: * THe relative solution tolerance.
    • abstol: * The absolute solution tolerance.
    • dtmax: * The maximum allowable step size in the independent variable.
    • maxiters: * The maximum number of iterations to allow.
    • modelparams: * A container for arguments to pass to fcn.

Additional Optional Arguments

    • pointspercycle: * Defines the number of solution points for each cycle

of the solution.

    • cyclecount: * Defines the number of solutions to use in the calculation

of the amplitude and phase.

    • transientcount: * Defines how many cycles to let pass for the transient

to sufficiently decay.

    • windowsize: * Defines the size of window to utilize for the discrete

Fourier transform.

    • windowfunction: * Defines the window function to utilize for the

discrete Fourier transform. The function should take the same form as the windowing functions from the DSP.jl library.

Example

The following example illustrates how to use the function to obtain both an ascending and descending solution to Duffing's equation.

# Duffing's Equation
function duffing_eqn(du, u, p, omega, t)
    du[1] = u[2]
    du[2] = gamma * sin(omega * t) - delta * u[2] - alpha * u[1] - beta * u[1]^3
end

# Frequency Points
f1 = Array(0.5:0.05:2.0)
f2 = Array(2.0:-0.05:0.5)

# Perform the sweep
sweep1 = frf_sweep(
    duffing_eqn,
    [0.0, 0.0],
    f1
)
sweep2 = frf_sweep(
    duffing_eqn,
    [0.0, 0.0],
    f2
)