I am using the IMSL DIVPAG routine for numerical integration of an initial value problem (Gear option for stiff equations). I've hit a bit of a snag. I want to integrate for, say, 100 steps, and save the final system state vector. Then run the code again, read the system state vector from the first run to use as initial conditions,and integrate another 100 steps. I've done this and it seems to work, but results do not match a run of 200 steps in one go. The first 100 match, but not the second 100 (close, but not close enough). I'm guessing DIVPAG uses internally stored values that make the 1st step of the second 100 step run get calculated differently than the 101st step of the 200 step run. Unless I missed it, the IMSL documentation says nothing about how to smoothly restart an integration with a saved system state vector. Has this issue ever come up before?
Which version of IMSL did you use?
Please read these comments about IVPAG/DIVPAG by Dr. Richard Hanson, former chief scientist of IMSL/VNI/Roguewave:
Thanks for the reply (on the weekend, too!). I had read that same forum post, and looked at the DAESL documentation. I am presently using IMSL 6. In April last year I bought Parallel Studio XE for Fortran and C++, with IMSL 7. But at that time, after reading gotchas in the installation notes, I decided to stick with my current setup. Anyhow, that makes switching from IVPAG to DAESL a really big job. I would do it if DEASL looked like it would solve my integration start-stop-restart problem. But the two routines have a lot in common, including the lack of a documented way to do start-stop-restart. At least that's how it appears to me.
I tried hand-feeding IVPAG the only values in its PARAM vector offering influence over the next integration step, but that didn't help (i.e. 31,32,33 holding time steps).
I wonder if saving all the workspace data from I2PAG would allow a smooth start-stop-restart. But that's a lot of extra work, and I2PAG requires YTEMP — Array of size NMETH, and I don't know what NMETH is.
I see nothing about workspace in DAESL documentation other than IDO=3 to release it.
As it stands, I get small jolts artificially injected into the system by IVPAG on a restart, and I expect DAESL will also do this..
The annoyances that you describe regarding a restart are common to many ODE integrators with step-size control. Often, the integrator finds the solution at a value tj just beyond the user's requested tend, and uses some interpolation scheme over the interval (tj-1, tj).
I think that you should try the Intel ODE integrators in https://software.intel.com/en-us/articles/intel-ordinary-differential-equations-solver-library . They are of high quality, and can be called from Fortran and C.
Other stiff ODE integrators to consider:
The advantage of using the source code of an integrator is that you can intercept the solver at the largest tj that is short of tend, and use that as a restart point.
I'm lousy at reading legalese, but Intel's license says "to make only the minimum number of copies of the Materials reasonably necessary for your internal testing and development of your products". That leads me to think I should not use it. But just in case, IVPAG is being used without an A matrix or RHSJ jacobian function. Can you say if dodesol_mk52lfn is the right replacement? It has exactly one double and one integer workspace arrays. Would saving these allow for a smooth restart?
Would using these intel integrators require I install the version of IVF I got on April 2017?
I understand what you say about extended end points and interpolation. So it makes sense that data in addition to the TEND state vector would be needed for a good restart. Perhaps integrators for initial value problems should offer the option deliver a TEND state vector which allows for a smooth restart if the user wants that.
Can you provide a description of the ODE system that you are integrating, along with time-interval, parameter ranges of interest and initial conditions? In other words, the mathematical problem statement?
I am certainly unqualified to give legal advice, but it seems to me that the injunctions against copying, etc., would come into play only if you build an application that uses the Intel ODE solver and then try to distribute it. The users of your software would need to dowload the ODE solver DLLs. In fact, the legal issues are more complex for redistributing applications that use IMSL-7. See
Why would there be any problems during the development work that you are doing?
It is a calculation of response of a structural dynamic system modeled with Finite Elements. The specific sub-field is called rotordynamics. It's all about vibration of rotating machinery ranging from very small very high speed like automotive turbochargers, to very large low speed rotors like hydroelectric power generation turbines and generators. Initial conditions are often just zero, but not always. Simulation times might range from a few 10ths of a second to 20 or 30 seconds, give or take. A few thousand output points is often adequate. Suitable time step sizes are often tied to the rotational period of the rotor, typically around 1/20th of that if the model is linear, and smaller if nonlinear. Does that help?
There's no inherent problem with doing development work, except there aren't enough hours in the day. I was looking for a quick fix to the restart issue with the IMSL IVPAG integrator. A quick fix is all I have time for, anything more goes on a back burner along with a lot of other things.
I will ask my sponsor to look at that web site about the IMSL license. Q7 seems to apply. He won't do that himself, but he should have access to someone qualified to do so. My sponsor is the Turbomachinery Laboratory of Texas A&M University. Texas A&M may already have a license worked out with Rogue Wave. The end users are graduate students and industrial members of the Turbo Lab's industrial research consortium. So users are both concentrated at the university, and scattered around the world.
Today I just heard from a user in China who says the code is not running. That's on the front burner. He's doing something other than time integration.