Alejandro Jakubi

MaplePrimes Activity


These are answers submitted by Alejandro Jakubi

Look e.g. at this thread on a somewhat similar integral. As Jacques has described there, what is wrong is relying on the FTOC method, even when the indefinite integration stage may produce discontinuous functions and there are known problems for detecting discontinuities and computing one-sided limits.

See also the imaginary part of the primitive function (in blue):

> F:=int(ln(sec(x)+ tan(x)),x);
> plot([Re,Im](F),x=0..10,color=[red,blue]);

It looks like that there is a numerical problem, deciding on which side of the branch cut going along the real axis to evaluate it (choosing the closure). Compare it with the plot produced by the numerically evaluated integral, in the interval (0,2):

J:=u->evalf(Int(ln(sec(x)+ tan(x)),x=0..u)):
plot([Re@J,Im@J],0..2,color=[red,blue]);

which should be equal up to an additive (complex) constant.

Probably, the best workaround is avoiding indexed names. Something like:

> DGsetup([seq(cat(x,i), i = 1 .. 4)], M, verbose);

It happens that the code in the procedure DifferentialGeometry:-DGsetup:-coord_init checks that the independent variables are not indexed names for Maple 17 and 18 (see the elif line):

> stoperror('expected'):
> with(DifferentialGeometry):with(Tools): with(Tensor): with(LinearAlgebra):
> DGsetup([seq(x[i], i = 1 .. 4)], M, verbose);
Error, expected %-1 independent variable to be an unassigned name, received: %2, 1, x[1]
DifferentialGeometry:-DGsetup:-coord_init:
  24       error "expected %-1 independent variable to be an unassigned name, received: %2", i, indep_var_list[i]
DBG> list
DifferentialGeometry:-DGsetup:-coord_init := proc(indep_var_list, opt1, opt2, opt3, {complexconjugatepairs := [], grading := [], jetnotation := ""})
local L, ans1, ans2, ans3, ans4, ans5, coframe_labels_list, prOrder, jetnotation1, prolongationOrder, contact_coframe_labels_list, dep_var_list, frame_labels_list, frame_name, horizontal_coframe_labels_list, i, jet_var_list, n, proc_name, variable_names_list, form_degrees, properties;
       ...
  19     frame_name := opt2;
  20     prolongationOrder := opt3
       end if;
  21   for i to nops(indep_var_list) do
  22     if not type(indep_var_list[i],'name') then
  23       error "expected %-1 independent variable to be an unassigned name, received: %2", i, indep_var_list[i]
         elif indep_var_list[i]::{indexed, 'Matrix', 'Vector', '`module`', 'string', 'function', 'procedure'} then
  24 !     error "expected %-1 independent variable to be an unassigned name, received: %2", i, indep_var_list[i]
         end if
       end do;
  25   for i to nops(dep_var_list) do
         ...
       end do;
       ...
end proc

In Maple 16 there is no such error as the type check is different:

> showstat(DifferentialGeometry:-DGsetup:-coord_init,21..23);
DifferentialGeometry:-DGsetup:-coord_init := proc(indep_var_list::list, {grading := [], jetnotation := ""})
local L, ans1, ans2, ans3, ans4, ans5, coframe_labels_list, prOrder, jetnotation1, prolongationOrder, contact_coframe_labels_list, dep_var_list, frame_labels_list, frame_name, horizontal_coframe_labels_list, i, jet_var_list, n, proc_name, variable_names_list, form_degrees, properties;
       ...
  21   for i to nops(indep_var_list) do
  22     if not type(indep_var_list[i],'name') then
  23       error "invalid arguments; expected independent variable to be a name, not assigned and not a constant, received: %1", indep_var_list[i]
         end if
       end do;
       ...
end proc

Which one is the best? This is a question that I have listened a lot when I was ambassador for Maplesoft in the 90's. As is, it is an ill-defined question. All of these, Maple , Matlab or Mathematica, are "general" systems with (relative) strengths and weaknesses. So, it much depends on what you want to do.

There has been several threads here dealing with comparisons. Some of them:

5th ncrunch math program review

Definite Double Integral = Definite Integral

Why with Mathematica, not with Maple?

Why is Maple last?

maple vs mathematica on wikiVS

maple 16 mathematica 9

continuous or not continuous use style=point

Mathematica potshots continue

Also, for the record, this is an unmaintained page were I have collected (in the 90's) links to sites providing comparisons of CAS, as a partial answer to such questions. Fortunately, a few links still work.

It works directly:

> coeff(p, x^n);
                                         m
                                      3 y
> coeff(coeff(p, x^n), y^m);
                                       3

> coeff(coeff(p, x^(n+2)), y^(m+2));
                                     4 + k

The operation requested by taro and the procedure in acer's answer is basically the same as in several previous threads. Examples:

How to extract a factor from an expression?

Plot of ln(cosh(x)) stops unexpectedly (yank and around)

So, my suggestion is about providing more uniform answers (as oposed to reinventing the wheel every time). In this sense, Preben's procedure in the first linked thread seems general enough to deal with most cases. And when not, improvements might be built upon it. In this case, it works fine:

> FactorOut(expr1,p^(-1/(theta-1)));
            /      1    \ / /    1    \     /  n       /  theta  \\\
            |- ---------| | |---------|     |-----     |---------|||
            \  theta - 1/ | \theta - 1/     | \        \theta - 1/||
           p              |p            p + |  )   q[i]           ||
                          |                 | /                   ||
                          |                 |-----                ||
                          \                 \i = 1                //
> FactorOut(expr1,p^(-1/(theta-1)),simplify);
             /      1    \ / /  theta  \   /  n       /  theta  \\\
             |- ---------| | |---------|   |-----     |---------|||
             \  theta - 1/ | \theta - 1/   | \        \theta - 1/||
            p              |p            + |  )   q[i]           ||
                           |               | /                   ||
                           |               |-----                ||
                           \               \i = 1                //

One way to identify conditions on the parameters is using (differential) elimination tools through PDEtools:-casesplit. The form of the output (how cases are presented) much depends on the choice of the ranking of variables, see ?PDEtools,casesplit. So, this is just one posibility:

> sys:={v*(1-a) = alpha[3], a*u*v*b = alpha[2], v*a*u*(1-b) = alpha[1], (1/2*(-1+a))*a*p*v = alpha[33], 
(-1+a)*(-1+b)*a*p*v*u = alpha[13], (1-a)*a*p*v*u*b = alpha[23], (1/2)*a*u*v*b*(-b*u*p+b*a*u*p+b*g-g) = alpha[22], 
-a*u*v*b*(p*u-u*p*a-g-b*u*p+b*a*u*p+b*g) = alpha[12], (1/2)*a*u*v*(-1+b)*(-b*u*p+b*a*u*p+b*g+p*u-u*p*a) = alpha[11]}:

> PDEtools:-casesplit(sys
,[[a, b, g, p, u]
,[alpha[1], alpha[2], alpha[3], alpha[11], alpha[12], alpha[13], alpha[22]
,alpha[23], alpha[33],v]]); v - alpha[3] alpha[23] [a = ------------, b = ---------------------, g = ( v alpha[13] + alpha[23] 2 4 alpha[13] alpha[22] alpha[33] - alpha[13] alpha[23] 3 + 4 alpha[22] alpha[23] alpha[33] - alpha[23] )/(alpha[13] alpha[23] 2 v alpha[33] (v - alpha[3])), p = - -----------------------, alpha[3] (v - alpha[3]) ... &where [alpha[33] <> 0, v <> 0, -alpha[3] (v - alpha[3]) <> 0, -alpha[13] alpha[23] (v - alpha[3]) <> 0, alpha[13] + alpha[23] <> 0], [ v - alpha[3] 2 v alpha[33] alpha[23] a = ------------, b = 1, p = - -----------------------, u = -1/2 ---------, v alpha[3] (v - alpha[3]) alpha[33] ... alpha[13] = 0, alpha[23] = 0, alpha[33] = 0] &where [ v <> 0, alpha[1] alpha[2] <> 0, alpha[1] + alpha[2] <> 0, alpha[3] - v <> 0 ...

Thus, 13 cases arise. A graphical representation of them as a tree of pair of branches, for each parametric factor equal or different to 0, may be obtained by adding the option caseplot:

                      ========= Pivots Legend =========
                                     p1 = v
                             p2 = u (-v + alpha[3])
                             p3 = alpha[2] alpha[3]
                                 p4 = alpha[23]
                                 p5 = alpha[33]
                                 p6 = alpha[13]
                                 p7 = alpha[1]
                                 p8 = alpha[3]
                             p9 = alpha[2] alpha[1]
                                 p10 = alpha[2]
                              p11 = -v + alpha[3]

It looks like a regression bug introduced in Maple 13. It works as expected in Maple 12.02:

> solve(tmp1);
{eta[p2] = 0.1311308668, eta[p3] = 0.1959062400, eta[p4] = 0.1981391032,
    eta[p5] = 0.2401043444, eta[p6] = 0.8162806605, mu[p] = -2.848639721,
    tau[p3] = 1.841148534, tau[p4] = 2.707274625}

Generically, the best approach is delaying the introduction of floats to the end. Mixing symbols and floats brings frequently a mess both in hand and computer calculations. So, by "veiling" them from your expressions, a non-null result arises:

> with(inttrans):
> 
> vout_fourier_num := fourier(phi[3](t), t, w) = 6.63569999999998*10^(-15)*w^2*fourier(V(t), t, w)/(-5.69875218358308*10^(-40)*w^4+(9.19473390627057*10^(-29)*I)*w^3+2.15219369729956*10^(-18)*w^2-(4.14691648617110*10^(-8)*I)*w-700.8):

> v:=evalindets(vout_fourier_num,float,LargeExpressions:-Veil[V]);
v :=
                                          2
                                         w  fourier(V(t), t, w) V[6]
    fourier(phi[3](t), t, w) = ------------------------------------------------
                                 4              3      2
                               -w  V[3] + V[4] w  I + w  V[5] - V[2] w I - V[1]
> drive:=5.70000000000000*10^(-6)*exp(-3.18877551020408*10^18*(t-2.0*10^(-9))^2)*cos(4.8*10^9*Pi*(t-2.0*10^(-9)))*10^(19/20):

> d:=evalindets(drive,float,LargeExpressions:-Veil[DD]);
                                                                /19\
                                                                |--|
                           2                                    \20/
      d := exp(-(t - DD[3])  DD[2]) cos(Pi (t - DD[3]) DD[1]) 10     DD[4]

> v2:=subs(V(t)=d, v):
> 
> vi:=invfourier(v2, w, t):
> LargeExpressions:-Unveil[DD](vi):
> LargeExpressions:-Unveil[V](%);
                                   19
                                   --
                             -25   20                 -8
phi[3](t) = - 0.1365215404 10    10   (0.1120000000 10   exp(9.600000000 I Pi)
                       -19                 10                       11 2
    exp(0.7840000000 10    (0.4800000000 10   I Pi - 0.1275510204 10  ) )
      (3/2)                   29                  29
    Pi      (-(0.2045857272 10   + 0.1952744220 10   I)
                         10                  11
    exp((-0.1206633323 10   + 0.2027674769 10   I) (t + _U1))
    
...

Yet, this result should be checked (note the names _U1 and _U2, escaped locals probably).

Yes, one way is going into the jet form:

> depfun:=indets(eq[1],And(function,Not(specfunc(anything,'diff'))));
                     depfun := {a[1](x), a[2](x), a[3](x)}
> depvars:=map(x->op([0],x),depfun);
                         depvars := {a[1], a[2], a[3]}
> indepvar:=map(x->op([1],x),depfun);
                                indepvar := {x}
> 
> jeteqn:=PDETools:-ToJet(eq[1],convert(depfun,list), notation = jetnumbers);
jeteqn := a[1][1] a[2][1] + a[1][1] a[2][1, 1] + a[1][1, 1] a[2][1, 1]
     + a[1][1, 1] a[3][1] + a[1][1, 1] + a[3][1]

Then using PDEtools:-Library:-DifforderJet to compute the order of each derivative and adding these orders for each term you get their total order. Here I show a list for each term, back in diff form, together with the term order, so that you can select what you wish:

> map([PDETools:-FromJet,(`+`@op)@(rhs~)@(PDEtools:-Library:-DifforderJet)],
[op(jeteqn)],convert(depvars,list)(op(indepvar))); / 2 \ /d \ /d \ /d \ |d | [[|-- a[1](x)| |-- a[2](x)|, 2], [|-- a[1](x)| |--- a[2](x)|, 3], \dx / \dx / \dx / | 2 | \dx / / 2 \ / 2 \ / 2 \ |d | |d | |d | /d \ [|--- a[1](x)| |--- a[2](x)|, 4], [|--- a[1](x)| |-- a[3](x)|, 3], | 2 | | 2 | | 2 | \dx / \dx / \dx / \dx / 2 d d [--- a[1](x), 2], [-- a[3](x), 1]] 2 dx dx > (x->op(1,x))~(select((x->op(2,x)=2),%)); 2 /d \ /d \ d [|-- a[1](x)| |-- a[2](x)|, --- a[1](x)] \dx / \dx / 2 dx

What you tell sounds to me similar to the issue described in this thread. Are you using the Standard GUI? If so, it may be a problem with the Java sector of the system.

You may use the two-argument command arctan to "wrap" the angle to its principal range:

angle:=x->arctan(sin(x),cos(x));#-Pi..Pi
angle3:=x->arctan(sin(x-Pi),cos(x-Pi))+Pi;#0..2*Pi

plot(angle(x),x=0..4*Pi,discont=true);
plot(angle3(x),x=0..4*Pi,discont=true);

See also this thread.

The question on how to get the result in terms of ln(abs(...)), it is not about "assuming" ln(0)=-infinity (which does not hold), but about expanding ln(z) in its real and imaginary parts ln(z) = ln(|z|) + I*argument(z), see ?ln . It can be done easily with a transformation rule, and then, the "otherwise" case can be simplified with an assumption on the parameter n:

> expandln:=ln(z::algebraic)=ln(abs(z))+I*argument(z):
> int(1/(x+n), x= 10..100);
           {         undefined                And(-100 < n, n < -10)
           {
           { -ln(n + 10) + ln(n + 100)              otherwise

> J:=applyrule(expandln,%);
J := { undefined , And(-100 < n, n < -10)
    -ln(| n + 10 |) - argument(n + 10) I + ln(| n + 100 |)
     + argument(n + 100) I , otherwise

> evalindets(J,argument(algebraic), u->`assuming`([simplify(u)],[n<-100])); 
       {             undefined                    And(-100 < n, n < -10)
       {
       { -ln(| n + 10 |) + ln(| n + 100 |)              otherwise

> evalindets(J,argument(algebraic), u->`assuming`([simplify(u)],[n>-10])); 
       {             undefined                    And(-100 < n, n < -10)
       {
       { -ln(| n + 10 |) + ln(| n + 100 |)              otherwise

I will answer first your original question. The difference in handling between the cases of the integrands 1/(x+n) and 1/x^n, is that, as implemented, discontinuities are searched in the integrand and the primitive functions for the integration variable (here x), not for the parameters (here n). Compare the output produced in these two cases by tracing discont:

> trace(discont):
> int(1/(x+n), x= 10..100);
{--> enter discont, args = 1/(x+n), x
                            _EnvAllSolutions := true
                                  disr := {-n}
                                   disc := {}
                                  disc := {-n}
                                      {-n}
<-- exit discont (now in GetPoles) = {-n}}
{--> enter discont, args = ln(x+n), x
                            _EnvAllSolutions := true
                                  disr := {-n}
                                   disc := {}
                                  disc := {-n}
                                      {-n}
<-- exit discont (now in FindDisconts) = {-n}}
           {         undefined                And(-100 < n, n < -10)
           {
           { -ln(n + 10) + ln(n + 100)              otherwise

> int(1/x^n, x= 10..100);
{--> enter discont, args = 1/(x^n), x
                            _EnvAllSolutions := true
                                  disr := {0}
                                   disc := {}
                                  disc := {0}
                                      {0}
<-- exit discont (now in GetPoles) = {0}}
{--> enter discont, args = -x/(n-1)/(x^n), x
                            _EnvAllSolutions := true
                                  disr := {0}
                                   disc := {}
                                  disc := {0}
                                      {0}
<-- exit discont (now in FindDisconts) = {0}}
                                (-n)   n         (-n)
                         10 (100     10  - 10 100    )
                         -----------------------------
                                     n - 1

So, in the first case a discontinuity is found at x=-n, which may fall within 10<x<100, but in the second case, the discontinuity that is found, x=0, falls outside. Of course, this is a quite different method than looking for the discontinuities of the definite integral for the parameters, as Carl did. And it seems to work in the first case just because of the symmetry between x and n.

On the other hand, the only actual facility implemented in the definite integration sector for handling parameter case splitting is that provided by the option AllSolutions ( see ?int,details ):

The 'AllSolutions' option, if set, forces int to return the entire set of solutions for the specified definite integral. This applies only to parameters in the endpoints of the interval, not to parameters in the integrand.

Consequently, as far as I know, there is nothing else that you could get by making assumptions or defining properties for the parameters.

About Mathematica, this is WolframAlpha result.

First 8 9 10 11 12 13 14 Last Page 10 of 29