## 5816 Reputation

19 years, 104 days
Munich, Bavaria, Germany

## formally...

For example like this (upload does not work ...):

NumericEventHandler(invalid_operation = `Heaviside/EventHandler`(value_at_zero = 0)):
assume(u::real, v::real);

[Int(Dirac(u-v), [u, v]) , min(u,v)]; value(%);
convert(%, piecewise, u);
%[1] - %[2]; # = 0

[Int(Dirac(u+v-1), [u, v]), max(u+v-1, 0)]; value(%);
convert(%, piecewise, u);
%[1] - %[2]; # = 0

## hm...

I think that limit(expr, +oo) = 24 is false

## Proof...

 > # https://www.mapleprimes.com/questions/237741-Integral-Of-Dirac restart; kernelopts(version);
 (1)
 > g:= x -> x^2+y^2-1; #g:= x -> x^2-a^2; 'y^2-1 = eval(-a^2, a=sqrt(1-y^2))'; is(%);
 (2)
 > Int(Dirac('g'(x)), x=-infinity ... infinity): '%'= value(%);
 (3)

This is correct using a quite common definition for composing Dirac and appropriate functions, see Maple's help

or https://en.wikipedia.org/wiki/Dirac_delta_function#Composition_with_a_function referencing to I M Gelfand

 > 0='g'(x); Zeros:={solve(%,x)}; x1,x2:=op(Zeros);
 (4)

The zeros of g are simple iff  . Therefore the composition can be defined as

 > delta('g'(x)) = Sum('delta(x-xi)/abs(D(g)(xi))', xi in 'Zeros'); subs(Sum=add, %): %; Dirac_of_g(x):=eval(rhs(%), delta=Dirac);
 (5)

Then the integral works out as asserted (since  ):

 > Int(Dirac('g'(x)), x=-infinity ... infinity); ``=Int(Dirac_of_g(x), x=-infinity .. infinity); value(%);
 (6)

## formula...

In Gradshteyn+Ryzhik 3.196.3 you find the following formula

Here B = Beta function and EH is a reference to Erdelyi's Bateman manuscripts.

From that you may derive your desired formula by changing the variable through x = f(y)
and renaming variables

## you can stay with Sum ......

... if you change the range of integration

 > # https://www.mapleprimes.com/questions/237601-Message--Error-in-Assuming-When
 > restart; kernelopts(version); N:='N':
 (1)
 > f := (n, x) -> 1/4*JacobiP(n, 1/2, 1/2, x)*4^(1/2)/   (GAMMA(n + 3/2)^2/((2*n + 2)*GAMMA(n + 2)*n!))^(1/2):
 > # change range of integration Int(sqrt(1 + x)*phi(n, x)*sqrt(-x^2 + 1), x = -1 .. 0) +   Int(phi(n, x)*sqrt(-x^2 + 1), x = 0 .. 1); ``=IntegrationTools:-Change( op(1,%), 1+x=xx,xx) +   op(2,%): subs(xx=x,%): combine(%); subs(phi=f, rhs(%)): c:=unapply('%',n);
 (2)
 > S:=(N,t) -> Sum(c(n)*f(n, t), n = 0 .. N);
 (3)
 > 'S(1,1)': '%'=  simplify(value(%)); `` = evalf(rhs(%));
 (4)

## Roughly ......

Some general explanation is given in  https://simple.wikipedia.org/wiki/Compiled_language

Roughly evalhf uses compiled code - but only for a small set of numerical function, results are limited w.r.t. "valid digits of the result"

For the full strength of Maple however one needs more, evalhf ist not enough.

## Acer's solution as numerical function an...

 > # https://www.mapleprimes.com/questions/237170-How-Can-I-Resolve-This-Numeric-Integral
 > restart; Digits:=15;
 (1)
 > E:= c -> 1/erfi(c); #c = sqrt(alpha/2)*b; g:= (x, c) -> exp(-x^2 - erfi(x)^2*E(c)^2/2); ``; h:=    (c) -> (x -> (g(x, c))); hNum:= (c) -> (x -> evalf[15](g(x, c))); # Acer ... ``; F:= (alpha,b) -> 2* sqrt(2)*alpha* E(sqrt(alpha/2)*b)/Pi*   Int( h(sqrt(alpha/2)*b)(x), x=0..infinity); P:= (alpha,b) -> 2* sqrt(2)*alpha* E(sqrt(alpha/2)*b)/Pi*   Int( hNum(sqrt(alpha/2)*b), 0..infinity,method = _d01amc);
 (2)
 > # One can show that F equals the original task using the change # sqrt(2)*sqrt(-alpha)*t/2 = I*x # erfi is a real valued function (convert to an integral)
 >
 > # after a first call for a 'warm up' the function P is ~ 20% slower than acer's direct version forget(evalf): CodeTools:-Usage(   evalf(P(0.1, 0.1)) );
 memory used=9.70MiB, alloc change=0 bytes, cpu time=124.00ms, real time=122.00ms, gc time=0ns
 (3)
 >
 > 'F(a,b)= a*F(1, sqrt(a)*b)', 'F(a,b)= 1/b^2*F(a*b^2,1)'; map(is, [%]) assuming 0
 (4)
 > # That identity means that F originates from an univariate function: # Let f(x) be given. Then define Phi(x,y) := f(x*y^2)/y^2 and f(x) = Phi(x,1)
 > plot(x -> P(x,1), 0 .. 20, numpoints=6, smartview=false);
 > # with some hand waving one sees: F(0,1) =, F(+oo, 1) = 0
 >
 > #a0, b0:= 0.1, 0.2; # example #['P(a0, b0) = a0*P(1, sqrt(a0)*b0)', 'P(a0, b0) = P(b0^2*a0, 1)/b0^2']; evalf(%);
 >

## Some remarks...

It depends on the user how Excel displays numbers in the GUI, especially for double precision floats. Note that this is only for  the display, not the number itself, internally the number does not change.

Also note that Excel (always) uses 15 decimals - using a "longer" display will be just be filled by trailing zeros. Also note that Excel internally may use more 'decimals' but does not allow direct access.

Note that there may be different presentations in Base 10 = decimal notation for the same Base-2 number = IEEE 754

Moreover Excel is not IEEE compliant.

This holds for the GUI, but also for Visual Basic code.

## "no way" ......

There is no way for a "formal result": essentially you have 2 equations which are like sums of polynom(A,B)*log(involving A,B).

You may try fsolve(..., complex). This may give 1 pair (A,B). However your system of intersecting (non-algebaic) curves may have many such points.

## a polynomial approach...

Here is the polynomial approach given in Higham's book "Functions of Matrices - Theory and Computation".

For the given matrices it is quite simple to use:

 > # https://www.mapleprimes.com/questions/236323-Why-Is-The-Performance-Of-LinearAlgebraMatrixPower
 > restart; #Digits:=15; kernelopts(version); # Maple 2023.1 with(LinearAlgebra):
 (1)

Find the "square root" of the matrices given below. This is a done using a polynomial approach given in Higham, Functions of Matrices - Theory and Computation (2008), Section 1.2.2. Polynomial Interpolation. It turns out that for these matrices M it is quite "simple" to compute f(M).

 > f:=t -> sqrt(t);
 (2)
 Matrix data

 > theMatrices:='[M1, M2, M3, M4]'; `Dimensions` = map('q -> Dimensions(q)[1]', theMatrices); `Eigenvalues` = map('q -> convert(Eigenvalues(q), set)', theMatrices);
 (3)

The matrices are medium sized but have very few different eigenvalues. Compute the minimal polynomials for the matrices

 > "MinimalPolynomial"; map('q -> MinimalPolynomial(q,x)', theMatrices): factor(%): map(evala@AFactor, %): # does not factor over Rationals ...   minPols:=convert(%, radical); #map(degree,%,x);
 (4)

The minimal polynomials are very simple: linear factors with multiplicity = 1 (BTW: that means that the blocks of the Jordan Normal Form are of size at most 1 or in less educated words: the matrices are diagonalizable).

Then it is easy to find a polynomial p such that p(M) = f(M) following the approach given in Higham's book.

As an example take the 3rd matrix:

 > Example:=3;
 (5)
 > M:='theMatrices'[Example]; EV:=Eigenvalues(M): EV:=convert(%, set); m:=minPols[Example]; deg:=degree(m,x):
 (6)

One has to find a (minimal) polynomial p such that p(j) (lambdai) = f(j)(lambdai) for i = all eigenvalues and derivates j = up to degree(factors of m) - 1, see Higham.

The degree of all the factors is 1. So it means: find the (minimal) approximating polynomial for the points f(lambda), lambda the eigenvalues. Either use a Maple package for that or do it directly, it is to solve a linear for the coefficients:

 > sum(a[i]*t^i, i=0 .. deg-1): p:=unapply(%, t); '[seq( p(lambda)=f(lambda), lambda in EV)]'; solve(%): # this is a linear system ... evala(%): eval(p(t), %): p:=unapply(%, t);
 (7)

Then  f(M) = p(M)

 > sqrtM:='p(M)'; # p to be understood in the algebra of square matrices
 (8)
 > 'sqrtM.sqrtM - M';  evala(%): #fnormal(%): simplify(%, zero): convert(%, list): #nops(%); `set of all entries`=convert(%, set);
 (9)
 >

## another way...

Another way to work around that bug

 > # https://www.mapleprimes.com/questions/236851-Is-LinearAlgebraMatrixFunction-Numerically
 > restart; Digits:=10; #kernelopts(version); # Maple 2023.1 with(LinearAlgebra): alias(MF=MatrixFunction):
 (1)
 > f:= x -> x/(exp(x)-1);
 (2)
 > M:= <<1, 0, 0> | <1, 1, I> | <3, 0, 2>>; MM:=convert(M, Matrix, datatype=complex(sfloat));
 (3)
 > id:= A -> MF(A, 1, x); # ~ identity Matrix F:= A -> A . (1/(MF(A, exp(x),x) - 1*id(A))); # ~ re-write f #F:= A -> A . (MF(A, exp(x),x) - MF(A, 1, x))^(-1);
 (4)
 > #F(MM); 'Norm(MatrixFunction(M, f(x), x) - F(MM))': '%'= evalf(%);
 (5)
 >

## brackets...

Use (,) not {,},[,] for computational brackets

Maple seems to use different approches. Rewriting both terms as sum one sees that they are equal

 > # https://www.mapleprimes.com/questions/236788-Limitations-Of-sum
 > restart;
 > Sum(binomial(3*n, n)*x^n/(2*n + 1), n = 0 .. infinity): value(%) assuming (abs(x) <= 4/27): subs(x=4/27*t, %): S:=simplify(%); ``; eval(binomial(3*n, n)*x^n/(2*n + 1), n = 0): # = 1 %+Sum(binomial(3*n, n)*x^n/(2*n + 1), n = 1 .. infinity): value(%) assuming (abs(x) <= 4/27): H:=subs(x=4/27*t, %); ``; `H=S  ?`;
 (1)
 > H: (%-1)/4*27/t: A:=%*t^(3/2)*sqrt(t); S: (%-1)/4*27/t: B:=%*t^(3/2)*sqrt(t): B:=expand(B); ``; `A=B  ?`;
 (2)
 > A: convert(%, Sum, dummy=k) assuming abs(t)<1: convert(%, GAMMA): combine(%): AS:=simplify(%);
 (3)
 > B: convert(%, Sum, dummy=k) assuming abs(t)<1: convert(%, GAMMA): BS:=simplify(%);
 (4)
 > 'AS-BS'; ``= combine(%);
 (5)

Here is how one could do for the 1st task

 > # https://www.mapleprimes.com/questions/236788-Limitations-Of-sum
 > restart;
 > 1/(sqrt(n*(n + 1))*(sqrt(n) + sqrt(n + 1))); ``=simplify(%) assuming 1<=n; evala(%);   expand(%); simplify(%); Sum(rhs(%), n=1 .. infinity):   %=value(%);
 (1)
 >