vv

12453 Reputation

19 Badges

9 years, 282 days

MaplePrimes Activity


These are answers submitted by vv

There are several facts to consider.

Not all special functions are standardized. Different authors / software may use distinct conventions (notations).

 

In Maple, EllipticPi(nu, k)  i.e. the complete elliptic integral of the third kind  is defined by the integral

 

Ell := Int(1/((-nu*t^2 + 1)*sqrt(-t^2 + 1)*sqrt(-k^2*t^2 + 1)), t = 0 .. 1)

Int(1/((-nu*t^2+1)*(-t^2+1)^(1/2)*(-k^2*t^2+1)^(1/2)), t = 0 .. 1)

(1)

(Principal Cauchy value (CPV) for nu>1, not mentioned in the help file)

EllipticPi(4,1/3);evalf(%);

EllipticPi(4, 1/3)

 

-0.2328295621e-1-.9197641499*I

(2)

but Mathematica (or WolframAlfa)  uses  k^2  instead of k.

 

So, for EllipticPi(4, 1/3)

WolframAlfa gives   -0.0810909715...

 

but for EllipticPi(4, 1/9)   it gives gives   -0.0232829562...

 

As you see, now the result is the same as in Maple's real part.

 

Let's compute the CPV for Ell in this case. We must do some manipulations because in general the numerical integrators in a CAS do not like CPV.
(Maple knows CPV for symbolic computations only).

 

a:=eval(Ell,[nu=4,k=1/3]);

Int(1/((-4*t^2+1)*(-t^2+1)^(1/2)*(-(1/9)*t^2+1)^(1/2)), t = 0 .. 1)

(3)

evalf(a);     # no CPV  
evalf(Im(a)); # obviously, the integrand is real

Float(undefined)

 

0.

(4)

f:=unapply(op(1,a), t);

proc (t) options operator, arrow; 1/((-4*t^2+1)*(1-t^2)^(1/2)*(-(1/9)*t^2+1)^(1/2)) end proc

(5)

ff:=simplify(f(t)+f(1-t)) assuming t>0,t<1;

-3/((4*t^2-1)*(-t^2+1)^(1/2)*(-t^2+9)^(1/2))-3/(t^(1/2)*(-t+2)^(1/2)*(-t^2+2*t+8)^(1/2)*(4*t^2-8*t+3))

(6)

evalf[15](Int(ff, t=0..1/2));

-0.232829562094268e-1

(7)

So, the imaginary part is 0.

The presence of the imaginary part seems to be a bug; maybe the Maplesoft programmer simply considered that it is not important (?).

Yes, `evalf/EllipticCPi`(nu,k) is buggy. E.g.

evalf(EllipticCPi(0.985, 0.4));                 # crash
plot(EllipticCPi(n, 0.91), n=0.97..0.99); # discontinuous

Obvious workaround: use  EllipticPi(nu, sqrt(-k^2 + 1))
You may redefine

`evalf/EllipticCPi` := (nu,k) -> `evalf/EllipticPi`(nu, sqrt(1-k^2));
f:=CurveFitting:-Spline(my_array, x, degree=1):
b:=0.7: a:=solve(f=b,x):
linevert:=plot([a,y,y=0.1..0.9]):
display(my_graph, line, linevert);

restart;
NumPartStrict:=proc(n::integer[8], m::integer[8])::integer[8];
  option autocompile;
  local k;
  if m=0 then return `if`(n=0,1,0) fi;
  if m>=n then return 1 + thisproc(n,n-1) fi;
  add(thisproc(n-m+k-1,m-k),k=1..m)
end:

NumPartStrict(15,8);
                               13
NumPartStrict(125,75);
                            3179161
NumPartStrict(150,50);
                            14682366
 

You want the complex integral, not the line integral.

restart;
f:=exp(I*z)/z; Z:=r*exp(I*t):
J:=Int(eval(f, z=Z) * diff(Z,t), t=0..Pi);
value(J) assuming r>0;
evalf(eval(%,r=1/2)) = evalf(eval(J,r=1/2)); #check

 

For a non-polynomial expression we must define a custom order and use an inert version. For this example:

f:=x^2+x^y+1-z:
`%+`(sort([op(f)], (u,v) -> y in indets(u))[]);
#        x^y - z + 1 + x^2
restart;
ex:=int(diff(u(x,y,t),x)*v(x,y,t)+diff(v(x,y,t),x)*u(x,y,t),x): 
simplify(eval(ex, v=F/u)): eval(%, F=u*v);
#                     u(x, y, t) v(x, y, t)

 

It is always better to use explicit (or parametric) plots when possible.

plots:-animate(plot, [[sqrt(x^3/(2*a - x)), -sqrt(x^3/(2*a - x))], x = -5 .. 5, view=-5..5, color=red], a=-5..5);

If you really want implicitplot, just add the option signchange = false.

You could use inert Sum instead of sum and then, when needed do e.g.:
value(eval(a, N=4));

This is obviously a cooked up example.
The minimum is 0, but it is attained not only for x=y=z=1,  but also for  any x=1, y=a, z=a  (a>=1).

int(sin(log(x+1)),x=-1..1, method=FTOC);

Should work in Maple 17 (the function has an elementary antiderivative).

So, you have

restart;
alias(a1 = a1(r), a2 = a2(r), a3 = a3(r));
#                           a1, a2, a3
array1 := [a1, a2, a3];
#                     array1 := [a1, a2, a3]
array2 := [seq(cat(a, i), i = 1 .. 3)];
#                    array2 := [a1, a2, a3]
array1 - array2;
#                             [0,0,0 ]
diff(array1, r),   diff(array2, r);
#              [ d       d       d    ]           
#              [--- a1, --- a2, --- a3], [0, 0, 0]
#              [ dr      dr      dr   ]           

To understand (and see) what happens, let's use macro instead of alias
 

restart;
macro(a1 = a1(r), a2 = a2(r), a3 = a3(r));
#                           a1, a2, a3
array1 := [a1, a2, a3];
#                array1 := [a1(r), a2(r), a3(r)]
array2 := [seq(cat(a, i), i = 1 .. 3)];
#                     array2 := [a1, a2, a3]
array1 - array2;
#            [-a1 + a1(r), -a2 + a2(r), -a3 + a3(r)]
diff(array1, r),  diff(array2, r);
#          [ d          d          d       ]           
#          [--- a1(r), --- a2(r), --- a3(r)], [0, 0, 0]
#          [ dr         dr         dr      ]           

So, as we see, array1 and array2 are totally distinct, array2 does not contain arguments.
The fact that array1 - array2  evaluates to  [0,0,0]  (for alias only!) is (most probably) due to (automatic?) simplification where the alias acts.

In general, when using alias (or macro), it is not wise to use "special" constructs such as your cat (involving aliased names). But at least, with macro, we see how Maple interprets our expressions.
 

An alternative for Carl's solution.

delind:=proc(T::table, S::set, M::set(list):={}) 
  local i,X; 
  unassign(subs(X=T, [seq(X[i],i=S),seq(X[i[]],i=M)])[])
end proc:

T:=table([11,22,33,44,55,66]): T[a,b]:=c: eval(T);
#           TABLE([1 = 11, 2 = 22, 3 = 33, 4 = 44, 5 = 55, 6 = 66, (a, b) = c ])

delind(T, {2,3});
delind(T, {1,5}, {[a,b],[6]});  eval(T);
#           TABLE([4 = 44])

Edit. @Carl. I think that you prefer:

delind:=proc(T::table, S::set, M::set(list):={}) 
  local X; 
  unassign(subs(X=T, `?[]`~(X,`[]`~(S)) union `?[]`~(X,M))[])
end proc:

 

Statistics is a modern package, full of modules and generated modules, not easy to debug.
For a workaround, as you noticed, it's possible to use (for the moment)

HR:=(X, t) -> (PDF(X, t)/(1 - CDF(X, t)));

HR(1/2*X1 + 1/2*X2, 0.4); 
        2.352941176

Use Tabulate(DF):  instead of  print(Tabulate(DF));
(notice the colon).

First 18 19 20 21 22 23 24 Last Page 20 of 111