@"By the way, on what grounds 'mp' inherits not 'double' ?"

This is very good question. I am asking it myself everyday.

From one point of view it would make existing code porting to mp much easier at some cases.

However, in the same time, this might cause hard-to find errors of accuracy loss.

For example, MATLAB scripts include many use cases of:

zero(m,n,class(a))

ones(..,class(a))

eps(class(a))

If class(a) return "double" for mp-objects - all the commands will create double objects silently, which is obviously will lead to hard-to-find errors. Especially with "eps".

Perfect way would be to overload the built-ins "zeros, ones, eps, etc." completely, for all arguments (not only mp). And add support for 'mp' class to all these functions. However MATLAB screams on this overloads.

Do you know, or could you suggest on how to overload all the system functions, so that MATLAB will not go crazy with all the warnings?

***

Cast is in my TODO list - thank you for sharing your code!

I will make this thread public, if no objections.

Actually it might indicate few things (in case of ordqz) - problem is too ill-conditioned and transformed matrices might be too far from generalized Schur or just indicate wrong inputs (mismatching matrix dims or else). In any case toolbox will just return 1 in case of extra output parameters.

Thank you for reminding, I almost forget to add this to toolbox - will add right now.

Yes, absolutely! Just forgot to add it to the list.

Actually I saw the syntax in GCARE too :).

Not sure if I get the meaning of the output - how re-ordering could fail? Probably will just return 1.

The issue was probably related to the fact that we didn't digitally sign all components of toolbox.

Now we do sign all the components and probably this is resolved.

Let me know if it still occurs.

This bug has been resolved in Windows version. You can download it from our website.

GNU Linux & Mac OSX will follow today.

In power function we handle inputs like complex numbers by default. Semantic for various singular cases in complex arithmetic is not defined uniformly and implementation dependent (e.g. see here: http://www.advanpix.com/2015/10/19/devnotes-1-multiplicative-operations-with-complex-infinity/)

That is the reason for such misbehavior. We will fix this today.

Thank you for reporting the issue

I am not sure, probably you need to use better precision (from numerical point of view).

But why don't you compute the integral analytically?

Try using Maple or Mathematica. I have just tried Maple - it gave long but sensible expression for the any r,p,q,s.

Hi again,

I couldn't find your last comment (did you delete it?).

Just few comments on the code.

The function should be able to accept the vector-parameters (quadqk is a vectorized method). So instead of

f = @(x)(cos(p*x)+cosh(p*x)+q*sin(p*x)+q*sinh(p*x))*(cos(r*x)+cosh(r*x)+s*sin(r*x)+s*sinh(r*x));

use element-wise multiplication:

f = @(x)(cos(p*x)+cosh(p*x)+q*sin(p*x)+q*sinh(p*x)).*(cos(r*x)+cosh(r*x)+s*sin(r*x)+s*sinh(r*x));

Also be careful with precision and interval counts. My example shows extreme values, please use more reasonable at first, e.g.:

mp.Digits(34);
...
[q,errbnd] = quadgk(f,mp(0),L,'RelTol',100*mp.eps(),'AbsTol',100*mp.eps(), 'MaxIntervalCount', 650)

Quadgk will warn you if it needs more levels of recursion (MaxIntervalCount)