Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

Extending the above procedure to arbitrary integral powers:

CombineLikePowers:=proc(u::algebraic,L::{posint,set(posint),list(posint)}:=infinity) local S0,S,k,u0,u1,u2,u3,u4;
    if type(u,`+`) then
       map(procname,u,L)
    elif type(u,`*`) then
       u0:=evalindets(u,function,freeze);
       if not member(true,map(type,{op(u0)},{name^posint,name^negint})) then return u end if;
       S0:=map(abs,map2(op,2,indets(u0,`^`(anything,integer)))) minus {1};
       if L=infinity then
          S:=S0
       elif type(L,posint) then
          S:={L} intersect S0
       else
          S:=convert(L,set) intersect S0
       end if;
       for k in S do
       u1:=evalindets(u0,{algebraic^k,algebraic^(-k)},x->``(root[k](x,symbolic)));
       u2,u3:=selectremove(type,u1,specfunc(algebraic,``));
       u4:=``(eval(u2,``=(()->args)))^k*u3;
       u0:=evalindets(u4,specfunc({name,`^`},``),expand);
       end do;
       thaw(u0)      
    else
       u
    end if
end proc;

You may try it on something like

w:=expand(randpoly([x,y,sin(x^2)],degree=4,dense)/a^4/b^3);

You could do

CombineLikePowers(w);

or if you only want powers 2 and 3:

CombineLikePowers(w,{2,3});

@Alex Smith You are missing the point, which (as I understod it) was to integrate term by term, which is NOT done when you do

value(Int((sum(x^(2*k+1), k = 1 .. infinity))*sqrt(1-x), x = 0 .. 1));

because

Int((sum(x^(2*k+1), k = 1 .. infinity))*sqrt(1-x), x = 0 .. 1);

is not the integral of an infinite sum as I pointed out in my answer above (In the remark: "This is not really an infinite summation example since Maple first computes the sum").


@Alex Smith You are missing the point, which (as I understod it) was to integrate term by term, which is NOT done when you do

value(Int((sum(x^(2*k+1), k = 1 .. infinity))*sqrt(1-x), x = 0 .. 1));

because

Int((sum(x^(2*k+1), k = 1 .. infinity))*sqrt(1-x), x = 0 .. 1);

is not the integral of an infinite sum as I pointed out in my answer above (In the remark: "This is not really an infinite summation example since Maple first computes the sum").


@hirnyk You are quite right. Something seems to suggest that the exact answer is correct, but the numerical evaluation of it is rather sensitive to the value of Digits. According to Maple the sum

sum(int(x^(2*k+1)*sqrt(1-x), x = 0 .. 1), k = 1 .. infinity);

has the value

(1/8)*sqrt(Pi)*sqrt(2)*MeijerG([[1], [11/4, 13/4]], [[5/2, 2, 1], []], -1)

If you do

evalf(sum(int(x^(2*k+1)*sqrt(1-x), x = 0 .. 1), k = 1 .. infinity),6);

you are using Digits = 6 (regardless of the value of Digits).

Interestingly, you get a better (and roughly correct) result than if you Digits = 10 or 20 (in the latter case the result is 4.0448194072012987198*10^11-3.4793277886925613872*10^10*I, which obviously is far off.

@hirnyk You are quite right. Something seems to suggest that the exact answer is correct, but the numerical evaluation of it is rather sensitive to the value of Digits. According to Maple the sum

sum(int(x^(2*k+1)*sqrt(1-x), x = 0 .. 1), k = 1 .. infinity);

has the value

(1/8)*sqrt(Pi)*sqrt(2)*MeijerG([[1], [11/4, 13/4]], [[5/2, 2, 1], []], -1)

If you do

evalf(sum(int(x^(2*k+1)*sqrt(1-x), x = 0 .. 1), k = 1 .. infinity),6);

you are using Digits = 6 (regardless of the value of Digits).

Interestingly, you get a better (and roughly correct) result than if you Digits = 10 or 20 (in the latter case the result is 4.0448194072012987198*10^11-3.4793277886925613872*10^10*I, which obviously is far off.

@hirnyk Thanks for the example. Luckily I titled my answer "A quick fix" and ended with "I am not necessarily saying that it would be a good idea in general, though."

My point was to call attention to the last line in the procedure 'discont', which in its entirety reads like this:

showstat(discont);

discont := proc(f::algebraic, x::name)
local disr, disc;
   1   if nargs <> 2 then
   2     error "exactly 2 arguments expected (f::algebraic,x::name)"
       end if;
   3   _EnvAllSolutions := true;
   4   disr := traperror(`discont/discontR`(f,x));
   5   if disr = lasterror then
   6     error "cannot determine discontinuities"
       end if;
   7   disc := traperror(`discont/discontC`(f,x));
   8   if disc = lasterror then
   9     error "cannot determine discontinuities"
       end if;
  10   disc := `discont/RootOf`(`discont/RemoveIm`(`discont/duplicates`(`union`(disr,disc))));
  11   remove(type,disc,nonreal)
end proc

Another quick fix (definitely not a good idea in general) is to replace line 3 with _EnvAllSolutions := false;

which can be done by

Discont:=subs(true=false,eval(discont));

This will handle both f and g, but would destroy discont(1/sin(x),x);

@hirnyk Thanks for the example. Luckily I titled my answer "A quick fix" and ended with "I am not necessarily saying that it would be a good idea in general, though."

My point was to call attention to the last line in the procedure 'discont', which in its entirety reads like this:

showstat(discont);

discont := proc(f::algebraic, x::name)
local disr, disc;
   1   if nargs <> 2 then
   2     error "exactly 2 arguments expected (f::algebraic,x::name)"
       end if;
   3   _EnvAllSolutions := true;
   4   disr := traperror(`discont/discontR`(f,x));
   5   if disr = lasterror then
   6     error "cannot determine discontinuities"
       end if;
   7   disc := traperror(`discont/discontC`(f,x));
   8   if disc = lasterror then
   9     error "cannot determine discontinuities"
       end if;
  10   disc := `discont/RootOf`(`discont/RemoveIm`(`discont/duplicates`(`union`(disr,disc))));
  11   remove(type,disc,nonreal)
end proc

Another quick fix (definitely not a good idea in general) is to replace line 3 with _EnvAllSolutions := false;

which can be done by

Discont:=subs(true=false,eval(discont));

This will handle both f and g, but would destroy discont(1/sin(x),x);

@Ratch The Java based interface came with Maple 9. And I believe the 2D-input was introduced at the same time (not earlier for sure). Then Maple 9.5, 10, 11, 12, 13, and now 14. So that amounts to 7 levels.

@Ratch The Java based interface came with Maple 9. And I believe the 2D-input was introduced at the same time (not earlier for sure). Then Maple 9.5, 10, 11, 12, 13, and now 14. So that amounts to 7 levels.

@acer By clearing the remember table of `is/internal` and also removing option remember, FAIL is the result when Digits is raised to 300, but not with Digits = 30. But the procedure still remembers!

If only the remember table is removed the answers remain true.

restart;
z:=3.1415926535897932384626433:
is(Pi<z);
                              true
lprint(op(4,eval(`is/internal`)));
table( [( Pi < 3.1415926535897932384626433 ) = true, ( 0., RealRange(Open(-3.1415926535897932384626433+Pi), infinity) ) = true ] )
forget(`is/internal`);
op(4,eval(`is/internal`));
`is/internal`:=subsop(3=NULL,eval(`is/internal`)):
Digits:=30:
is(Pi<z);
                              true
Digits:=300:
is(Pi<z);
                              FAIL
lprint(op(4,eval(`is/internal`)));
table( [( Pi < 3.1415926535897932384626433 ) = FAIL, ( 0., RealRange(Open(-3.1415926535897932384626433+Pi), infinity) ) = FAIL ] )

@acer By clearing the remember table of `is/internal` and also removing option remember, FAIL is the result when Digits is raised to 300, but not with Digits = 30. But the procedure still remembers!

If only the remember table is removed the answers remain true.

restart;
z:=3.1415926535897932384626433:
is(Pi<z);
                              true
lprint(op(4,eval(`is/internal`)));
table( [( Pi < 3.1415926535897932384626433 ) = true, ( 0., RealRange(Open(-3.1415926535897932384626433+Pi), infinity) ) = true ] )
forget(`is/internal`);
op(4,eval(`is/internal`));
`is/internal`:=subsop(3=NULL,eval(`is/internal`)):
Digits:=30:
is(Pi<z);
                              true
Digits:=300:
is(Pi<z);
                              FAIL
lprint(op(4,eval(`is/internal`)));
table( [( Pi < 3.1415926535897932384626433 ) = FAIL, ( 0., RealRange(Open(-3.1415926535897932384626433+Pi), infinity) ) = FAIL ] )

There is still an rtable inside the procedure that you want to translate.

And also convert is not recognized:

doesnotworkeither := proc(p) convert(p,string) end proc;

CodeGeneration[C](doesnotworkeither);
Warning, the function names {convert} are not recognized in the target language
double doesnotworkeither (double p)
{
  return(convert(p, string));
}

There is still an rtable inside the procedure that you want to translate.

And also convert is not recognized:

doesnotworkeither := proc(p) convert(p,string) end proc;

CodeGeneration[C](doesnotworkeither);
Warning, the function names {convert} are not recognized in the target language
double doesnotworkeither (double p)
{
  return(convert(p, string));
}

@Joe Riel My last solution doesn't cover cases like

aa:=a:
plot(aa*sin(x),x=0..b) assuming a=-2, b=7;

which clearly is not acceptable.

So I returned to the first version handling temporary assignments by actually assigning the parameters before evaluating the first argument to `assuming`.

However, for the first version to handle vectors and matrices  op(eval(res,par)) in the last line before 'end proc' replaces op(eval(res)). The eval was inserted to get procedural output (e.g. from dsolve/numeric) evaluated.

Probably more problems will turn up.

p:=subsop(3=overload,eval(`assuming`)):
`assuming`:=overload(
       [
          proc(x::uneval,par1::list({equation,set(equation),list(equation)})) option overload;
           local par,res;
           par:=ListTools:-Flatten(evalindets(par1,set,[op]));
           assign(par);
           try
            res:=eval(x);
           catch:
            error;
           finally  
            map(unassign@lhs,eval(par,2));
           end try;
           op(eval(res,par))
          end proc,
          eval(p)
        ]
):

@Joe Riel The following doesn't split mixed assumptions and temporary assignments, but does handle both classes.

restart;
p:=subsop(3=overload,eval(`assuming`)):
`assuming`:=overload(
     [
        proc(x::uneval,a::list({equation,set(equation),list(equation)}),$) option overload;
            op(eval(x,ListTools:-Flatten(evalindets([a], set, [op]))));
        end proc,
        eval(p)
      ]
):
plot(a*sin(x),x=0..b) assuming {a=-2,b=7};
param := {a=2,b=7}:
plot(a*sin(x),x=0..b) assuming param;
plot(a*sin(x),x=0..b) assuming a=-2,b=7;
plot(a*sin(x),x=0..b) assuming [a=-2],{b=Pi};
simplify(sqrt((x*y)^2)) assuming x>0,y>0;
int(exp(a*t),t=0..infinity) assuming a<0;

#Interestingly, this one works (with or without overloading), but in the result b is not replaced by 2:

int(exp(b*a*t),t=0..infinity) assuming b=2,a<0;

#This one doesn't work, because the main procedure doesn't accept sets or lists of properties.

int(exp(b*a*t),t=0..infinity) assuming {a<0,b=2};

First 212 213 214 215 216 217 218 Last Page 214 of 223