Hallo Ulrich! How does method behave for stiff systems? Danke!

Hallo Eduardo. Good question! Stiff systems require small timesteps, adjusted to the 'stiffest of the degrees of freedom'. Since the algorithm is simple and fast, small time steps are acceptable. There is a Wolfram demonstration contributed by me that allows to experiment with the method and which also has a stiff example.

@Andrés. Thanks a lot for your comment which led me to some interesting public sources describing Gragg's work (which was not known to me). Actually we (i.e. Gragg's community on one side and I on the other side) seem to live in different worlds. My main interest is in algoritms that are straightforward and streamlined enough that it is not too unnatural a idea for me to think that Nature follows a similar method in letting systems evolve in time. This excludes implicit methods and, very probably, also sophisms like higher order Runge Kutta and those you mention. Experience taught me that when a wide range of physics applications is under consideration then sophisticated methods work well for one example and very badly for another. Nature can't work that way. For your information, I added a link to the promotional description of my method. Whenever you come accross something similar, please let me know again.

http://www.ulrichmutze.de/cpmgraphics/graphics1209.png

@Eduardo.

Let me come back to behavior for stiff problems. I agree with the emphasis that you put onto this issue. Early this year I discovered that the method performs extremely well when coupled to a very simple step control method. See the link for details.

http://demonstrations.wolfram.com/TheAsynchronousLeapfrogMethodAsAStiffODESolver/

Claude Massé wrote in the thread

post/Which_is_the_most_interesting_and_widely_used_data_structure?ev=home_person_add_comment_question_target

which, as far it refers to my integrator proposal, I would like to re-direct into the present thread:

H.E. Lehtihet wrote in the thread

post/Which_is_the_most_interesting_and_widely_used_data_structure?ev=home_person_add_comment_question_target

which, as far as it refers to my integrator proposal, I would like to re-direct into the present thread:

Halim (is this OK for you?),

a partial answer to your question concerning performance for stiff problems is given in the Wolfram Demonstration appended to my note to Eduardo. It would be great if you would try to activate the interactivity of this program (you had to download from Wolfram Research a free 'CDF-reader-addon' to your browser.

In the nice book by Paul Taylor: Practical Foundations of Mathematics , I found (on p. 470):

>

The Wolfram Demonstrations together with their 'computable document format' CDF seem to be presently the easiest way to realize Taylor's vision.

Before I prepared this Wolfram Demonstration, I worked out as a C++ program the example of a non-linear electric circuit (a ring modulator) which I found in a collection of stiff problem test cases in the internet. It reproduced the 15 extremely structure-rich current versus time diagrams (and value lists) with perfect accuracy. So I'm quite sure that stiff systems are not among the weak points of the method.

Ulrich, you haven't understood what is the order of a method. The method may yield the exact solution in certain cases, and still be of first order. It is accidental cancelation. The matter isn't who is right, I don't care, and then this discussion is of no value for me.

Claude, in my last argument I intentionally avoided to speak about order, and so I can't see a basis for your diagnosis about my lack of understanding of this concept. The discussion would have had a value for you if it would have led you to absorb the basic numerical math insight

( f(x + h) - f(x) )/h = f ' (x) + O(h)

( f(x + h) - f(x - h) )/(2h) = f ' (x) + O(h²)

which, by the way, is nearly evident geometrically.

(for the second of my formulas see for instance Abramowitz/Stegun, p. 883,25.3.21).

Ulrich, I think I can get the time this week-end to take a look at this method.

Perhaps you should edit your above post and add a space between f and the prime < ' >I symbol to make the relations more readable.

BTW, the first relation can be written also as follows :

(f(x+h) - f(x))/h = f ' (x+h/2) + O(h²)

This could be relevant to your method but I am not sure...

Halim

:

Halim,

this would be great!

I would be glad seeing unbiased people experimenting with 'Asynchronous Leapfrog'.

Ulrich: (Page 1/3)

Unfortunately, I could not get to download the application. Nonetheless, I have done some programming. This is definitely a second-order low-cost jewel, with however some potential problems as will be discussed later. But first, I would like to make sure that I have understood properly your method, so please correct me if I am wrong.

For simplicity, I write f(k) instead of f(t(k), y(k));

Initialization:

h = dt/2

y(0) = y0

S(0) = f(0) // initial slope

Iteration:

Do (k = 0 … N)

--- Step 1) y(k+ ½) = y(k) + h*S(k)

--- Step 2) S(k+ ½) = f( k+½ )

--- Step 3) S(k+1) = 2*S(k+ ½) – S(k)

--- Step 4) y(k+1) = y(k+ ½) + h*S(k+1)

End do

Step 1: Gives the approximate solution at mid-interval

Step 2: Gives the estimated slope at mid-interval via a call to f.

This is the only call to f per iteration.

Step 3: Gives the estimated slope at the end of the interval.

This is done via a linear extrapolation with no call to f

Step 4: Gives the approximate solution at the end of the interval

These steps have a nice straightforward visual interpretation (see figures). The top figure shows the time evolution of y and the bottom figure shows that of the slope. The colors correspond to the above steps. Their meanings are as follows:

Yellow = intermediate y, used only to estimate the slope via f.

Red = calculated slopes (including initial slope)

Blue = extrapolated slopes

Green = y outputs

Ulrich: (Page 2/3)

In the following

E2H = Explicit 2-order Heun's method.

E2M = Explicit 2-order Midpoint method;

E2LF = Explicit 2-order Leap-Frog (or modified midpoint) method;

First iteration (k = 0).

The first iteration is particular in the sense that the value calculated in Step 4 will always coincide with that of E2M. Indeed, Step 1 and step 2 are already the same as those of E2M. Step 3 is nothing more than a simple linear extrapolation (in the derivative space) to estimate the slope S(k+1) at the end of the interval. Finally, Step 4 (when combined with the previous step 1) can be rewritten either as:

y(k+1) = y(k) + h*( S(k) + S(k+1) )

or, if we account for Step 3 (extrapolation),

y(k+1) = y(k) + 2*h* S(k+ ½).

The first form (trapezoidal rule) is similar to that of E2H. Except that, here, S(k+1) would be slightly more accurate than that of E2H. This can be seen directly in the second form, which turns out to yield exactly the same result as that of E2M.

In fact, if we remove the instruction S(0) = f(0) from the initialization and if we insert instead a new Step 0 : S(k) = f(k) in the iterations, then we will end up with the code of E2M (albeit an awkward code for E2M, I must concede).

Subsequent iterations (k > 0).

For k > 0, the iterations will differ significantly from those of E2M precisely because Step 0 is missing. In other words, the fundamental difference between E2M and this jewel is that E2M does update the slope at the end of the interval whereas the jewel skips this update. Skipping the update not only offers an obvious advantage in terms of calculation cost but also, interestingly, does not downgrade the order of the method. The reason is that, for k > 0, the iterations use a second-order process very similar to that of E2LF. Is this a Free Lunch?

Ulrich: (Page 3/3)

Is this a Free Lunch? Well, there is a slight difference with E2LF.

First, the original E2LF process would be written as follows:

y(k+1) = y(k –1) + 2*dt*f(k)

whereas the jewel uses instead two interleaved E2LF-like processes:

one for the middle (yellow) points

y(k+ ½) = y(k – ½) + 2*h*S(k).

and another one for the green points

y(k+ 1) = y(k) + 2*h* f( k+½ ).

This is not directly visible in the algorithm but can be derived easily by examining for two successive iterations, respectively, Step 1 (yellow) and Step 4 (green). In any case, the above expressions yield a third-order local error and a second order global error (thus my remark of the other day, has turned out to be relevant). However, the S in one of these expressions is always a "blue" slope (extrapolated) whereas the original E2LF makes instead a call to f to get an updated slope. I don't know if this is a serious problem but imho it deserves to be carefully analyzed.

Second, perhaps more problematic is the fact that the E2LF process is known for its so-called "week instability". This is why this process always incorporates a smoothing operation for damping the oscillating error. This operation is done every N iterations (with N even, usually N about 20 or 30). In the case of the jewel, this operation is not done and this could be a serious issue that needs to be checked.

Halim,

you got the method perfectly right. The differences to your E2LF are discussed also in the appended articles (the original, and a further developed version on my homepage) where E2LF is taken as a starting point. Then a modification is made that allows one to change the time step dt from iteration to iteration. I recognize that this changes the method. In the test worked out there (and all the tests I did later) the performance of the method grew by the modification! In the second of the linked articles also including 'relaxation' is discussed and found not to have much influence. So there are indications that the method is more stable than E2LF. Of course, it would be desirable to make this more precise. It would be great if you could propose a suitable test equation.

Thanks a lot for your insightful analysis.

http://www.ma.utexas.edu/mp_arc/c/08/08-197.pdf

http://www.ulrichmutze.de/articles/leapfrog4.pdf

As far as your original question goes, all Adams methods use a "single" function evaluation per time step, but they also store a time history to attain higher-order.

Is this the same method as on "page 2" of these slides? Given Halim's description, this looks an awful lot like a two-stage Runge-Kutta method ...

http://www.win.tue.nl/casa/meetings/seminar/previous/_abstract020130_files/talksem.pdf

@David Seal.

As you have already pointed out, the Adams-Bashforth-Moulton multi-step schemes do store a time history of the evaluations of f. Because of that, it is not straightforward to change the time step during the computations. This contrasts with the single-step method proposed by Ulrich (which he calls ALF for Asynchronous Leap Frog).

Furthermore, if I am not too wrong, the 2-stage Runge-Kutta methods you are referring to are using two distinct evaluations of f per time step. Therefore, the ALF method does not seem to belong to this category. As described in my precedent post, ALF appears to use two interleaved leap-frog steps, one of which uses an evaluation of f while the other uses an extrapolated slope. This is the striking particularity of this method and, at the same time, it might be a source of serious difficulty as I have pointed out in my above post.

In any case, I believe that ALF rather belongs to the category of leap-frog methods and not to 2-stage Runge-Kutta methods. However, please note that I have never heard of ALF until this week. So I will be waiting to see what Ulrich has to say about all this.

@H.E. Lehtihet

A couple of comments.

1.) Very good. It looks like the method is indeed second-order, and I see that you only need one function evaluation per time step because you're evolving not only y, but the second quantity that you're calling S as well.

2.) If you walk through the steps for doing a stability analysis, it's interesting that the stability region is identical to the explicit midpoint method.

@David Seal.

I am glad to hear that you've already checked the stability region of ALF. I intended to do that next week-end. I am not surprised to hear that the stability region is the same as that of E2LF. However, as I have pointed out earlier, E2LF is known for its weak instability and thus requires a smoothing process every N steps with N even (N about 20 or 30). Therefore, if ALF belongs to the same category, it is reasonable to think that it also has the same inherent problem. Fortunately, Ulrich has already proposed a couple of techniques to circumvent this potential difficulty (see his last post).

@David,

excellent observation. It's exactly the fact that together with y also S evolves (and thus can be considered part of the 'state') that allows to have second order at a single input of new information (evaluation of f).

@all

It is a remarkable property of the ALF-algorithm that appart from y and S it needs no further memory: Let me write it down in the form I love most and which includes the relaxation parameter r

To all readers of my previous post:

Once, the post was longer, and in an attempt to add one more link to a list of links it was cut down to the size of the edit window. Silly software or silly user?

So I continue.

... from Wolfram Research can enjoy playing with ALF in the programs for which links are given subsequently:

http://demonstrations.wolfram.com/LeakingBucketEquationAndReversibility/

http://demonstrations.wolfram.com/IntegratingBeyondInfinityAndBack/

http://demonstrations.wolfram.com/IntegratingAcrossSingularities/

http://demonstrations.wolfram.com/ComparingLeapfrogMethodsWithOtherNumericalMethodsForDifferen/

http://demonstrations.wolfram.com/OnTheStabilityLimitOfLeapfrogMethods/

http://demonstrations.wolfram.com/RelativisticQuantumDynamicsIn1DAndTheKleinParadox/

http://demonstrations.wolfram.com/TheAsynchronousLeapfrogMethodAsAStiffODESolver/

The boost.numeric.odeint C++ library is "for solving initial value problems (IVP) of ordinary differential equations". Here is an overview of the implemented algorithms: http://www.boost.org/doc/libs/1_53_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/getting_started/overview.html. The "Modified Midpoint" algorithm is Order 2, and you may like to look up these other algorithms with "configurable" Order, whether they can be second order.

@Thomas

I think this was the link that you intended to post:

http://www.boost.org/doc/libs/1_53_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/getting_started/overview.html

Hallo Thomas,

good to hear from you.

Unfortunately the link works not for me. I'll try it again tomorrow, perhaps the boosters have temporary problems. I would be totally surprised to find ALF in a family of configurable order algorithms since I was, despite some effort, unable to generalize ALF for higher order.

Crank nicolson is the method of choice

@Azzeddine. This was not the question. ALF is explicit, CN implicit, quite different worlds!

Halim gave me a predator-prey-problem (Lotka-Volterra system) which turned out to be really hard for ALF. With my present insight I could foresee the problem. Whenever the exact solution is close to constant, a reversible method (which ALF is) has problems. In the Lotka-Volterra system under consideration there is a phase where the two populations are nearly constant at a low level. There is a rather natural way to destroy reversibility without changing order and overall accuracy and this method blew away the problems seen before. Let me write the modified method down here:

We recall ALF (without the relaxation factor r which we don't need)

t+=h; y+=hS; S+=2 ( f(t, y) - S); y+=hS; t+=h; ( h=dt/2) .

The modification has two such steps in sucession and an averaging of two S's at the end:

t+=h; y+=hS; S+=2 ( f(t, y) - S); y+=hS; t+=h; ( h=dt/4)

Sm=S;

t+=h; y+=hS; S+=2 ( f(t, y) - S); y+=hS; t+=h;

S=(S+Sm)/2

This was the form to show the structure, what follows is the form for coding:

t+=h; y+=hS;

S+=2 ( f(t, y) - S); y+=2hS; t+=2h;

Sm=S; (* extra memory for the velocity *)

S+=2 ( f(t, y) - S); y+=hS; t+=h;

S=(S+Sm)/2

This is the 'averaged densified asynchronous leapfrog' method: ADALF

One ADALF-step has two function evaluations such as Heun and explicit midpoint (RK2) so that it makes a fair comparison to use ADALF, Heun, and RK2 with equal time step. In Halim's Lotka-Volterra system, ADALF is much more accurate than Heun and RK2.

Ulrich, That is an interesting paper. Thanks for sharing.

I think Von Neuman (spelling ?) stability analysis (which is for PDEs) would be helpful for your analysis/understanding. In particular, as you have noticed, leapfrog and your method are really for solving a second order ODE. You can apply them to a 1st order ODE, but they still have two 1st order ODE solutions lurking in them. For a 1st order ODE one of those solutions the right one and one of them is spurious and typically unstable. The spurious solution is excited by the initial conditions (as you noticed) and round-off error (at long times). One fix is to never use old information (at time n-1) for a 1st order ODE. RK methods for example do not look back. They do not have spurious solutions of this type. The well posed initial conditions for a 1st order ODE are really telling you to use 1 value of initial data (not 2).

Blair,

you bring up interesting ideas that I understand only to a small part. I suspect that I misguided you by not mentioning the initialization step S = f (t0,y0) in my previous post. With that step done, I see no room for a lurking second solution. I never observed suszeptibility to round-off errors. Actually ALF is reversible and so a solution that had 'gone half to infinity' goes strictly back after reversion of dt (which it would not if the trajectory would be spoilt by numerical noise).

Please understand that I'm here explaining my view and that this does not mean that I'm not interested in your view. The contrary is true! For an efficient communication (that I would very much like to see happen) it would probably help if you could scan through the paper for which a link is given yet in the question that starts the thread. Here it is again:

http://www.ulrichmutze.de/articles/leapfrog4.pdf

This mentions an (at least historically ) important application of the basic leapfrog method to the time-dependent Schroedinger equation (1st order).

Perhaps also the Wolfram demonstrations listed in an earlier post of mine are of interest. One of them is on stability.

I would be glad to get your second oppinion.

Ulrich, your paper "leapfrog4" is what I read. I liked it because I have also wondered what is the equivalent of Stormer-Verlet for a 1st order ODE. You can guess that it is probably leapfrog (or implicit midpoint), but I never saw a direct connection until your paper.

What I am telling you about the fact that you are solving two ODES (one correct and one spurious) is not my invention at all. It is pretty classic stuff in textbooks on numerical methods. If you google leapfrog and von Neumann you will probably get it (maybe even on Wikipedia). This is true of all methods that "look back" when applied to 1st order ODEs. Simple 2nd-order accurate Adams-Bashforth uses (n-1) information and has exactly the same issues.

When you solve a 2nd order ODE (on paper) you have two solutions and 2 initial conditions. For example: solutions of sin(at) and cos(at) for d2u/dt2 = -a2 u .

(those 2's are all superscripts) The initial conditions will determine how much of each solution is in the final answer.

A 1st order ODE (on paper) has one (homegeneous) solution that has an amplitude fixed by one initial condition.

The same happens with Ordinary Difference Equations (which is what we deal with).

A difference equation with two previous times in it (like midpoint, and Adams-Bashforth) has two fundamental solutions.

It also needs two things (essentially initial conditions) to start the method. Those two things set the amplitudes of the two solutions.

u^(n+1) - u(n-1) = -2*dt*lambda*u^(n) (the ^ means superscript) has a very simple right hand side that should produce exponential decay. It has the exact solution.

u^n = A*(sigma_1)^n + B*(sigma_2)^n

where sigma^2 - 1 + 2*dt*lambda*sigma = 0 has two solutions (sigma_1 and sigma_2).

The values of A and B (above) are set by the initial conditions. Usually A is large

and B is quite small.

You will find that sigma_1 is less than 1 and so for large n this part of the solution goes to zero (decays as it should) even when A is large.

You will find that sigma_2 is less than -1, so for large n this part of the solution ends up winning eventually (even when B is small). You will also find this second solution essentially corresponds to solving a second 1st-order ODE (coupled to the first and correct ODE you are trying to solve).

I know this is all very fast. But it is in lots of books. Lax (Stanford University) may have figured all this out (30 years ago) and may have written a book. I learned it it grad school but don't know the book (maybe Moin).

Modern Physics loves 2nd order systems (Particles, Quantum, E and M, etc), so they obsess on symplectic and Stormer-Verlet and its many extensions. But engineers are simpler beasts. They deal with 1st-order ODEs a lot. And you will notice, they never user simplectic methods (you have already figured out why, symplecity doesn't mean much). Have fun.

Blair,

your lesson helped me to read my 'bible' Stuart/Humphries: 'dynamical systems and numeric analysis' with different eyes and this made me aware that you are completely right.

Concerning your last paragraph on physicists vs. engineers and 2nd order vs. 1st order I have a different perception. Engineers are used to work with 2nd order Newton/Lagrange equations wheras physicists work with their 'holy' 1st order Hamilton equations. Only an observation!

@All

I thank you for your making me aware that the kind of stability that is found both in normal leapfrog and in ALF is not really satisfactory for a general purpose method. (Halim's predator-prey problem was particularly effective in convincing me.) As communicated in an earlier post, without adding to the computational burden of ALF one is able to increase stability considerably by adding an averaging step (one addition gets saved by doing two ALF steps always in combination). This defines ADALF ('averaged densified asynchronous leapfrog' ). What I have proved today is that in the standard stability analysis of ADALF the 'region of absolute stability' contains the unit half circle {z \in C| Re(z)

To all who would like to see stability analysis working interactively, here is a link to my most recent 'computable document' on the topic.

http://www.ulrichmutze.de/misc/spectralAnalysis5.cdf

ALF and ADALF working in predator prey population dynamics

http://www.ulrichmutze.de/misc/predatorPrey2.cdf

The newest status of my leapfrog integrators is now summerized in an electronic article on my homepage (see link below).

Unfortunately it has a bit invonvenient size of 14 MB, so be patient in loading it down.

The size comes from many Mathematica-created figures in PDF format.

Please, let me know any remarks, corrections, critique, ... and, again, many thanks to those who contributed already.

http://ulrichmutze.de/articles/leapfrog11.pdf

I'm happy to announce that a comprehensive description of the new integrator is now publicly available in the e-print arXiv.

http://www.arxiv.org/abs/1311.6602

@Ulrich, thanks for this thread! I found Your comment under other thread, based on Your answer!

Thanks Ulrich, downloaded arXiv version. You don't perhaps have a C source code implementation?

Ludwig,

what I presently have ready, are codes in Mathematica and in C++. Unfortunately I'm extremely bussy the next week. Then, however, I will send you a C version within a few days.

Ulrich, thanks that will be great, C++ is fine as well.

Ludwig,

good that C++ is fine for you. This makes things much easier and allows close to instant delivery. See the appended file, which has all required information in it.

Notice that the class template valarray, available in C++, makes the main difference to a mere C implementation. The present thread emerged from another one:

https://www.researchgate.net/post/Which_is_the_most_interesting_and_widely_used_data_structure

in which I claimed that arrays are by far the most important data structures for supporting algorithms in physics. Asked for an example I gave a definition of the ALF algorithm there.

Thanks Ulrich.

Ludwig, you are welcome.

The arxiv article 1311.6602 mentioned in my contribution of November 27, 2013 now got a new version (replacing the old version, same URL) which was triggered by valuable suggestions that I got from H. E. Lehtihet and from J.M. Sanz-Serna. These suggestions are gratefully acknowledged!

The fist scientific paper which bases its numerical part on the method mentioned in my previous contribution has now been published:

O. Jonasson, F. Karimi, I. Knezevic

Partially coherent electron transport in terahertz quantum cascade lasers based on a Markovian master equation for the density matrix

Journal of Computational Electronics/Springer US    (August 2016)

The preprint version is easy to find in arXiv.

I asked author Jonasson via RG-message about his experiences with the method, and this was his answer:

"...My experience was very good. With other methods (simple Euler and some second order methods), the simulations all diverged after about 1000 time steps, while your method did not for at least 100,000 time steps. ..."

@all

Recently I found these citations of my ALF-method:

arXiv:2202.02435v1 [cs.LG] 4 Feb 2022

arXiv:2102.04668v2 [cs.LG] 3 Mar 2021

On my webpage www.ulrichmutze.de there is now an animation section in which most of the animated graphics are created by ALF-programs

(C++20, making use of 'concepts').

In retrospect I feel obliged more than ever to acknowledge the furthering

contributions of H.E. Lehtihet and J.Blair Perot.

The most recent version of my ALF-article is

arXiv:1311.6602v4 [math.NA] 25 Apr 2016