The numeric results computed above seem to agree with the following results from symbolic dsolve.

exact0.mw

acer

Apart from reducing Digits (ie. keeping it at the default, 10) there may be another easy way to speed up some of this kind of calculation.

The idea is simple, if the parameter(s) hasn't changed very much -- from the previous iteration to the current one -- then the numerical solver might sometimes be able to converge to a new root quite quickly using the previous iteration's root as the new starting guess.

Usually, fsolve will try and pepper the multivariable solution space with random starting points, and we hope that one will converge. (This is also why it takes long to fail completely, by the way.)

So an ajustment to the methodology could be like this: first try the previous iteration's root as the new starting point. If that fails to converge then call fsolve without specifying a starting point, and let it proceed as usual.

At Digits=10 the old methodology takes 40 seconds to get all those 45 discovered roots for `n1` between -3.99 and 7.01. domain_try1.mw

Using the above methodology it takes about 3 seconds.  domain_try2.mw

Computations run in 64bit Maple 16.01 on Windows 7.

acer

Apart from reducing Digits (ie. keeping it at the default, 10) there may be another easy way to speed up some of this kind of calculation.

The idea is simple, if the parameter(s) hasn't changed very much -- from the previous iteration to the current one -- then the numerical solver might sometimes be able to converge to a new root quite quickly using the previous iteration's root as the new starting guess.

Usually, fsolve will try and pepper the multivariable solution space with random starting points, and we hope that one will converge. (This is also why it takes long to fail completely, by the way.)

So an ajustment to the methodology could be like this: first try the previous iteration's root as the new starting point. If that fails to converge then call fsolve without specifying a starting point, and let it proceed as usual.

At Digits=10 the old methodology takes 40 seconds to get all those 45 discovered roots for `n1` between -3.99 and 7.01. domain_try1.mw

Using the above methodology it takes about 3 seconds.  domain_try2.mw

Computations run in 64bit Maple 16.01 on Windows 7.

acer

@jschulzb Can we suppose that your newer question is your focus for this issue now, and that this earlier thread is not so pressing?

By the way, you can use the "branch" item and the botton of a post/question/answer is you would like to ask some new related question and have mapleprimes automatically put cross-referencing links on both of old and new.

The combination of dsolve events and parameters is very powerful. That, along with a good numeric DAE solver, makes the cornerstone of MapleSim, IMO. But the documentation and examples for events/parameters in Maple are a little thin relative to its power and usefulness.

PatrickT has asked some questions here about `events` and `parameters` in the dsolve numeric solver, which might help shed some light.

Without seeing you full code, I might hazard a guess that you are going for some implicitplot, and just need a little help figuring out the beast that could be plotted (since it will be a procedure that sets the parameters according to the arguments, and then runs until events hit/miss, and then returns some indicative scalar... or something like that).

@jschulzb Can we suppose that your newer question is your focus for this issue now, and that this earlier thread is not so pressing?

By the way, you can use the "branch" item and the botton of a post/question/answer is you would like to ask some new related question and have mapleprimes automatically put cross-referencing links on both of old and new.

The combination of dsolve events and parameters is very powerful. That, along with a good numeric DAE solver, makes the cornerstone of MapleSim, IMO. But the documentation and examples for events/parameters in Maple are a little thin relative to its power and usefulness.

PatrickT has asked some questions here about `events` and `parameters` in the dsolve numeric solver, which might help shed some light.

Without seeing you full code, I might hazard a guess that you are going for some implicitplot, and just need a little help figuring out the beast that could be plotted (since it will be a procedure that sets the parameters according to the arguments, and then runs until events hit/miss, and then returns some indicative scalar... or something like that).

Could you upload the whole working code that you have so far, including the DEs. It makes it so much easier to assist meaningfully and avoid wasting time on the wrong goals.

acer

@jschulzb Yes, it is easy to adjust it to discover points that attain some other value. Simply create a new procedure whose result is adjusted by the target value -- thus becoming a root-finding (zero-finding exercise).

If you are trying to discover where a solution from dsolve/numeric attains some specific value then the fastest way might be to properly use dsolve,events (which were implemented for just this kind of thing). Done right, that should be able to outperform rootfinding "after the fact" via fsolve or whatever else, both in terms of accuracy and speed.

@jschulzb Yes, it is easy to adjust it to discover points that attain some other value. Simply create a new procedure whose result is adjusted by the target value -- thus becoming a root-finding (zero-finding exercise).

If you are trying to discover where a solution from dsolve/numeric attains some specific value then the fastest way might be to properly use dsolve,events (which were implemented for just this kind of thing). Done right, that should be able to outperform rootfinding "after the fact" via fsolve or whatever else, both in terms of accuracy and speed.

@barefoot1980 Are you using Maple 13? If so, then you could try it instead as,

sol := solve({expand(eq1) < 1, expand(eq2) < 1}, {q});

Incidentally, the RootOfs of cubics that may appear for this particular example may be expressed explicitly and exactly,

evalc(allvalues([sol]))[];

@barefoot1980 Are you using Maple 13? If so, then you could try it instead as,

sol := solve({expand(eq1) < 1, expand(eq2) < 1}, {q});

Incidentally, the RootOfs of cubics that may appear for this particular example may be expressed explicitly and exactly,

evalc(allvalues([sol]))[];

@jschulzb Note that DEplot and DEplot3d can do animations too.

That may well be an effective and reasonably efficient way to show evolution in time from a/some initial conditions (ICs).

If, on the other hand, you intend on animating with respect to some other parameter that appears in an IVP's DEs or in its ICs and if you need it quicker than DEplot offers then you might take a glance at this.

@jschulzb Note that DEplot and DEplot3d can do animations too.

That may well be an effective and reasonably efficient way to show evolution in time from a/some initial conditions (ICs).

If, on the other hand, you intend on animating with respect to some other parameter that appears in an IVP's DEs or in its ICs and if you need it quicker than DEplot offers then you might take a glance at this.

@maple fan I'm afraid that a bug in the memory management of high precision LinearAlgebra causes the calculation to use up prohibitively too much memory in Maple 13. It is fixed in Maple 16, which is what I used above.

You can do the calculation quickly at hardware double precision at Digits=15, even in Maple 13, though, and see the relationship between accuracy and the reported condition numbers. Basically, there is an approximate near linear relationship between Digits and the error of a computed eigenvalue and the log[10] of its condition number.

For your example, the worst condition numbers are about 1e-13 and so when the calculation is done using fast hardware double precision (about Digits=15 in base 10) then there are only a couple of accurate decimal digits in those most troublesome of the computed eigenvalues. The best eigenvalues have condition number of about 1e-2 and the double precision results are likely accurate to about 12-13 decimal digits, etc (I haven't measured them precisely for your example, but it should be near that ballpark). So to get even the worst conditioned eigenvalues accurate to over 10 decimal digits the computation can be done with a working precision of Digits around 25-30 or so. I used Digits=35 above and seemed to get all roots accurate to 15 decimal places.

Unfortunately Digits=15 is as high as you can go before needing "software, high precision", where the mentioned bug exists before Maple 16. (Actually, I recall that the bug may have only appeared about Maple version 10 or 11.)

If you need it quick then in Maple 13, you can probably compute all the roots to about 10 correct places by using EigenConditionNumbers at Digits=15 and then using something like RootFinding:-Analytic on suitably restricted complex ranges to get the badly conditioned eigenvalues/roots (ie. those for which the condition number is worse that about 1e-4 or 1e-5).

In my experience all this is often not necessary, and fsolve can compute all roots of even badly conditioned univariate polynomials. I consider you example out of the ordinary, hitting some bug in the internal routine `fsolve/polyill` where it may be getting into a nonterminating loop or something.

If the condition numbers are all at least better than, say, 1e-14 then you may get about 1 correct digit in even the worst eigenvalue in the results from executing EigenConditionNumbers at Digits=15 and fast hardware double precision. If that is true then even if not fully accurate those eigenvalue estimates may be close to giving you a workable estimate of the spectral radius. Maybe.

@maple fan I'm afraid that a bug in the memory management of high precision LinearAlgebra causes the calculation to use up prohibitively too much memory in Maple 13. It is fixed in Maple 16, which is what I used above.

You can do the calculation quickly at hardware double precision at Digits=15, even in Maple 13, though, and see the relationship between accuracy and the reported condition numbers. Basically, there is an approximate near linear relationship between Digits and the error of a computed eigenvalue and the log[10] of its condition number.

For your example, the worst condition numbers are about 1e-13 and so when the calculation is done using fast hardware double precision (about Digits=15 in base 10) then there are only a couple of accurate decimal digits in those most troublesome of the computed eigenvalues. The best eigenvalues have condition number of about 1e-2 and the double precision results are likely accurate to about 12-13 decimal digits, etc (I haven't measured them precisely for your example, but it should be near that ballpark). So to get even the worst conditioned eigenvalues accurate to over 10 decimal digits the computation can be done with a working precision of Digits around 25-30 or so. I used Digits=35 above and seemed to get all roots accurate to 15 decimal places.

Unfortunately Digits=15 is as high as you can go before needing "software, high precision", where the mentioned bug exists before Maple 16. (Actually, I recall that the bug may have only appeared about Maple version 10 or 11.)

If you need it quick then in Maple 13, you can probably compute all the roots to about 10 correct places by using EigenConditionNumbers at Digits=15 and then using something like RootFinding:-Analytic on suitably restricted complex ranges to get the badly conditioned eigenvalues/roots (ie. those for which the condition number is worse that about 1e-4 or 1e-5).

In my experience all this is often not necessary, and fsolve can compute all roots of even badly conditioned univariate polynomials. I consider you example out of the ordinary, hitting some bug in the internal routine `fsolve/polyill` where it may be getting into a nonterminating loop or something.

If the condition numbers are all at least better than, say, 1e-14 then you may get about 1 correct digit in even the worst eigenvalue in the results from executing EigenConditionNumbers at Digits=15 and fast hardware double precision. If that is true then even if not fully accurate those eigenvalue estimates may be close to giving you a workable estimate of the spectral radius. Maybe.

@Cody With a few edits and corrections, your PSeries can do both, sure.

> PSeries:=proc(cfn,var,minreg::integer,maxreg::integer)
>   local r,imaxreg,iminreg;
>     if minreg>maxreg then
>       imaxreg:=minreg; iminreg:=maxreg;
>     else
>       imaxreg:=maxreg; iminreg:=minreg;
>     end if;
>      add(cfn(r)*var(r),r=iminreg..imaxreg);
>   end proc:
 
> PSeries(i->1/i!, R->x^R, 2, 0);

                                             2
                                1 + x + 1/2 x

 
> PSeries(i->(-1)^(i-1), R->1/R^2, 4, 1);

                                      115
                                      ---
                                      144

 
> h:=proc(n::posint)
>   local i;
>   add( (-1)^(i-1) * 1/i^2 , i=1..n );
> end proc:
 
> seq(h(i),i=1..6);

                                  31  115  3019  973
                          1, 3/4, --, ---, ----, ----
                                  36  144  3600  1200

> seq( PSeries(i->(-1)^(i-1), R->1/R^2, i, 1), i=1..6);

                                  31  115  3019  973
                          1, 3/4, --, ---, ----, ----
                                  36  144  3600  1200