Advanced Search

Axel Vogt

5906 Reputation

20 Badges

19 years, 305 days
Munich, Bavaria, Germany
Contact Axel Vogt

MaplePrimes Activity


These are answers submitted by Axel Vogt

formula...

Axel Vogt 5906
January 04 2024
2 10

@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

you can stay with Sum ......

Axel Vogt 5906
December 30 2023
1 0

... 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

Roughly ......

Axel Vogt 5906
December 27 2023
0 0

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...

Axel Vogt 5906
October 22 2023
1 0
> 

# 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

Some remarks...

Axel Vogt 5906
September 14 2023
0 0

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" ......

Axel Vogt 5906
August 26 2023
1 0

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...

Axel Vogt 5906
August 10 2023
2 0

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...

Axel Vogt 5906
August 05 2023
1 4

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

brackets...

Axel Vogt 5906
August 02 2023
1 0

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

Task2...

Axel Vogt 5906
July 25 2023
0 0

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

Task 1...

Axel Vogt 5906
July 24 2023
1 0

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

one way...

Axel Vogt 5906
June 24 2023
1 0
> 

# 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

and here is my answer...

Axel Vogt 5906
June 07 2023
2 1


 

> 

# https://www.mapleprimes.com/questions/236532-Workaround-To-Integrating-This
restart;

> 

integrand:=(ln(-x + sqrt(x^2 - 1))*x + sqrt(x^2 - 1))*x*sqrt(-x^2 + 1)/sqrt(x^2 - 1):
Int(integrand,x);
IntegrationTools:-Expand(%);
value(%):
IntegrationTools:-Parts(%, ln(-x + sqrt(x^2 - 1))):
res:=simplify(%);: #evala(%);

Int((ln(-x+(x^2-1)^(1/2))*x+(x^2-1)^(1/2))*x*(-x^2+1)^(1/2)/(x^2-1)^(1/2), x)

  

Int(x^2*(-x^2+1)^(1/2)*ln(-x+(x^2-1)^(1/2))/(x^2-1)^(1/2), x)+Int(x*(-x^2+1)^(1/2), x)

  

(1/9)*(-x^2+1)^(1/2)*(3*x^3*ln(-x+(x^2-1)^(1/2))+4*(x^2-1)^(1/2)*x^2-(x^2-1)^(1/2))/(x^2-1)^(1/2)

 (1)
> 

'diff(res, x) - integrand': '%'= simplify(%);

diff(res, x)-integrand = 0

 (2)
> 

 


 

Download MP_236532_some_integral.mw

well ......

Axel Vogt 5906
June 03 2023
0 0

Maple may decide for a "solving approach" depending on (pattern of) the input. And evala(expr) resp. expr are different.

May be that for *this* evala(expr) the algorithm quickly decides that it will not try to succeed.

If you really want to know (and suffer) in detail you can use "infolevel".

c[2] = ln(a)...

Axel Vogt 5906
May 30 2023
1 2

So this is a bug report for "int" ...

writing c[2] = ln(a) is a workaround for that bug

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