Axel Vogt

5816 Reputation

20 Badges

19 years, 104 days
Munich, Bavaria, Germany

MaplePrimes Activity


These are answers submitted by Axel Vogt

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

 

 

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


 

# https://www.mapleprimes.com/questions/237741-Integral-Of-Dirac
restart; kernelopts(version);

`Maple 2023.2, X86 64 WINDOWS, Nov 24 2023, Build ID 1762575`

(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(%);

proc (x) options operator, arrow; x^2+y^2-1 end proc

(2)

Int(Dirac('g'(x)), x=-infinity ... infinity): '%'= value(%);

Int(Dirac(g(x)), x = -infinity .. infinity) = 1/abs(y^2-1)^(1/2)

(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);

0 = g(x)

 

{(-y^2+1)^(1/2), -(-y^2+1)^(1/2)}

 

(-y^2+1)^(1/2), -(-y^2+1)^(1/2)

(4)

 

The zeros of g are simple iff y <> -1 . 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);

delta(g(x)) = Sum(delta(x-xi)/abs((D(g))(xi)), `in`(xi, Zeros))

 

delta(x^2+y^2-1) = (1/2)*delta(x-(-y^2+1)^(1/2))/abs(y^2-1)^(1/2)+(1/2)*delta(x+(-y^2+1)^(1/2))/abs(y^2-1)^(1/2)

 

(1/2)*Dirac(x-(-y^2+1)^(1/2))/abs(y^2-1)^(1/2)+(1/2)*Dirac(x+(-y^2+1)^(1/2))/abs(y^2-1)^(1/2)

(5)

 

Then the integral works out as asserted (since Int(Dirac(x+constant), x = -infinity .. infinity) = 1 ):

 

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

Int(Dirac(g(x)), x = -infinity .. infinity)

 

`` = Int((1/2)*Dirac(x-(-y^2+1)^(1/2))/abs(y^2-1)^(1/2)+(1/2)*Dirac(x+(-y^2+1)^(1/2))/abs(y^2-1)^(1/2), x = -infinity .. infinity)

 

`` = 1/abs(y^2-1)^(1/2)

(6)
 

 


 

Download MP_237741-Integral-Of-Dirac.mw

@Zeineb I can not read your local pdf ...

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

... if you change the range of integration


 

# https://www.mapleprimes.com/questions/237601-Message--Error-in-Assuming-When

restart; kernelopts(version); N:='N':

`Maple 2023.2, X86 64 WINDOWS, Nov 24 2023, Build ID 1762575`

(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);

Int((1+x)^(1/2)*phi(n, x)*(-x^2+1)^(1/2), x = -1 .. 0)+Int(phi(n, x)*(-x^2+1)^(1/2), x = 0 .. 1)

 

`` = Int(x*phi(n, x-1)*(-x+2)^(1/2)+phi(n, x)*(-x^2+1)^(1/2), x = 0 .. 1)

 

proc (n) options operator, arrow; Int(x*f(n, x-1)*(-x+2)^(1/2)+f(n, x)*(-x^2+1)^(1/2), x = 0 .. 1) end proc

(2)

S:=(N,t) -> Sum(c(n)*f(n, t), n = 0 .. N);

proc (N, t) options operator, arrow; Sum(c(n)*f(n, t), n = 0 .. N) end proc

(3)

'S(1,1)': '%'=  simplify(value(%)); `` = evalf(rhs(%));

S(1, 1) = (1/210)*(105*Pi-536+704*2^(1/2))/Pi

 

`` = 1.19665354425139

(4)
 

 


 

Download MP_237601.mw

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.

# https://www.mapleprimes.com/questions/237170-How-Can-I-Resolve-This-Numeric-Integral

restart; Digits:=15;

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);

proc (c) options operator, arrow; 1/erfi(c) end proc

 

proc (x, c) options operator, arrow; exp(-x^2-(1/2)*erfi(x)^2*E(c)^2) end proc

 

``

 

proc (c) options operator, arrow; proc (x) options operator, arrow; g(x, c) end proc end proc

 

proc (c) options operator, arrow; proc (x) options operator, arrow; evalf[15](g(x, c)) end proc end proc

 

``

 

proc (alpha, b) options operator, arrow; 2*sqrt(2)*alpha*E(sqrt((1/2)*alpha)*b)*(Int((h(sqrt((1/2)*alpha)*b))(x), x = 0 .. infinity))/Pi end proc

 

proc (alpha, b) options operator, arrow; 2*sqrt(2)*alpha*E(sqrt((1/2)*alpha)*b)*(Int(hNum(sqrt((1/2)*alpha)*b), 0 .. infinity, method = _d01amc))/Pi end proc

(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

 

0.999002158629825e-1

(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<b: op(%);

F(a, b) = a*F(1, sqrt(a)*b), F(a, b) = F(a*b^2, 1)/b^2

 

true, true

(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(%);

 


 

Download MP_237170_avoid_overflow.mw

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.

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.

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):

`Maple 2023.1, X86 64 WINDOWS, Jul 07 2023, Build ID 1723669`

(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);

proc (t) options operator, arrow; sqrt(t) end proc

(2)

Matrix data

   

 

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

[M1, M2, M3, M4]

 

LinearAlgebra:-Dimensions = [155, 85, 85, 36]

 

LinearAlgebra:-Eigenvalues = [{0, 3}, {0, 28}, {0, 20, 100, 252}, {0, 24, 29-181^(1/2), 29+181^(1/2)}]

(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);

"MinimalPolynomial"

 

[x*(x-3), x*(x-28), x*(x-100)*(x-20)*(x-252), (x-29+181^(1/2))*(x-29-181^(1/2))*x*(x-24)]

(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;

3

(5)

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

theMatrices[3]

 

{0, 20, 100, 252}

 

x*(x-100)*(x-20)*(x-252)

(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);

proc (t) options operator, arrow; t^3*a[3]+t^2*a[2]+t*a[1]+a[0] end proc

 

[seq(p(lambda) = f(lambda), `in`(lambda, EV))]

 

proc (t) options operator, arrow; t^3*(-1/121600+(1/185600)*5^(1/2)+(1/1481088)*7^(1/2))+t^2*(17/7600-(11/5800)*5^(1/2)-(5/61712)*7^(1/2))+t*(-63/1520+(63/464)*5^(1/2)+(125/92568)*7^(1/2)) end proc

(7)

 

Then  f(M) = p(M)

 

sqrtM:='p(M)'; # p to be understood in the algebra of square matrices

p(M)

(8)

'sqrtM.sqrtM - M';  evala(%): #fnormal(%): simplify(%, zero):
convert(%, list): #nops(%);
`set of all entries`=convert(%, set);

sqrtM.sqrtM-M

 

`set of all entries` = {0}

(9)

 


 

Download MP_236323_MatrixSquareRoot_-_Polynomial.mw

 

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):

10

(1)

f:= x -> x/(exp(x)-1);

proc (x) options operator, arrow; x/(exp(x)-1) end proc

(2)

M:= <<1, 0, 0> | <1, 1, I> | <3, 0, 2>>;
MM:=convert(M, Matrix, datatype=complex(sfloat));

Matrix(3, 3, {(1, 1) = 1, (1, 2) = 1, (1, 3) = 3, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = I, (3, 3) = 2})

 

Matrix(%id = 36893488148092049636)

(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);

proc (A) options operator, arrow; MF(A, 1, x) end proc

 

proc (A) options operator, arrow; A.(1/(MF(A, exp(x), x)-id(A))) end proc

(4)

#F(MM);
'Norm(MatrixFunction(M, f(x), x) - F(MM))': '%'= evalf(%);

Norm(MF(M, f(x), x)-F(MM)) = 0.683013562280646e-9+0.*I

(5)

 


 

Download MP_236851_MatrixFunction_-_short.mw

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  ?`;

S := 3*sin((1/3)*arcsin(t^(1/2)))/t^(1/2)

``

H := 1+(4/27)*(t*hypergeom([1, 4/3, 5/3], [2, 5/2], 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  ?`;

A := hypergeom([1, 4/3, 5/3], [2, 5/2], t)*t^2

B := (81/4)*(t^(1/2)*sin((1/3)*arcsin(t^(1/2))))-27*t/4

``

`A=B  ?`

(2)

A: convert(%, Sum, dummy=k) assuming abs(t)<1: convert(%, GAMMA): combine(%):
AS:=simplify(%);

AS := Sum((27/2)*(t^(2+k)*3^(1/2)*GAMMA(4/3+k)*GAMMA(5/3+k)*4^k/(Pi*GAMMA(4+2*k))), k = 0 .. infinity)

(3)

B: convert(%, Sum, dummy=k) assuming abs(t)<1: convert(%, GAMMA):
BS:=simplify(%);

BS := Sum(3*27^(-k)*4^k*t^(2+k)*GAMMA(3+3*k)/((4*k+6)*GAMMA(2+k)*GAMMA(2+2*k)), k = 0 .. infinity)

(4)

'AS-BS'; ``= combine(%);

AS-BS

`` = 0

(5)


 

Download MP_236788_unable_to_sum_Task2.mws

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/(((n*(n+1))^(1/2)*(n^(1/2)+(n+1)^(1/2))))

`` = 1/(((n+1)^(1/2)*n^(1/2)*(n^(1/2)+(n+1)^(1/2))))

`` = (-n*(n+1)^(1/2)+n^(3/2)+n^(1/2))/(n*(n+1))

`` = -1/(((n+1)^(1/2)))+n^(1/2)/(n+1)+1/((n^(1/2)*(n+1)))

`` = ((n+1)^(1/2)-n^(1/2))/(n^(1/2)*(n+1)^(1/2))

Sum(((n+1)^(1/2)-n^(1/2))/(n^(1/2)*(n+1)^(1/2)), n = 1 .. infinity) = 1

(1)

 


 

Download MP_236788_unable_to_sum_Task1.mws

# https://www.mapleprimes.com/questions/236653-Can-The-Output-Of-Int-Options-Be-Suppressed
Int(1/x, x = a .. 2):
%=int(op(%), 'AllSolutions');

Int(1/x, x = a .. 2) = piecewise(a < 0, undefined, a = 0, infinity, 0 < a, -ln(a)+ln(2))

(1)

 

Download MP_236653.mw

1 2 3 4 5 6 7 Last Page 1 of 92