Detecting many apparently non-existent BP and FP
Ok, thanks a lot!
Hi Jonathan, The phase lives on S^1, but since Matlab works in R^1, I map S^1 onto [-pi,pi] and identify the end points. The jump from pi to -pi just brings the phase back to the beginning of the interval. There is no relationship to directions. The phase always increases at a constant rate. I hope that helps. Harry
Ok, now I understand the code, thank you. But regarding the meaning of this: Why should there be a phase jump from pi to -pi? Is it because of the change of direction of the system? Because if yes, I think I don't need that event function in my system.
Thank you for your question. The po toolbox constructors build problems associated with periodic trajectory segments, i.e., those for which the two end points coincide. The hspo toolbox constructors can handle other types of relationships between the end points of successive segments, including jumps. In the non-autonomous encoding of the linear oscillator, the two state variables are identical at the two end points, and so po works fine. But if you augment the state to include a phase variable,...
Hello Harry, In order to better understand the functionality of the 'po' toolbox, I attempted to reproduce the results of the linode example by encoding the system as an autonomous set of equations. However, when I run the continuation problem with variations in the stiffness of the oscillator, I don't find a solution manifold as the COCO output under type shows MX. One step I took towards troubleshooting this issue was verifying that the Initial trajectory from the autonomous system that I was feeding...
Hi Jonathan. The phase reset event triggers when x(3)-pi=0 and at that point resets x(3) to x(3)-2*pi, in other words so that while the phase is pi at the terminal point of the incoming segment, it is -pi at the initial point of the outgoing segment. Similarly, the impact event triggers when x(1)-p(5)=0 and at that point resets x(2) to -p(6)*x(2), in other words so that while the displacement is p(5) at the terminal point of the incoming segment, the outgoing velocity is -p(6)*x(2). Here, p(6) is...
Hi Jonathan. The phrase reset event triggers when x(3)-pi=0 and at that point resets x(3) to x(3)-2*pi, in other words so that while the phase is pi at the terminal point of the incoming segment, it is -pi at the initial point of the outgoing segment. Similarly, the impact event triggers when x(1)-p(5)=0 and at that point resets x(2) to -p(6)*x(2), in other words so that while the displacement is p(5) at the terminal point of the incoming segments, the outgoing velocity is -p(6)*x(2). Here, p(6)...
Thank you, Harry! I think I get the concept of the impact demo. But the part with the phase reset event is not clear to me. I understand it like if omega * t = pi then the reset map changes it to omega * t = 2 * pi But I still search for the reason for that.
Oh, and since the hspo constructors assume an autonomous vector field, you can either replace your harmonic forcing with a scaled version of one of the components of the Hopf normal form dynamics or introduce a phase variable that gets reset at a phase reset event (like in the impact demo.
Hi Jonathan, The hspo constructors were designed to handle problems with piecewise-defined vector fields, so that discontinuities would not affect the accuracy of the results. There are two issues here: the lack of smoothness (or even continuity) of the vector field and the lack of smoothness of the solution. The former means that your problem Jacobian will change discontinuously as mesh points move across the discontinuity. It is still possible (perhaps even likely, but surely problem-dependent)...
Error when following TR
fixes
Hello Harry This works perfectly, thank you!! Marc
Harmonic oscillator with piecewise excitation force
Thank you, Marc. I very much appreciate you bringing this to my attention. You have attempted to do something very logical, namely restart continuation of torus bifurcations from a point on a previously computed branch of torus bifurcations. As you note, the same functionality appears to work for period-doubling bifurcations, and it seems very reasonable that po constructors should support this. Indeed, I see no reason why they shouldn't. Obviously, I did not test this functionality or else I would...
Hello, thank you for your answer. I could not find a coding error so far. However, I was able to reproduce the error by modifing the file demo.m in the tor folder. 1. PD %% Continuation of period-doubling bifurcations PDlabs = coco_bd_labs(bd2, 'PD'); prob = coco_prob(); prob = coco_set(prob, 'cont', 'NAdapt', 5); bd4 = coco(prob, 'pd1', 'ode', 'PD', 'PD', ... 'po2', PDlabs(1), {'nu' 'po.period' 'be'}, [-0.65, -0.55]); %addendum prob2 = coco_prob(); prob2 = coco_set(prob2, 'cont', 'NAdapt', 5, 'ItMX',...
Hi Marc, The tor demo in the po toolbox folder shows the continuation of families of saddle-node, period-doubling, and torus bifurcations of periodic orbits. The script syntax you have shown above seems correct, so I am inclined to assume a coding error. Feel free to share your full set of files for debugging, or consider comparing your code against the demo (which I hope runs without difficulty for you). Best, Harry
Error when following TR
Hi Harry, that is helpful, thank you! Best, Jonathan
Hi Jonathan, Your question is both general to numerical methods and specific to COCO. For the general answer, I refer you to a multitude of other sources that explore the computational cost and benefits of explicit derivatives. For COCO, and given the use of derivatives in detecting special points (e.g., bifurcations), my general practice is to encode derivatives of vector fields, at least to first order and in many cases to higher order, including directional derivatives, as appropriate. The symcoco...
Thanks a lot, I think there was a problem with the derivatives! Maybe a general question regarding that topic: is it sufficient to just let Coco calculate the derivatives numerically (easiest way to implement in my eyes) or is it better to always give analytical derivatives if possible? Best, Jonathan
Dear Harry, many thanks for looking into our code and your advice, we will try to adapt our code accordingly. Kind regards, Joseph
Hello, Jonathan. This is difficult to debug in the absence of your actual code. My gut guess is that you have made errors in the encoding of derivatives. That's where I usually start, rather than to assume that the error is due to something inherently remarkable about the problem at hand. In case you are able to make use of the symcoco routines for automatic generation of vector fields and their derivatives (see symcoco-doc.pdf), consider comparing the results. Best, Harry
Hello, Jonathan. This is difficult to debug in the absence of your actual code. My gut guess is that you have made errors in the encoding of derivatives. That's where I usually start, rather than to assume that the error is due to something inherently remarkable about the problem at hand. In case you are able to make use of the symcoco routines for automatic generation of vector fields and their derivatives (see symcoco-doc.pdf, consider comparing the results. Best, Harry
Dear Joseph, I see that your problem involves use of the ode_isol2hspo constructor. To construct the corresponding adjoint contributions, use the adjt_isol2hspo constructor, as described in Section 12.3 of the PO_Tutorial.pdf file. For the monitor function, I note that it uses the trapz function to perform numerical integration using the trapezoid rule. Since optimization will require derivatives of your monitor function, you will need to provide these explicitly, as first-order numerical differentiation...
Hi Harry, it took a while for me to have time to go back to that topic... First, thanks for your last tip, I managed to implement those. Nevertheless, for some "more complicated" parameters in my system I get the following warning by Matlab during continuation: Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.252075e-17. I'm not sure how to deal with this warning. I have the feeling that the results are also not correct. Do you have any tips for me? Thanks...
Thanks Harry, your advice in this matter is greatly appreciated.
Hi Joseph, This kind of functionality is made available in COCO using the staged construction of adjoint conditions. You can find several demos in the release that rely on this functionality., e.g., duffing_optim. The methodology is described in several publications, e.g., "Staged Construction of Adjoints for Constrained Optimization of Integro-Differential Boundary-Value Problems" and "Methods of continuation and their implementation in the COCO software platform with application to delay differential...
Detection and continuation of extrema of user-defined monitor functions
When lambda equals 0, the two equilibrium branches are x=0 and mu=x^2, which intersect at x=mu=0 at a branch point (not an SN). This is a singular point of the continuation problem. When you include explicit derivatives, points at some distance away from the singular point can be located using Newton iterations within the allowed number of iterates and the continuation algorithm typically jumps across the singular point along either of the two branches. When you omit explicit derivatives, convergence...
Hi Harry, Thank you very much for the quick reply and for pointing out the error in encoding the Jacobian. That was silly of me... Everything is, of course, now working as expected. Thank you also for pointing out the tutorial. If I can ask a follow up question: if I pass to ode_isol2ep() only a handle to the function, and do not also provide the Jacobian, coco() terminates with MX near the SN at lambda=0. Is this to be expected? This is with no changes to default options. Many thanks again, Dav...
Issues continuing simple two-parameter ODE
Thank you, David, for reaching out. First of all, I recommend that you have a look at the EP_Tutorial.pdf documentation, which includes a demo that concerns precisely the dynamical system you are investigating. Secondly, it appears to me that there may be an error in the encoding of the Jacobian with respect to the problem parameters. Finally, I recommend against using coco_set to change default options, especially of the corrector, as more than likely unexpected behavior is not due to their settings...
Issues continuing simple two-parameter ODE
If you are looking to perform bifurcation analysis of periodic orbits in piecewise-smooth or hybrid dynamical systems, please see the tutorial documentation for the 'po' toolbox. If you are looking to perform bifurcation analysis of equilibria in piecewise-smooth systems, with emphasis on border collision bifurcations, you will need to extend the 'ep' toolbox with your preferred monitor functions and events, per the methodological guidance in Recipes for Continuation. I hope that helps.
Please see the responses in ticket #51. Kind regards, Harry
Hello ，How to use COCO for bifurcation analysis in non-smooth systems
Sorry, I missed to attach the input file in the first post. Here it comes. With kind wishes Alois
Continuation of a P2-orbit in an non-autonomous ODE: How to deal with the gluing conditions
Dear Harry I truly appreciate your responses. I got it working! Best, Ali
The function I had in mind you using was po_add_func and I suspect that's what you, in fact, did. This function is a wrapper to coco_add_func that also includes a preceding call to coco_get_func_data. Together these commands extract required information from the associated 'coll' instance and then impose a further constraint on the corresponding continuation variables. The reason for introducing po_add_func as part of the GettingStartedwithCOCO tutorial was to avoid the need to require i) explicit...
Dear Harry, Thank you for your fast response I looked at section 7.3 and 7.4 as you said. I saw that it uses "po_add_bddat". I tried it with my problem and I managed to monitor the orbits. Thank you for your guidance. However, I have used coco for some time now and I have always used "coco_add_func" for all my monitor and zero (constraints) functions. Since I know how to apply constraints using "coco_add_func", I was wondering if there is a way to feed "xbp" into "coco_add_func" (type zero) and define...
Dear Ali, Thank you for reaching out and for contributing to this thread. I am wondering if the information in the Getting StartedwithCOCO tutorial document and demos (new with the August release of COCO) may be helpful to you. For example, Section 7.3 shows how you can monitor properties of periodic orbits during continuation and Section 7.4 shows how to track and constraint orbit maxima. There are other ways to do this as well, but perhaps this material is a good start. Best, Harry
Hi Harry I am using PO toolbox to calculate the frequency response of a vibration system with nonlinear stiffness. I have one problem that I do not know how to solve. In the first step, I calculate the frequency response (just like the example for Duffing in the tutorials). However, I do not know how to access the time response (orbits) in each step. I need them because I need to apply constraints on the amplitude and/or the phase of the time responses while relaxing another parameter in the system....
As described in Section 2 of ATLAS-Tutorial.pdf, you can see the settings associated with the 'atlas_1d' atlas algorithm by typing atlas_1d_settings on the command line. When you write 'cont' in the call to coco_set, you are specifying that the following name-value pair defines a setting associated with the atlas algorithm. In your example, you are setting the 'NAdapt' setting to 1, implying that adaptive updates will be made to the zero problem before each new continuation step. Individual toolboxes...
Those coll_settings etc. commands are very useful. This solves this issue, thank you!
As described in Section 2 of ATLAS-Tutorial.pdf, you can see the settings associated with the 'atlas_1d' atlas algorithm by typing atlas_1d_settings on the command line. When you write 'cont' in the call to coco_set, you are specifying that the following name-value pair defines a setting associated with the atlas algorithm. In your example, you are setting the 'Nadapt' setting to 1, implying that adaptive updates will be made to the zero problem before each new continuation step. Individual toolboxes...
Hello, Marc. Continuation may terminate because of a lack of convergence of the nonlinear solver or because of the detection of a user- or toolbox-defined terminal event. In the former case, the label is defined by the atlas algorithm. In the latter case, the label used to describe the event is also user- or toolbox-defined. For example, the MXCL label describes a terminal event that is triggered when the 'coll' toolbox estimates a truncation error of the discretization of a trajectory segment that...
Additional parameters : prob = coco_set(prob, ’cont’, ’NAdapt’, 1); adaptive setting of the Atlas algorithm (see [PO-Tutorial.pdf, p. 5])
Different continuation fails and what to do
required function not found' with '-var' option in ode_isol2po
Thank you, Jacob, also for the greetings! Although option #2 provides a fix, the behavior you have identified should count as a bug, since a construction error is thrown for a syntactically correct construction. I will aim to encode a transparent and well-documented solution to this in a future COCO update. Kind regards, Harry
In case anyone from the future sees this, I've attached the example code with the fixes implemented.
Kia ora Harry, Thank you so much for the response. The second option has worked perfectly. On an unrelated note, Bernd and Hinke say hi! Cheers, Jacob
Kia ora Harry, Thank you so much for the response. The second option has worked perfectly. On an unrelated note, Bernd and Hinke say hi! Cheers, Jacob
Thank you, Jacob, for this question. I am very happy to see that you are exploring COCO and attempting to reproduce existing results using independently developed code. There are several parts to this that I will address below. 1. Even if the error you encountered hadn't been thrown, your code would not have executed, since the introduction of the variational problem added 9 to the dimensional deficit. To prevent this from happening, you would need to also introduce 9 constraints, e.g., using the...
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...
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...
required function not found' with '-var' option in ode_isol2po
Amplitude definitions
Hopf bifurcation continuation with 1 parameter only
Thank you!
Hi Jonathan, This is Matlab syntax for an anonymous function of three arguments (p, d, and u) that returns two arguments, one of which evaluates to the input argument d and one of which evaluates to the expression u(1)-2*pi/u(2). Look at the Matlab manual for more information about anonymous functions. The reason to use this particular syntax here is the requirement that zero and monitor functions added using coco_add_func must be of this form. An alternative is to provide a function handle to a...
Hey Harry, my further question is how to understand the part @(p,d,u) deal(d, u(1)-2*pi/u(2)). Which function is meant with @(p,d,u)? I get that we somehow have to define u(1)-2*pi/u(2)=0.
Dear Harry,thank you very much for your kind response and the fast treatment of the problem, I am really impressed.I guess, that it is actually a fault in Matlab on Windows, because the help entry for strsplit in Matlab states that "The substrings specified in delimiter do not appear in the output newStr." so I think I should send a bug report.I am learning coco and could already find a family of periodic solutions containg BPCs and LPCs. I still have to get familiar with the details, but coco makes...
Hello Harry, thanks for this extremly detailed and fast answer! I'll give my best to get along with my problem applying your tips. Have a nice weekend, Jonathan
coco_close_efunc: Added text to warning message
One further observation: however you choose to define your vector field, you can always perform the conversion between excitation frequency and excitation period in your post-continuation analysis, e.g., when graphing the frequency-response diagram using coco_plot_bd or your own plotting routines.
A longer answer is as follows: Since COCO is an object-oriented development platform, the software is agnostic to the way you choose to construct your zero functions and monitor functions. The COCO toolboxes package a limited set of problem constructors for problems that are common in basic bifurcation analysis. These are each wrapped around one or several calls to the COCO core constructor coco_add_func. As you work with particular classes of problems, you may realize that you repeat the same actions...
Frequency Response Plots using COLL/PO toolbox
Hello, Jonathan. Thank you for your interest in using COCO. The short answer to your question is as follows: Preferred toolbox For the continuation of periodic orbits in a smooth dynamical system, the 'po' toolbox is preferred to the 'coll' toolbox, since the former makes use of the latter and then appends additional constraints imposing periodicity and, as appropriate, non-degeneracy using a phase condition. The demo in po/examples/bistable shows how to use the 'po' toolbox for continuation of periodic...
symcoco: syntax error in Matlab 2023b
Frequency Response Plots using COLL/PO toolbox
Thank you Alois for flagging this and Jan for fixing it! An updated release (coco_2023October26) has been posted with this update.
fixes
Updated sco_sym2funcs to ensure proper function endings for both Unix and Windows-based Matlab systems.
You are guessing right, I use Matlab in Windows. I just tried your suggestion with strncmp and that also works. Thanks again for your efforts! Alois
Mhm, it seems that the behaviour is operating-system dependent? What operating system are you using? (I had tried linux and Matlab online, can't check it on windows). In Matlab online the behaviour is as I expected: K>> strsplit(sco_symcode({xy},x,y),newline())' ans = 9x1 cell array {0x0 char } {'function out1 = sys_0(in1,in2)' } {'%SYS_0' } {'% OUT1 = SYS_0(IN1,IN2)' } {'% This function was generated by the Symbolic Math Toolbox version 23.2.'} {'% 23-Oct-2023 15:33:38' } {'out1 = in1.in2;' } {'end'...
Dear Jan, thank you very much for your patient treatment of the problem! I put a breakpoint at the end of funend and could observe, that the treatment of the newline character in string functions might be the problem: After splitting the string at the newlines, the newline character remains in the strings, so one gets 'end\n' (actually Matlab 2023b distinguishes between the placeholder '\n' and the actual newline character <nl>, so the second last entry in str is 'end<nl>'. Therefore strcmp gives...
That's strange. I downloaded from sourceforge, applied your modelprep and got a valid file. Double-check that, if you open sco_sym2funcs, at the bottom of the file is a test function, looking like below (it checks if "end" is appended). Make sure that you have the new sco_sym2funcs in your path. %% test if output creates fucntions with 'end' keyword function funend=sco_test_for_end() x=sym('x'); y=sym('y'); if sco_isoctave() newline=@()sprintf('\n'); %#ok<sprintfn> else newline=builtin('newline');...
Dear Jan, thank you very much for your quick and helpful response! I downloaded and installed the new coco version on August 19; today I downloaded it again and the zip-files are identical. Still I get the error message with both Matlab 2023a and 2023b. I will try to find out what happens during the day. With kind regards Alois
Dear Alois, This is a known problem with the old version (Mar2020) of symcoco [1]. I just downloaded your example files and the problem no longer occurs for symcoco attached to the the August 2023 version. Attached is the (renamed) output with the new version, which behaves (hopefully) correctly. Best, Jan [1] Unfortunately, from one Matlab version to another the symbolic toolbox function generator matlabFunction changed its output format, appending the "end" keyword. The new version performs a test...
symcoco: syntax error in Matlab 2023b
Dear Marc, Thank you for your query. I'm happy to hear that you are giving COCO a try. Learning by modifying demo code is a great way to explore a package like COCO. The short answer to your question (but see below) is as follows: In the demo, I first detect a Hopf bifurcation under variations in one parameter ('B') and then continue along a curve of such bifurcations under simultaneous variations in two parameters ('A' and 'B'). If you don't allow 'A' to also vary, then you will at best detect a...
Hopf bifurcation continuation with 1 parameter only
fixes
The proposed change in the nullspace function is included as commented code as of today. Speed test suggest minimal difference for the po/examples/tor demo.
This took many years, but a recent change to atlas_1d scales the step size by norm(x,inf) which appears to produce consistent behavior.
Potential bug in coco_remesh
Add function for adding integral constraints to COLL.
See po/examples/int_optim
Enable adaptation by default.
Due to backward compatibility, mesh adaptation is disabled by default and must be enabled by the use of the coco_set command. Error estimation for trajectory segments and termination when the error estimate exceeds a threshold is enabled by default. The latter may be disabled by the use of the coco_set command. As a consequence, continuation of solutions to boundary-value problems, e.g., periodic orbits from a Hopf bifurcation, may result in premature termination at an MXCL point, corresponding to...
Classification of HB points
The 'ep' tutorial documentation and demos includes an example (chemosc) that accomplishes this for a quadratic vector field. The GettingStartedwithCOCO tutorial documentation and demos shows a complete encoding for an arbitrary vector field.
Check if PD->PO restarter adapts mesh.