FRF for 2DOF system

[Paw]

Dear all,

I would like to generate a nonlinear frequency response function(NFRF) for a 2DOF system.
A NFRF is a plot of the stationary amplitude for varying excitation frequncy.

Example for the 1DOF duffing oscillator:
To find the stationary response for the equation x''+dx'+bx+ax^3=gcos(omega*t), we use an expansion technique, here multiple scales, and end up with modulation equations governing amplitude and phase around a selected resonant frequency.
These eqs. are only valid in the vicinity of the selected resonant frequency, thus we need to derive a new set of modulation equations for each resonant frequency of interest( sub- and super harmonic)

The modulation equations are solved used continuation, using the frequency as the free variable.

Here is an example for the duffing system exited near the first primary resonance.

#!/usr/bin/python2

from PyDSTool import *
import matplotlib.pyplot as plt

"""
Multiple scales analysis of a duffing eq.
Excited near primary external resonance
See eq. 3.188 in 'Vibration and Stability', Jon J Thomsen
"""

# plot settings
plt.close('all')

# fixed simulation parameters
maxNumPoints = int(1e4)
maxStepSize = 1e-2
stepSize    = 1e-3
minStepSize = 1e-5

# fixed physical parameters
O = 0.0   # starting point for excitation
o = 1.0   # omega0 - linear eigenvalue
g = 0.5   # gamma - nonlinear parameter
q = 0.2   # q - external forcing level
b = 0.05  # beta - damping

# fixed equilibrium curve
# use that O = o + eps*s <=> s = O-o
ustr = '-b*o*u+q/(2*o)*sin(v)'  # amplitude
vstr = '(O-o)*u-3*g/(8*o)*u*u*u+q/(2*o)*cos(v)'  # phase

# predictor-corrector arguments
PCargs = args(name='EQ1', type='EP-C')
PCargs.freepars = ['O']
PCargs.StepSize = stepSize
PCargs.MaxNumPoints = maxNumPoints
PCargs.MaxStepSize = maxStepSize
PCargs.MinStepSize = minStepSize
PCargs.LocBifPoints = ['LP', 'BP', 'B']
PCargs.StopAtPoints = ['B']
PCargs.verbosity = 2
PCargs.SaveEigen = True
# Note: 'MoorePenrose' is much more stable than'Natural' corrector:
PCargs.Corrector = 'MoorePenrose'

DSargs = args(name='duffing')
DSargs.pdomain = {'O': [0, 4.0]}  # domain boundaries
DSargs.pars = {'O': O, 'o': o, 'g': g, 'q': q, 'b': b}
DSargs.varspecs = {'u': ustr, 'v': vstr}
DSargs.ics = {'u': 0.01, 'v': 1.0}
testDS = Generator.Vode_ODEsystem(DSargs)

# loop beta's
bvec = [0.20, 0.10, 0.05]

fig1 = plt.figure(num=1)
for b in bvec:
    testDS.set(pars={'b': b})
    PyCont = ContClass(testDS)

    PyCont.newCurve(PCargs)

    print('Computing equilibrium curve...')
    start = clock()
    PyCont['EQ1'].forward()
    print('done in %.3f seconds!' % (clock()-start))
    PyCont.display(('O', 'u'), stability=True, linewidth=0.5, label='beta=%f' % b)

# Plots
plt.legend(loc=2)
plt.grid(True)
plt.title('Frequency response')
ymin, ymax = plt.ylim()
plt.ylim([0, ymax])
fig1.show()

plt.show()

Giving the attached plot, where each curve is for a different level of damping.
The y-axis show the amplitude and the x-axis is detuning away from the resonant frequency.
For low levels of damping, we see birfucation indicated with stability.

That was a long introduction.
My question is now:

Can I obtain the same NFRF without deriving the modulation equations, eg. just plug the ODE into PyDSTool?

How do I get the NFRF for a two DOF system, for instance

x1'' + (0.02x1' - 0.01x2') + (2x1 - x2) + 0.5x1^3 = Fcos(omega*t)
x2'' '(0.002x2' - 0.01x1') + (2x2 - x1) = 0

?

Best regards,
Paw

[Rob Clewley]

Off the top of my head, I'd say you have one primary option in the
absence of modulation equations. Bearing in mind the potential for
numerical loss of precision near regions in which solutions "blow up"
in some way, you could simply forward integrate to the steady-state
response frequency over an adaptive 1D mesh of forcing frequencies (or
whatever parameter). You need to define a reasonable tolerance for
defining "steady state" oscillatory behavior that would program a
terminating event. You would also use an appropriate integration
scheme for potentially multi-scale behavior in x1 and x2 (or maybe
adapt some integration parameters as your code sense proximity to
resonance). You'd have to interpolate over the mesh, and you wouldn't
be able to calculate the unstable branches (although you might try
infering them under some assumptions if you have a refined enough
mesh).

This kind of general algorithmic question would
probably be well suited for one of the specialized Stack Exchange
forums on computational methods / numerical algorithms. There might
well be some more sophisticated ways you could approach this.

If you find an interesting approach that works, please report back here with some code!