I am running a source code using the Intel® Parallel Studio XE Composer Edition for Fortran with Rogue Wave IMSL 7.0 for Windows*
I am using a solver called the IVPAG and it returns an error-"Maximum number of iterations exceeded"
I have investigated the IMSL manual and documentation to obtain the parameter that will allow me to increase the number of steps, but there is no such parameter available to edit. Several other solvers allow to modify this parameter , but I do not see any option to change this parameter in IVPAG.
How can I overcome this problem?
You may find it worthwhile to use DAESL instead of IVPAG; see the comments in #2 of https://forums.roguewave.com/showthread.php?1542-Ivpag .
A pertinent question: are the ODEs that you are integrating stiff?
Hello. Thank you for your response. Yes, they are stiff. However, I have had good success with using IVPAG in terms of the quality of the solution with the other input files. So, ideally I would have liked to continue using IVPAG if that is possible.
Typically, when I have come across this error, I was able to solve it using a trial and error combination of step size and error criteria. However, it doesn't work this time and I had to see if there is any other parameter that I can change toward fixing this error? For instance, I've seen other solvers that were able to set the parameter for the
I could try using the new solver you suggested and see how it compares to the previous results. I will keep things posted in the forum. Thank you for your suggestions.
What do you hope to find out from the user supplied function? The only reason I am asking this is because as you may have already guessed, there are several lines of code and some of the parameters have been passed on from another subroutine where the IVPAG is called. If I know why you may want to look at the user supplied function, I can figure out which sections of the code to show here?
It is characteristic of many algorithms that they work well on most problems, but fail on a small set of pathological problems. In such cases, a different solver may succeed with no special effort required. All that we are interested in is to obtain from you the exact specification (ODE, initial conditions, range of independent variable) of one pathological problem instance.
When the ODE contains some parameters, we need just the values of the parameters, even if some long calculation had to be done in order to obtain those values.
The first call is made to DIVPAG as follows inside a DO and CONTINUE loop.
DO 1000 K = 1,NSTEPS
TOUT=TOUT+SSIZE ------SSIZE is the inputted step size, in this case 0.4
TIME(K)=TOUT -----------TOUT=0 before it enters the DO loop
CALL DIVPAG (IDO,N,DIFFUN,FCNJ,A,T0,TOUT,EPS,PARAM,YO)
These parameters are:
A=A(1,1) //Not very important I believe. It uses default value.
PARAM is defined as
PARAM(4) = 500000
PARAM(6) = 5
PARAM(10)=MABSE //This value is 1
PARAM(12)=MBDF //This value is 2
PARAM(13)=MSOLVE //This value is 2
PARAM(19) = 0
User supplied function=DIFFUN. DIFFUN is a separate subroutine and has parameters as
YO is the matrix that contains the initial values for the ODEs and the new values at t+delta t are stored in YDOT. Then, it continues until convergence.
I do want to mention that in some of the cases, the program never seems to come back out of the solver. Normally, I would expect that it either gives the solution or "maximum number of iterations exceeded" error. However, there are times when neither of these show up. And it seems stuck in the solver. The reason I know this is because I have a print statement after the program comes out of the solver. I wonder why this happens.
There is something seriously wrong in the way that you use IVPAG, or I am misreading what you wrote in #7.
Please study the example ODE problems and codes in the IMSL documentation and construct your program in a similar manner. Many IMSL subroutines have rather complex interfaces, and care is needed to use them correctly.
In hindsight, had a means been provided in IVPAG to increase the number of iterations, you could have found yourself in a worse situation.
It would be useful if you gave a statement of the mathematical problem at this juncture. That would set the stage for how to proceed with obtaining a numerical solution.
I will provide two separate comments. In this comment, I will provide responses to the issues you have mentioned in your comment. In the second comment, I will describe the problem we are trying to solve.
The problem we are trying to solve is as follows-
We have a chemical system. This chemical system has different species namely C1,C2,C3, and C4. For now, we will use 4 species but we end up having around 150 species in our system. They interact with one another. What I mean by this is simply that,
C1+C2 -> C3 k1=rate constant
C2+C4 -> C1 k2=rate constant and so on. Therefore we take into account their interaction which changes the nature of the system. So, if we started with 10 of C1, after a certain time, t, we would have 6 of C1.
The goal is to quantify the values of C1,C2,C3 and C4 after a "time" desired by the user. We are able to describe their interaction by a parameter called "rate of reaction" given by dC1/dt = -k1*C1*C2 + k2*C2*C4. So that, here C1,C2,C3 and C4 are independent variables and t is the dependent variable.
Similarly, we write such simultaneous ODEs and solve for each of these variables. They are stiff likely because their values are in a wide range (I believe this is the reason for their stiffness but it could be anything else). For instance, C1 could be 0.001 and C3 could be 10^-11.
It is to solve these simultaneous ODEs that we use DIVPAG.
Ah, I see.
In #7 you wrote:
YO is the matrix that contains the initial values for the ODEs and the new values at t+delta t are stored in YDOT.
That statement is what caused me alarm. Likewise, I saw no updating of T0 in the code fragment in #7.
Simultaneous chemical reactions are a well-known source of stiff ODEs. Robertson's example is famous. For comprehensive coverage of algorithms, software and comparison of results, see Francesca Mazzia's Web page at https://archimede.dm.uniba.it/~testset/testsetivpsolvers/ and Hairer's Web page at https://www.unige.ch/~hairer/testset/testset.html .
Since the characteristic response times of the simultaneous reactions can span many orders of magnitude, using a uniform time step is not a good choice. See, for example, the time axes in the figures at https://archimede.dm.uniba.it/~testset/problems/rober.php .
Thankyou for the useful links.
I am not sure how to choose the right step size at every iteration, especially when solving computational problems. While solving it by hand, one can make a guess to achieve the solution because we know the values at each iteration. Of course, these were really small problems and therefore, less iterations, so it's feasible.
However, I am not sure how I can do this in a computer program, especially, if it has 1000+ iterations, so it is obviously not practical to look at the value after every iteration and provide a new step size. A uniform step size mostly works or has worked so far for all our systems. Any thoughts about this?
And again, are we able to go back and re-think on modifying the parameter that can change the maximum number of iterations?
>>However, I am not sure how I can do this in a computer program, especially, if it has 1000+ iterations, so it is obviously not practical to look at the value after every iteration and provide a new step size.
Assume you have a candidate guess at useable step size. You can partition the guess into 2 or 4 sizes that bound the guess and run those in parallel. Upon return of all partitions determine a new guess.
Please note, if the number of cores/HT's are underutilized when running without partitioning the step size guess, then the extra work run with partitioning may come mostly free. And, by improving the step size estimations, you can reduce the total number of steps, thus potentially reducing runtime.
Hi Divya K.,
You can change the number of steps at the array "param(4)". The standard value is 500. I usually put something very high, like 2d9. I believe this is the limit value. Also, another thing you should check it is the tolerance. If I were you, I would print the time variable in the FCN (DIFFUN in your case) subroutine in order to have an idea about the step size.
Hi Jim Dempsey,
Related to you last suggestion, it is very interesting. But, I'm not following you. :-( How can I do the partition in the guess step?