@Preben Alsholm

If you take an arbitrary polynomial, you're guaranteed to have only those relations between its roots that are unchanged under any permutation of the roots, such as a[1]*...*a[6]=free term. This polynomial has an extra relation ee=4, which is preserved only by some permutations. Exactly how this is tripping up evala, I have no idea.

To elaborate some more on acer's first comment, the relevant part of the code from LinearAlgebra:-Rank is:

svd := LinearAlgebra:-SingularValues(localM,('output') = ['S'],('conjugate') = true);
for i to min(m,n) do
  if svd[i] = 0. or 10^(Digits-1)*abs(svd[i]) < abs(svd[1]) then
    return i-1
  end if
end do;
return min(m,n)

So there doesn't seem to be any way to provide a non-default tolerance (except by lowering the value of Digits) or to use the function Testzero (which is user-settable).

Also, eval won't be able to recover the original expression if it contained lists or sets. That is, `[]`(1,2) is a true list, but op(0,[1,2])(1,2) isn't (indexing doesn't work, etc.).

@deniscr

You say that the procedure should work for other cases as well, but you still don't specify what types of matrices it should work with.

If your matrices are always real symmetric, then there is a criterion (Sylvester's law of inertia) to tell when they are congruent.

If they are real symmetric and congruent, then, writing C=B_C.D_C.B_C^t and M=B_M.D_M.B_M^t, we get

U=B_C.DiagonalMatrix(sqrt(Diagonal(D_M)/Diagonal(D_C))).B_M^t,

where B_C and B_M are matrices whose columns are normalized eigenvectors of C and M, and D_C and D_M are diagonal matrices of eigenvalues of C and M.

One difficulty is that if you want a real-valued U, then D_M/D_C should divide eigenvalues of the same sign, which is problematic for symbolic matrices.

Is anything more known about the matrices C and M? The solution doesn't always exist for arbitrary C and M. For instance, if C is symmetric, then U^t.C.U is symmetric, and if M isn't symmetric, the equation cannot hold. Do you know from the start that your matrices are congruent?

@vv

Yes, I think it could be useful to have series return more than one asymptotic "tail". So, for instance, one could obtain the asymptotics of BesselI(1/4, z)-BesselI(-1/4, z) at +infinity directly, without having to do a conversion into something else with a known expansion.

I don't know of any convenient way to turn off Executable Math. I'm using Ctrl-R/Ctrl-T to switch between the math mode and the text mode in a paragraph; it's not always clear which mode is on (so I have to look at the Context Bar), but other than that, it works reasonably well. But pressing Shift-F5 every time is just too much.

@carriewong

I think you should analyze the code to understand what it does. As a starting point, open the help page ?plot/parametric and read the first and the third bullet points in the description section.

@Jjjones98

But how are you going to write the code without knowing how the NC formulas are derived? Have you looked them up? (Starting with this, perhaps.)

Here are some hints.

The approximations that you want are called open Newton-Cotes formulas.

They're constructed in the same way as the closed NC formulas (which are available in Maple directly in Student:-Calculus1:-ApproximateInt), except the endpoints a and b are not used.

So, construct a list of data points [[a+h, f(a+h)], ..., [b-h, f(b-h)]].

Create the interpolating polynomial with CurveFitting:-PolynomialInterpolation.

Integrate the polynomial over the segment [a,b], simplify (alternatively, think about how to make everything simpler from the start by using integers instead of a+h, a+2h, ...), and you're done.

To estimate the error term, apply your code to f on the interval [0, h], subtract int(f(x),x=0..h) and expand the result into a series.

@vv

But then what about list[[3,2,1]] to do a permutation of elements? It's not redundant, as there is no  corresponding range syntax.

@dharr

list[range] is unambiguously documented in ?[]:

"The selection op(i..j,S) selects the sub-sequence (S[i],S[i+1],...,S[j]).

...

If S is a list then the selection S[i..j] selects the sublist [op(i..j,S)]."

Clearly this is useful.

@dharr

I'm not saying m[[1,2],2] is doing something random, the question is rather whether it should be consistent with the other three cases (and whether it's documented).

@dharr

It doesn't explain m[1..2, 2] though.

The documentation for rtables makes a distinction between integer and non-integer indices. i is likely to mean an integer index.

@ablmf

If evala@Norm works, then you can use algcurves:-puiseux:

poly := (sqrfree@evala@Norm)(x-nu3)[2, 1, 1];

algcurves:-puiseux(poly, u = 0, x, 0);
      {u*RootOf((w^2+2*w-1)*_Z^2+4*w^2*RootOf(_Z^2-2*_Z*w-2*w)+4*w*RootOf(_Z^2-2*_Z*w-2*w)+2*w^2-
      2*RootOf(_Z^2-2*_Z*w-2*w)-2*w)+RootOf(_Z^2-2*_Z*w-2*w), RootOf(_Z^4+(-4*w+8)*_Z^3+
      (4*w^2-28*w+16)*_Z^2+(40*w^2-16*w)*_Z-16*w^3+4*w^2)}

You can see that there are indeed two expansions: one that contains a linear term and one that doesn't (the next term in the second expansion will be quadratic in u).

I think this is fairly self-explanatory: nu3 can be viewed as a solution of a polynomial equation, and then the solution can be expanded into a Puiseux series.

@acer

I took this formula:

sum(Stirling1(n, i)*x^i, i = 1 .. n) = (-1)^n*GAMMA(n-x)/GAMMA(-x); # n>0

We need to set x=-1/a to match the sum in det2, and the rest is just shuffling around the terms.