#52 Amplitude definitions

[Marc]

Where can I find more information about the different available choices of amplitudes when plotting the bifurcation diagram, and where can I define my own amplitude if needed?

In GettingStartedwithCOCO.pdf (11/08/23), one finds

• x (p. 6)
• ||x||_2 : Euclidian norm (p. 17)
• ||x||_{2,MPD} : "suitably defined scalar magnitude for the periodic orbits" (p. 17).

Thank you

[Harry Dankowicz]

Hello Marc,

Thank you for exploring ways to extend COCO for your own needs. That is exactly what the platform was built to support. I am glad that you are feeling inspired to push the envelope.

The short answer to your question is as follows.

When constructing a continuation problem for bifurcation analysis of equilibria using the 'ep' toolbox, the constructor ep_add_bddat can be invoked to append a column to the output cell array consisting of scalars computed from the values of the state variable vector x and parameter vector p associated with each equilibrium point.

You can find an example of such a construction in Section 7.1 of GettingStartedwithCOCO.pdf in the command

prob = ep_add_bddat(prob, '', 'svds', @(d,x,p) min(svds(feval(F('x'),x,p))));

This appends a column 'svds' to the output cell array with the smallest singular value of the Jacobian of the vector field with respect to the state variable vector at each equilibrium point.

Similarly, when constructing a continuation problem for bifurcation analysis of periodic orbits using the 'po' toolbox, the function po_add_bddat can be used to append a column to the output cell array consisting of scalars computed from the values of the basepoint discretization xbp , the initial time T0 and period T, and the parameter vector p associated with each periodic orbit.

You can find an example of such a construction in Section 7.3 of GettingStartedwithCOCO.pdf in the command

prob = po_add_bddat(prob, '', 'cx1', @fourier, 'data', struct('n', 1));

This uses the function fourier defined earlier in the problem to append a column 'cx1' to the output cell array with the Fourier coefficient of the first harmonic evaluated for each periodic orbit.

A longer answer will follow as a separate post.

I hope this helps.

/Harry

[Harry Dankowicz]

Here's a longer answer:

For an arbitrary COCO-compatible continuation problem, columns may be added to the output cell array by defining slot functions that respond to the pre-defined 'bddat' signal and perform calculations using abstract and numerical data associated with each solution.

In the October 26, 2023, release of COCO, one such slot function is defined in the ep_add constructor in the 'ep' toolbox and is shown below.

function [data, res] = bddat(prob, data, command, varargin)
%BDDAT   Append EP bifurcation data to BD.

res = {};
switch command
  case 'init'
    xid   = data.ep_bddat.xid;
    nrmid = sprintf('||%s||_2', xid);
    maxid = sprintf('MAX(%s)', xid);
    minid = sprintf('MIN(%s)', xid);
    res   = { nrmid, maxid, minid, xid };
  case 'data'
    eqn   = data.ep_eqn;
    chart = varargin{1}; % Current chart
    uidx  = coco_get_func_data(prob, eqn.fid, 'uidx');
    u     = chart.x(uidx);
    x     = u(eqn.x_idx);
    res   = { norm(x,2), x, x, x };
end

end

When the 'bddat' signal is emitted with command equal to 'init', the function returns a cell array with four string labels that become headers of four columns added to the output cell array.

When the 'bddat' signal is emitted with command equal to 'data', the function returns a cell array containing the Euclidean norm of the state vector, and three copies of the state vector.

A similar construction is found in the following slot function defined in the po_add constructor in the 'po' toolbox.

function [data, res] = bddat(prob, data, command, varargin)
%BDDAT   Append PO bifurcation data to BD.

res = {};
switch command

  case 'init'
    xid   = data.po_bddat.xid;
    nrmid = sprintf('||%s||_{2,MPD}', xid); % mean + deviation plotting measure
    maxid = sprintf('MAX(%s)', xid);
    minid = sprintf('MIN(%s)', xid);
    res   = { nrmid, maxid, minid };

  case 'data'
    chart = varargin{1};
    [fdata, uidx] = coco_get_func_data(prob, data.cid, 'data', 'uidx');
    maps = fdata.coll_seg.maps;
    mesh = fdata.coll_seg.mesh;
    xbp  = reshape(chart.x(uidx(maps.xbp_idx)), maps.xbp_shp);
    iw   = (0.5/maps.NTST)*mesh.kas2*mesh.wts2*maps.W;
    xav  = sum(reshape(iw*xbp(:), maps.x_shp),2);
    dev  = bsxfun(@minus, xbp, xav);
    ndv  = reshape(iw*(dev(:).^2), maps.x_shp);
    mpd  = norm(xav)+sqrt(sum(sum(ndv,2),1));
    res  = { mpd , max(xbp,[],2) , min(xbp,[],2) };

end

end

The constructors ep_add_bddat and po_add_bddat referred to in the previous response were developed only recently in order to lower the barrier for users to extend the default functionality of the 'ep' and 'po' toolboxes in this way, while bypassing some of the repetitive syntax. The functions ep_HB_add_bddat and ep_SN_add_bddat allow output data to be computed not only from the values of x and p but also using the eigenvectors and eigenvalues of the corresponding Jacobian of the vector field.

[Marc]

Thank you for this detailed information.
I have not yet understood the syntax when using ep_add_bddat and po_add_bddat (for instance, why the argument d in the example of svds in your first answer above)?
Very practically speaking, I have calculated the continuation of a vector field as function of (position x, speed xdot). I now would like to draw the bifurcation diagram using

• either abs(max(x)-min(x))
• or abs(max(xdot)-min(xdot))
  as amplitude.
  When using po_add_bddat, how to indicate to use position x or speed xdot in the calculation?
  prob = po_add_bddat(prob, '', 'myamplitude', @(??) abs(max(??)-min(??)));

[Harry Dankowicz]

Hi Marc,

By design, the function defined in the fourth argument of ep_add_bddattakes three input arguments, the last two of which correspond to the state vector and vector of problem parameters for an equilibrium point. The first argument is not used in the example that you are referring to (but is required syntactically). We could have written

prob = ep_add_bddat(prob, '', 'svds', @(~,x,p) min(svds(feval(F('x'),x,p))));

to make that clear.

By analogy, the function defined in the fourth argument of po_add_bddat takes five arguments, the last four of which correspond to the basepoint discretization, initial time, period, and vector of problem parameters for a periodic orbit. The first argument is used in the example shown in my first response and contains a struct with fields provided by the po toolbox as well as the user-provided field n which here equals 1. The function @fourier makes use of this information.

As to your practical questions, I am interpreting them as follows: you are performing continuation of periodic orbits in a single-degree-of-freedom mechanical system (with position x and velocity xdot). For each periodic orbit, you would like to calculate, say, the absolute value of the difference between the maximum and minimum values of the position along the orbit and/or similarly for the velocity.

Of course, since we are continuing a discrete approximation of the orbit, you should first decide whether you are ok with limiting consideration to maximal and minimal values among the basepoints or wish to find these extreme values along the assumed continuous, piecewise-polynomial interpolant.

If it is the former (simpler), then you have two options. The first is to use bifurcation data that is already computed by the po toolbox, as in

coco_plot_bd('po_run', 'p1', 'p2', {'MAX(x)' 'MIN(x)'}, @(x,y) abs(x(1,:)-y(1,:)));

for the position and

coco_plot_bd('po_run', 'p1', 'p2', {'MAX(x)' 'MIN(x)'}, @(x,y) abs(x(2,:)-y(2,:)));

for the velocity.

Alternatively, at the time of construction, do

prob = po_add_bddat(prob, '', {'xamp' 'xdotamp'}, @marc);

where

function y = marc(data, xbp, ~, ~, ~)

mps = data.coll_seg.maps;
xbp = reshape(xbp, mps.xbp_shp)';
y = num2cell(abs(max(xbp)-min(xbp)));

end

and then

coco_plot_bd('po_run', 'p1', 'p2', 'xamp');

for position and

coco_plot_bd('po_run', 'p1', 'p2', 'xdotamp');

for velocity.

If you would like to find the maximal and minimal positions and velocities along the interpolant, then you will need to reconstruct this using the discretized data. This is obviously more complex. The demos in Section 7.4 of the Getting Started tutorial is a good place to start.

I hope that helps.

/Harry

[Marc]

Dear Harry

Thank you. For the practical case of a SDOF system with state (position x, velocity v), and where we only consider the values that were calculated during continuation, your following suggestion :

coco_plot_bd('po_run', 'p1', 'p2', {'MAX(x)' 'MIN(x)'}, @(x,y) abs(x(1,:)-y(1,:)));

and

coco_plot_bd('po_run', 'p1', 'p2', {'MAX(x)' 'MIN(x)'}, @(x,y) abs(x(2,:)-y(2,:)));

enables to plot the maximum amplitude of the position (resp. velocity).

I have a follow-up question: how to handle a 2-DOF system with state
(position x1, position x2, velocity v1, velocity v2)
where position and velocity are calculated from the modal contributions:
position x = x1phi1 + x2phi2 (phi1,2 being modal shapes)
velocity v = v1phi1 + v2phi2
?
Obviously, for the example of the position x,

coco_plot_bd('po_run', 'p1', 'p2', {'MAX(x)' 'MIN(x)'}, @(x,y) abs((phi1*x(1,:)+phi2*x(3,:)) - (phi1*y(1,:)+phi2*y(3,:)));

is wrong as it considers the maxima (resp. minima) for x1 and x2 separately:
abs(phi1MAX(x(1,:))+phi2MAX(x(3,:))
instead of
abs(MAX(phi1x(1,:)+phi2x(3,:))

Best regards
Marc

[Marc]

It seems that the asterisk is interpreted as italic text. I mean
position x = x1 . phi1 + x2 . phi2
abs(phi1 . MAX(x(1,:)) + phi2 . MAX(x(3,:)))
abs(MAX(phi1 . x(1,:) + phi2 . x(3,:)))

[Harry Dankowicz]

Hello, Marc.

You have several options.

For example, you can create a new bddatfunction that computes the quantities of interest, thereby allowing you to use them with coco_plot_bd. Just extract the xbp array as in the example above and then perform the necessary additions and application of the maximum and absolute value functions.

You can also add a monitor function to the continuation problem that does the same thing. In this case, you will need to include a remesh function, as may be found in several demos.

I hope that helps.

Harry

[Marc]

Hello Harry

Thank you. I now have understood how to use a bddat function during the continuation run to store quantities of interest to be used during the plotting, and am adapting your example above.
For this, I also need to access the current parameter values inside the bddat function.
I checked the content of data.coll_seg but did not find this, although data.coll_seg.maps has the field pdim (dimension of parameter vector).
Related to your comments above, I suppose that the parameter vector p has to be passed to the bddat function, but which argument of the 5?

function y = marc(data, xbp, ~, ~, ~)

[Marc]

Got it! It seems to be
function y = marc(data, xbp, ~, ~, p)

[Harry Dankowicz]

Yes, the arguments are the continuation variables associated with an instance of the trajectory zero problem, xbp, T0, T, and p.

Best,

Harry

[Marc]

Perfect! This question can be closed.