Pseudo Arc-length continuation

[Nidish96]

I've started adding an implementation of pseudo-arc length continuation to the code. I've checked it to be working for the SDOF Duffing oscillator. You can have a look at the code in the branch "Nidish" in my fork of the project (https://github.com/Nidish96/mousai/tree/Nidish). I'll raise a pull request once I'm done qualifying the rest. I'm raising this issue here for any suggestions you may give me. I'm having the following things on my agenda:

• Implement Jacobian through finite-differences for both Frequency domain as well as Time domain user supplied routines
• Qualify the code for MDOF systems
• Switch to a better solver
• Create a time-marching-based Monodromy stability estimation method

[josephcslater]

Wow. Sorry I didn't see this earlier.

So: I'm trying to refactor the code. I'm not happy with what I have now. The whole reason being that I want to adjust it to have a wrapper around the harmonic balance function. Right now I think that the version you forked at least outputs this- maybe only for the collocation harmonic balance method (time evaluation- hb_time), maybe the other. What would be better is a wrapper function that creates a generalized enough function with a single parameter for pseudo-arc length continuation. It's driving me nuts that my day-job is in the way of this.

Also, there are some stability issues to work on. Many of those are best worked on my simplifying the code (with this wrapper) and having the solver function content be mostly about the solution process, not the generalization of the build of the functions.

Further, I don't know how this happened, but the form of the function being used is not compatible with numerical integration in scipy. That's a big oversight on my part. This is the best way to encourage usage of this module and also provides a very simple way to validate solutions through numerical simulation. I'm very concerned about moving forward without fixing these points as it will make reworking that much more challenging. What are your thoughts?

[josephcslater]

I've finished a major bug fix (released as 0.3.1). I'm making a fork for 0.4 for refactoring per 0.3.1 release notes.

[Nidish96]

Hello, sorry for the late response.

I understand your concerns about the design, but I'm not really sure if making it compatible with numerical integration will offer any significant value addition. In most cases people are not interested in solving a system with a single frequency and looking at the time response but at the "frequency response" consisting of a sweep of frequencies. So it might even just be reasonable to leave it as is. I agree that adding a compatibility might encourage more users though.

I'll have a look at the changes in the newest release and see how we can incorporate continuation. In my fork I had pretty much written a wrapper function called "hb_res" which allows the user to specify the system either in time domain or in frequency domain itself, from which it constructs the frequency domain residuals and jacobians (also gives the option of inputting user-defined analytical jacobian or use an inbuilt finite difference scheme). Have a look when you're free and we can discuss in some more detail about tying both of these sets of ideas together.

[josephcslater]

None of the changes in the new release should be incompatible with what you've done. It was all a correction of the frequency method (Galerkin frequency) that was incorrectly finding (actually- portraying) higher harmonics. So, you should be able to update your code with a pull without problems. I'd prefer that you did before I look so that I can more easily see "real" differences.

That said- there are many places where numerical integration is indeed the current method of solution because of the lack of HB methods. So, attracting those people really requires commonality. It's also useful for validating a response. It shouldn't be much of a change to the code- two lines each maybe.

The other issue is that I was not aware of Python's ability to handle functions as first-class objects when I wrote this. It is much cleaner to pull out the error function into a separate external function and simply call it from the solver. This will provide the user direct access to it. I've found that to be of interest for what I'm doing as well.

I am branching all changes that may break your code for the 0.4.0 release. This way you can have some confidence moving forward and so can I. We can manually merge later.

[josephcslater]

So it seems that the new integrators to allow passing of parameters. I'm not a fan of that as design requires parameters that can be changed. (do you see otherwise?) One can use globals- but that's truly bad form. I'm trying to step back and make sure the next step is the right one. I am sure that I want the same function to be usable for numerical integration. This way Mousai can be used trivially for appropriate scenarios. I see

1. Use new integrators, commit to using globals (I would hate this...)
2. Use the old integrator- which is likely better/faster as it's fortran-based. This would not require globals.
3. Some sort of wrapper that can wrap functions from one to the other so that mousai can handle either. I have no idea how to do this.

I do think it might actually be best to logically separate out the pseudo-arc length continuation into a separate package as it has applications broader than just this. It can start here, but I think it's valid to stand it up on its own when it's ready, just like PyCont (which is frankly too dang complicated to work for real problems).

What are your thoughts?

[josephcslater]

I have just pushed to the 0.40 branch. It includes a new function that will wrap scipy integration functions (the functions called by the integrators) to the "new" Mousai form. For consistency, I still want to change the Mousai form, but for backward compatibility, I think (pretty confident) that I can detect the old form and convert it as well. The new format is "documented" best in _function_tomousai, which is intended to be private (thus the ''.). I'll be working on an "old mousai to new mousai" function, and don't know whether that should be private or not. I think long term it causes more confusion making it public than just switching, so I might make it private with a warning to adjust to the new format.