Dear Harry, many thanks for the soon reply, the code have been now adapted according to your advice. Kind regards and thanks, Joseph
Thank you, Joseph, for the question. The most recent release of COCO provides new general-purpose support for generic branch switching (that was originally motivated by and is compatible also with analysis of symmetry-breaking bifurcations in dynamical systems with symmetries). As a result, constructors of the form ===BP2TYPE have all been deprecated in favor of the corresponding ===TYPE2TYPE constructor and the assignment of 'switch' to the 'branch' property of the 'cont' settings. You can find...
Dear Harry, with the kind help of one of your former PhD students we were able to use the adjt_isol2hspo constructor to tackle the problem described in this post. However, during the continuation runs we get the following warning: Warning: Use of ADJT_BP2COLL is deprecated. Reset 'branch' property In adjt_BP2coll (line 45) In adjt_BP2bvp (line 58) In adjt_BP2hspo (line 53) And similarly for ADJT_BP2HSPO, ODE_BP2COLL and ODE_BP2BVP. Does the warning mean that the code may not run in future COCO releases?...
I am thinking of an example where the family of solutions lies on a circle with two fold points at either extreme. The default atlas algorithm then keeps going around and around, since it has no memory of the parts of the atlas that have been covered previously. An alternative is the atlas-kd algorithm, which does retain such a memory (it is similar to the advancing front algorithm in Recipes for Continuation). To switch atlas algorithm, use coco_set(prob, 'cont', 'atlas', 'kd'). Best, Harry
Ok, thank you. That looks quite complex but I will have a look. Besides the automatic stop, is there a more elegant way to prevent Coco to compute single branches over and over again without having to adjust PtMX in a 'correct way'?
Ah, I understand now that 'BP' in your question refers to the event string used by the atlas algorithm for branch points, not an event type of 'BP', i.e., a boundary. I am proposing a path below, which I believe addresses your need, but I haven't tried it and would probably need to experiment a bit to make it work. You would like the code to track the number of events of a certain type that have been encountered during continuation and take an action that depends on this number, e.g., terminate continuation...
Hi Harry, I already tried to understand the explanation of the coco_add_event function but I still have difficulties. I'm using the po toolbox. During continuation Coco detects an FP/BP. It continues at that point which is fine. But now, I want Coco to stop when such a FP/BP is detected. More detailed, when a FP/BP is detected the first time, it should continue and when it detects a second, it should end. Best, Jonathan
Hi Jonathan, As you can see by typing help coco_add_event on the command line, event signatures come in several different forms, depending on the desired function. The 'BP' or 'boundary' event type defines a computational boundary beyond which continuation will not proceed, although it may still go in other directions from any available chart. The 'MX' or 'terminate' event type, on the other hand terminates continuation entirely. I can't tell from your question whether you are relying on an event...
Stop Continuation at Special Point
Perfect, thank you!
Hi Marc, Yes, the problem is that when the solution type is a periodic orbit, the plotting routine by default assumes the existence of a po.test.USTAB column in the corresponding bifurcation data array. You can turn this off by changing thm = struct('special', {{'EP'}}); to thm = struct('ustab', '', 'special', {{'EP'}}); I hope that helps. /Harry
Or first a follow-up question, still in the file coco/po/examples/hopf/demo.m, where I have added following lines to plot po_run: figure(2); clf; thm = struct('special', {{'EP'}}); coco_plot_bd(thm, 'po_run', 'p1', 'p2', '||x||_{2,MPD}') When po_run has been calculated with 'bifus', 'off', the plotting function gives following error: Error using coco_bd_col>coco_get_col coco_bd_col: column not found: 'po.test.USTAB' Error in coco_bd_col (line 32) col = coco_get_col(bd, name, cat_flag); Error in coco_plot_bd...
Hi Thanks for confirming the syntax. A quick test with coco/po/examples/hopf/demo.m confirmed that order matters. prob = coco_prob(); %prob = coco_set(prob, 'po', 'bifus', 'off'); %test 1: OK prob = ode_HB2po(prob, '', 'ep_run', HBlab); %prob = coco_set(prob, 'po', 'bifus', 'off'); %test 2: too late fprintf(... '\n Run=''%s'': Continue periodic orbits from point %d in run ''%s''.\n', ... 'po_run', HBlab, 'ep_run'); bd3 = coco(prob, 'po_run', [], 1, 'p1', [-1 1]); The error was that I tried to add...
Yes, both syntaxes are correct. This is the meaning of 'log|on', as in logical (true/false) or on/off. A cursory glance suggests that this should work for all po constructors and I have no reason to believe this is not the case. If you are willing to share your script, I can have a look. Best, Harry
*all po constructors
Disable detection of bifurcations
Amplitude definitions
Perfect! This question can be closed.
Yes, the arguments are the continuation variables associated with an instance of the trajectory zero problem, xbp, T0, T, and p. Best, Harry
Got it! It seems to be function y = marc(data, xbp, ~, ~, p)
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...
Thank you so very much Dear Prof. Harry for your suggestion. I will surely look into it and let you know the update.
Hello Azmir, I recommend that you compare your implementation to one of the demos that include continuation of periodic orbits from a Hopf bifurcation. Without checking your code, I am inclined to suspect a typo somewhere in your Jacobian. As always, you need to be concerned about appropriate use of vectorization, differentiability, and division by 0. I hope that helps. Harry
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...
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,:)))
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...
How to solve the issues with "Warning: Matrix is singular"
January 28, 2025, release of COCO
fixes
Corrected typos in old covering demos and updated symbolic demos for compatibility with new branch switching algorithm
Corrected bug in demo to agree with 'bistable' demo, specifically use of -no-var option in ode_po2po.
fixes
fixes
fixes
fixes
Hi Harry, thank you, I think I'm understanding the theory behind your answer. Jonathan
Added remark regarding event detection using non-embedded monitor functions. This is not robust to step-size selection, since the step size algorithm only considers embedded variables.
Add handling of discrete monitor functions such that no value needs to be provided, contrary to continuous monitor functions
Fixed bug when discrete monitor function takes the value 0.
Adding support for discrete monitor functions
Adding support for discrete monitor functions and event detection
Corrected typo in Tutorial
Updated atlas tutorial to describe 'branch' optional setting.
Remove temporary UST1 solution type.
Updated tutorials for 'ep', 'coll', 'ep', and Getting Started to agree with code and demo updates to branch handling.
New demo for showing how core handles regularization of linearly dependent constraints.
Major update to branch handling for core and 'ep', 'coll', and 'po' toolboxes, including updates to all relevant demos (only CORE tutorial has been updated). In new version, branch handling depends on a 'branch' property of the 'cont' settings. This property takes three distincts values. CurveSegment and the default value are designed to ensure backward compatibility.
Hi Jonathan. Interesting question. It is true that the 'ep' toolbox only handles autonomous vector fields. Making a non-autonomous vector field autonomous through the introduction of a phase variable, however, violates the condition for an equilibrium, since the rate of change of the phase variable is nonzero. This is the reason for a failure to converge, since there is no solution to converge to. The 'ep' toolbox simply isn't designed for this case. My recommendation is to use calls to coco_add_funcand...
How to set perturbation in ode_PD2po
Just to make your question clearer: a MX error is raised at STEP -1, indeed. Furthermore, there is a Warning: Matrix is singular to working precision.
Just to make your question clearer: a MX error is raised at STEP -1, indeed. Furthermore, there is a Warning: Matrix is singular to working precision.
Just to make your question clearer: a MXerror is raised at STEP -1, indeed. Furthermore, there is a Warning: Matrix is singular to working precision.
Thank you for the information. I was confused. I finally understood how to use ode_PD2po with success in order to find the period-doubling branch between two PDs. This question can be closed.
I was confused. I finally understood how to use ode_PD2po with success in order to find the period-doubling branch between two PDs. This question can be closed.
Continuation of Equilibria / Formulation of Autonomous System
Frequency Response Plots using COLL/PO toolbox
Updated po utilities and demos to handle tangent vectors correctly. Fixed bug in mvdP files.
Updated coll restart constructors, utilities, and demos to handle tangent vectors correctly.
Updated ep_construct_tst with event handler to store SN eigenvector also when not detected as fold point: Modified ep_read_solution accordingly. Updated constructors to properly handle tangent vector stored with solution file. Tentatively updated restart constructors to allow for a -no-pars option in the event that one wishes to reconstruct without adding continuation parameters. This is necessary if constructing multiple instances from stored solutions.
Ok nice, I figured out how to solve that problem. Thanks a lot!
Ok nice, I figured out how to solve that problem!
Hi Jonathan, The tbp array contains the time stamps for the basepoint mesh. Assuming that you use adaptation (i.e., with NAdapt greater than 0), this is not uniformly distributed, as the algorithms redistribute the mesh to capture areas of greater variability in the solution. To obtain a uniformly distributed time array, use interpolation. I recommend interpolation using the continuous piecewise-smooth polynomial approximant that COCO assumes when discretizing the trajectory segment. Such interpolation...
Hi Harry, thank you. My next question is regarding the tbp vector. I recognized that the values are not equally distributed (which is in terms of efficency I guess). Is it somehow possible to get a equally distributed time vector? Best, Jonathan
Hello Jonathan, There are two ways to have user-defined data appear in the bd array: either by appending such data using a slot function that responds to the bddat signal or by adding continuation parameters. For example, the po toolbox defines several slot functions for the bddat signal that add the Floquet multipliers and norms of the solution to the bd array. The function po_add_bddat (as explained in the GettingStartedWithCOCO tutorial) can be used to create such a slot function. The function...
Hi Harry, I was using the coco_plot_bd feature for plotting amplitude over excitation frequency in the past and now I'm wondering if the same is possible for the phase response? I've found no fast way using a column in the bd array. Best, Jonathan
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...
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...
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...
coll and po: corrected use of VAR for restarting from branch of Hopf bifurcations.
CurveSegment: support for branch switching
covering/toolbox/@atlas_1d: introduced support for branch switching
Dear Harry, Thank you for your help. The first question is answered by setting PTMX. For the second question, I just rewrote the function coco_add_func. Best Regards, Bo
Dear Harry, Thank you for your help. The first question is answered by setting PTMX. For the senod question, I just rewrote the function coco_add_func. Best Regards, Bo
Dear Harry, Thank you for your help. The first question is answered by setting PTMX. For the senod question, I wil try again. Best Regards, Bo
Hi Bo, Continuation proceeds as long as solutions remain inside the computational domain AND as long as the maximum number of points (PtMX) has not been reached. The default value for PtMX is 100 in both directions along a one-dimensional solution manifold. You can change this to other values, including different values in each of the two directions, by using a command like prob = coco_set(prob, 'cont', 'PtMX', [200 50]); As to how much of the solution manifold that is covered by a particular setting...
Dear Harry, Thank you for your reply. I am not sure I understand what you mean. why can we not calculate the entire frequency response curve at conce? Reason 1: If I take the excitation period T=2pi/5 as the initial parameter, so that in the jump region of the frequency response curve, the actual excitation period T may be not continue with the frequency response curve (the excitation period has changed, not continued). Instead, we need to constrain the curve from the other point, such as T=2pi/20....
Hi Bo, As I read your questions, I think it is important to separate out different components of the code. Let's start with coco_plot_bd('freq_resp', 'T', @(T) 2*pi./T, '||po.orb.x||_{L_2[0,T]}') This command is executed once a run of the coco entry-point function is complete. It assumes that the run was named 'freq_resp' and that there was a continuation parameter named 'T' (or that a user-defined function responding to the bddat signal created a column with this header). It also assumes that the...
Dear Harry, Thank you for your reply. I have changed the codes based on your comments. I have a few questions I would like to ask for your help. In the line 23 of demo_po.m, cont_args = {1, {'po.period' 'T'}, [2*pi/200 2*pi/2]}; it appears that the period of the system is limited between 2 pi/200 and 2 pi/2 (gen_sym_po.m and demo_po.m) . The result is shown in Figure 1. If I change the following to read p0=[2*pi/20; 1]; [~,x0]=ode45(@(t,x)funcs{1}(t,x,p0),[0 20*pi],[0.01; 0; 0; 0; 0; 0; 0; 0; 0;...
Dear Harry, Thank you for your reply. I have modified the codes based on your comments. I have a few questions to ask you for help. In the line 23 of demo_po.m, cont_args = {1, {'po.period' 'T'}, [2*pi/200 2*pi/2]}; it appears that the period of the system is limited between 2 pi/200 and 2 pi/2 (gen_sym_po.m and demo_po.m) . The result is shown in Figure 1. If I change the following to read p0=[2*pi/20; 1]; [~,x0]=ode45(@(t,x)funcs{1}(t,x,p0),[0 20*pi],[0.01; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]); [t0,x0]=ode45(@(t,x)funcs{1}(t,x,p0),[0...
Dear Harry, Thank you for your reply. I have modified the codes based on your comments. I have a few questions to ask you for help. In the line 23 of demo_po.m, cont_args = {1, {'po.period' 'T'}, [2*pi/200 2*pi/2]}; it appears that the period of the system is limited between 2 pi/200 and 2 pi/2 (gen_sym_po.m and demo_po.m) . The result is shown in Figure 1. If I change the following to read p0=[2*pi/20; 1]; [~,x0]=ode45(@(t,x)funcs{1}(t,x,p0),[0 20*pi],[0.01; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]); [t0,x0]=ode45(@(t,x)funcs{1}(t,x,p0),[0...
Dear Bo, Thank you for your question. The MX label indicates a failure of the algorithm to converge to a solution within the desired tolerance. This is a warning to you, not necessarily that there is anything wrong with your continuation settings, but perhaps that there may be something amiss or at least needing attention with your problem. For example, I note that the initial value for the first problem parameter is (line 4 of demo_po.m) 2*pi/7. From line 5 of gen_sym_po.m, I conclude that p1 is...
Dear Harry, I have a problem plotting the frequency amplitude curve. I am following the example for plotting the frequency amplitude curve of bistable, the following codes are constructed. Unfortunately, MX is displayed. I do not see in the documentation how to adjust the continuation settings to avoid MX. The other equation is about the frequency amplitude curve in relation to the excitation frequency Omega https://sourceforge.net/p/cocotools/tickets/51/#021b. Are there details on how to plot this?...
Dear Harry, I have a problem plotting the frequency amplitude curve. I am following the example for plotting the frequency amplitude curve of bistable, the following codes are constructed. Unfortunately, MX is displayed. I do not see in the documentation how to adjust the continuation settings to avoid MX. The other equation is about the frequency amplitude curve in relation to the excitation frequency Omega https://sourceforge.net/p/cocotools/tickets/51/#021b. Are there details on how to plot this?...
Dear Harry, I have a problem plotting the frequency amplitude curve. I am following the example for plotting the frequency amplitude curve of bistable, the following codes are constructed. Unfortunately, MX is displayed. I do not see in the documentation how to adjust the continuation settings to avoid MX. The other equation is about the frequency amplitude curve in relation to the excitation frequency Omega https://sourceforge.net/p/cocotools/tickets/51/#021b. Are there details on how to plot this?...
Hi Marc, I am not sure what you mean by "one of the period-doubled branches." If I understand the situation correctly, there is only one branch that connects the two bifurcation points. It may show up as having two parts because of sampling at the original period, but it's the same orbit sampled twice (as when we plot intersections with a Poincare section). If you would like to sample twice per period, you would need to create an appropriate monitor function. For example, such a function can use...
coco_plot_bd many points
Indeed, Marc, there are surely (as you point out) other ways to address the problem of too many data files. I am encouraging of users developing alternatives to, or modifications of, the routines that ship with the release. I do want to stress, however, that the situation you describe is quite unusual, and very possibly a symptom of something amiss. I am assuming that you are tracking this down in your code. Best, Harry
Hello Harry Thanks for your comments. In any case, if plotting would become too slow because too many data files have to be accessed, one can find in fact easy workarounds such as extracting the relevant data and storing them in another form (.mat/.fig/...) to speed up manipulation in Matlab. This question can therefore be closed. Best regards Marc
How to set perturbation in ode_PD2po
Hi Marc, If you find yourself adjusting the minimum stepsize to 10^-16 in order to achieve convergence, then there is likely something amiss in your implementation (or the problem is highly unusual). Does convergence appear to be quadratic? Generally speaking, I would recommend against adjusting continuation or nonlinear solver settings until you are sure that the need is motivated by unusual aspects of your problem. I am inclined to suspect an error in your Jacobians. Harry
Hello I am following e.g. SN points as function of two parameters, for a quite stiff system. In some cases, even when setting a large maximum step size, and a minimum step size of 1E-16, the chosen step size gradually becomes smaller, so many points are calculated. 1. Sometimes, this works fine: you just have to wait long enough. 2. Sometimes, COCO gets stuck at the minimal step size of 1E-16, when the Newton algorithm is not converging (even having increased ItMX from 20 to say, 100). Stopping and...
Hi Marc, Let me begin by asking for the reason for such a large number of points (I assume this is from continuation along a one-dimensional manifold). Is this necessary or a result of setting the step size maximum to a small value? Best, Harry
coco_plot_bd many points
Hello NK, Sorry for the delayed response. I notice a couple of things when running your code: The periodic orbit found using ode15s includes a couple of variables that are of order 10^-7 and 10^-5 while all the other variables are of order 1. I still think that could cause numerical difficulties as the numerical routines use relative tolerances. The third component of your vector field equals 0. That seems to suggest that the periodic orbit that you compute is not isolated (since any value of y3...
Hello again Harry! I have normalized the ODE as you suggested and made coordinate transformations for certain variables (the resulting ODE is the best I can think of at the moment to meet your requirements, it can be found in attached new ODE_MMC_SIMP.m ). I performed bifurcation analysis with three different ranges of Kp_CC (see attached coco_bifurcation.m): TEST1: Kp_CC belongs to [0.05, 0.6]. In this case, the program detected a torus bifurcation at Kp_CC=0.1129. By observing the corresponding...
Hello NK. My recommendation is that you try to rebalance your variables so that they are of similar magnitude and vary on comparable scales. When I graph the solutions that COCO finds, the variables are of very different orders of magnitude and vary over very different orders of magnitude. This is not good practice when doing numerics. Since COCO uses various tolerances for identifying convergence and estimating errors, significant imbalance between variables renders the results insensitive to some...