Robert Israel

6577 Reputation

21 Badges

18 years, 218 days
University of British Columbia
Associate Professor Emeritus
North York, Ontario, Canada

MaplePrimes Activity


These are replies submitted by Robert Israel

It's not a bug.  When you ask for the integral with a float as upper or lower limit, Maple takes this as an indication you want a numerical answer, and calls the numerical integration methods, just as if you had entered

> evalf(Int(...));

On the other hand, if the upper and lower limits are rational, Maple tries to give an exact symbolic answer for the integral, in this case using polylogs.  Although the end result should be real, its computation involves complex numbers in the intermediate steps,  and roundoff error produces the small imaginary part.

 

Why not introduce c^2 + s^2 = 1 as an identity, and use solve (or Groebner bases)?

M := subs(cos(Theta)=c, sin(Theta)=s,M);
V:= <v1,v2,v3>;
eqs:= convert(M.V - lambda*V, set) union {c^2+s^2-1};
S:= remove(t -> (subs(t,{v1,v2,v3})={0}, [solve(eqs)]);



 

Why not introduce c^2 + s^2 = 1 as an identity, and use solve (or Groebner bases)?

M := subs(cos(Theta)=c, sin(Theta)=s,M);
V:= <v1,v2,v3>;
eqs:= convert(M.V - lambda*V, set) union {c^2+s^2-1};
S:= remove(t -> (subs(t,{v1,v2,v3})={0}, [solve(eqs)]);



 

> with(Statistics):
   CDF(FRatio(m,n),x) assuming x > 0, m::posint, n::posint;
 

x^((1/2)*m-1)*(m/n)^((1/2)*m)*n*(GAMMA((1/2)*m)*Pi*csc((1/2)*Pi*n)*x^(-(1/2)*m+1)*m^(-(1/2)*m+1)*n^((1/2)*m-1)-2*n^(-2+(1/2)*m+(1/2)*n)*x^(-(1/2)*m-(1/2)*n+1)*m^(-(1/2)*m-(1/2)*n+1)*hypergeom([(1/2)*n, (1/2)*m+(1/2)*n], [1+(1/2)*n], -n/(x*m))*GAMMA((1/2)*m+(1/2)*n)*GAMMA(1-(1/2)*n))/(m*GAMMA(1-(1/2)*n)*GAMMA((1/2)*m)*GAMMA((1/2)*n))

 

> CDF(FRatio(4,3),x) assuming x > 0;

((9+12*x)^(7/2)-2187-10206*x-9720*x^2)/(9+12*x)^(7/2)

 

> with(Statistics):
   CDF(FRatio(m,n),x) assuming x > 0, m::posint, n::posint;
 

x^((1/2)*m-1)*(m/n)^((1/2)*m)*n*(GAMMA((1/2)*m)*Pi*csc((1/2)*Pi*n)*x^(-(1/2)*m+1)*m^(-(1/2)*m+1)*n^((1/2)*m-1)-2*n^(-2+(1/2)*m+(1/2)*n)*x^(-(1/2)*m-(1/2)*n+1)*m^(-(1/2)*m-(1/2)*n+1)*hypergeom([(1/2)*n, (1/2)*m+(1/2)*n], [1+(1/2)*n], -n/(x*m))*GAMMA((1/2)*m+(1/2)*n)*GAMMA(1-(1/2)*n))/(m*GAMMA(1-(1/2)*n)*GAMMA((1/2)*m)*GAMMA((1/2)*n))

 

> CDF(FRatio(4,3),x) assuming x > 0;

((9+12*x)^(7/2)-2187-10206*x-9720*x^2)/(9+12*x)^(7/2)

 

Maple V Release 4, the earliest version I have running,  gives the integral a value of -1/2/r*arctan(1/r)+1/8*(I*ln(r^2+1)*(I*r-1)^(1/2)+2*arctan(1/r)*(I*r-1)^(1/2)-I*ln(r^2+9)*(I*r-1)^(1/2)-2*arctan(3/r)*(I*r-1)^(1/2)-4*I*arctan(2/(I*r-1)^(1/2)))/r/(I*r-1)^(1/2)

without any assumptions.  But that's not really progress...

All you need to do nowadays is

> int(1/(y^4-2*y^3+2*I*r*y^2-2*I*r*y-r^2),y=-1..1) assuming r > 0;

-1/8*(ln(r^2+9)*(-1+r*I)^(1/2)*I+2*arctan(3/r)*(-1+r*I)^(1/2)+4*I*arctan(2/(-1+r*I)^(1/2))-I*ln(r^2+1)*(-1+r*I)^(1/2)+2*arctan(1/r)*(-1+r*I)^(1/2))/r/(-1+r*I)^(1/2)

> map(normal, series(%,r)) assuming r > 0;

series((-1/2*Pi)*r^(-1)+(-1/8*ln(3)-1/8*I*Pi+1/2)+(3/32*Pi-3/32*I*ln(3)+5/24*I)*r+(-25/144+5/64*I*Pi+5/64*ln(3))*r^2(-1415/10368*I-35/512*Pi+35/512*I*ln(3))*r^3+(-63/1024*ln(3)-63/1024*I*Pi+12223/103680)*r^4+O(r^5),r,5)

 

A polynomial of degree 3 has only 3 roots.  Given any 6 distinct numbers, there is no polynomial (except the constant 0) of degree less than 6 that they all satisfy.

A polynomial of degree 3 has only 3 roots.  Given any 6 distinct numbers, there is no polynomial (except the constant 0) of degree less than 6 that they all satisfy.

seq(assign(a[n], a[n-1]+1/2^(n-1)), n = 2 .. 10);

will not work if a[2] to a[10] have already been assigned values.  You could start by unassigning a with

a := 'a':
a[1]:= etc

Or you could try

seq(assign('a[n]', a[n-1]+1/2^(n-1)), n = 2 .. 10);

 

seq(assign(a[n], a[n-1]+1/2^(n-1)), n = 2 .. 10);

will not work if a[2] to a[10] have already been assigned values.  You could start by unassigning a with

a := 'a':
a[1]:= etc

Or you could try

seq(assign('a[n]', a[n-1]+1/2^(n-1)), n = 2 .. 10);

 

It's really not hard to get examples.  The basic strategy is: take a random function F of one variable, involving some square roots and some transcendental functions, then try

f := simplify(diff(F,x));

There's a pretty good chance that int(f, x) will come back unevaluated.

For example:

int(1/(ln(x)+2*x)^(3/2)*(1+2*x)/x,x) = -2/(ln(x)+2*x)^(1/2)

Please note: I do not regard this as evidence of a bug in Maple.
It merely shows that the Risch algorithm is not completely
implemented.

It's really not hard to get examples.  The basic strategy is: take a random function F of one variable, involving some square roots and some transcendental functions, then try

f := simplify(diff(F,x));

There's a pretty good chance that int(f, x) will come back unevaluated.

For example:

int(1/(ln(x)+2*x)^(3/2)*(1+2*x)/x,x) = -2/(ln(x)+2*x)^(1/2)

Please note: I do not regard this as evidence of a bug in Maple.
It merely shows that the Risch algorithm is not completely
implemented.

Here's how I might do it.

 

y := 3*cos(4*Pi*x-13/10)+5*cos(2*Pi*x+1/2);

 

A convenient change of variables:

y:= expand(eval(y, x =solve(2*Pi*x+1/2=t,x)));

Express the derivative in terms of s = sin(t) and c = cos(t).

e1:= eval(diff(y, t), {sin(t)=s, cos(t)=c});

Here's the favourite trigonometric identity:

e2:= s^2 + c^2 - 1;

Solve for s and c:

S := solve({e1, e2}, {s,c});

Note that this involves the RootOf a quartic.  Looking at the graph of y, we see that all four roots will be real.

SS := [allvalues(S, implicit)];

Here are the corresponding y values.

yvalues:= map2(eval, eval(y, {sin(t)=s, cos(t)=c}), SS);

Use evalf to get the numerical values.

evalf(yvalues);

[-7.689610910, 5.781103186, 1.279771068, 2.016811697]

So the largest y value is the second.

simplify(yvalues[2]);

 

-1/72*(1032*sin(23/10)*RootOf(120*sin(23/10)*_Z^3+36*sin(23/10)^2-60*sin(23/10)*_Z-119*_Z^2+144*_Z^4,index = 2)^2-366*sin(23/10)-55*RootOf(120*sin(23/10)*_Z^3+36*sin(23/10)^2-60*sin(23/10)*_Z-119*_Z^2+144*_Z^4,index = 2)-540*RootOf(120*sin(23/10)*_Z^3+36*sin(23/10)^2-60*sin(23/10)*_Z-119*_Z^2+144*_Z^4,index = 2)*cos(23/10)^2+720*RootOf(120*sin(23/10)*_Z^3+36*sin(23/10)^2-60*sin(23/10)*_Z-119*_Z^2+144*_Z^4,index = 2)^3)/cos(23/10)/sin(23/10)

Because type is (usually) syntactic rather than semantic, i.e. it looks at what the Maple structure of an object is, rather than the mathematical nature of what it represents.  In this case, your expression is not literally an  integer or  a quotient of two integers, so it is not of type rational.  You can, however, use simplify to transform this expression into 2.

 

First 143 144 145 146 147 148 149 Last Page 145 of 187