Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@hirnyk 

This is in Maple 14. Comments about 13 at the end.

restart;
dlist := [0., .1666666667*Pi, .3333333333*Pi, .5000000000*Pi, .6666666667*Pi, .8333333333*Pi, 1.000000000*Pi, 1.166666667*Pi, 1.333333333*Pi, 1.500000000*Pi, 1.666666667*Pi, 1.833333333*Pi, 2.000000000*Pi]:
#convert~(dlist,rational) gives
#[0, (1/6)*Pi, (1/3)*Pi, (1/2)*Pi, (2/3)*Pi, (5/6)*Pi, Pi, (7/6)*Pi, (4/3)*Pi, (3/2)*Pi, (5/3)*Pi, (11/6)*Pi, 2*Pi]

#The following doesn't take much time:

S:=[seq(limit(csc(6*x),x=dlist[k]),k=1..nops(dlist))];
# There are 3 kinds of results, which are shown below.
Sn:=ListTools:-Enumerate(S):
#For these values of k the results are as they ought to be. Not surprising, since they
#correspond to integral multiples of Pi/2, so that no roundoff error is possible:
op~(1,select(has,Sn,Float(undefined)));
                       [1, 4, 7, 10, 13]
#For these values of k the limits are unevaluated:
op~(1,select(has,Sn,limit));
                      [2, 3, 5, 6, 11, 12]
#For the remaining two values of k the limits are evaluated to large absolute values:
op~(1,remove(has,Sn,{Float(undefined),limit}));
                             [8, 9]

#Now with these results in mind the following results are not surprising
seq(time(evalf(S[k])),k=1..nops(dlist));
 0.015, 17.472, 17.550, 0., 17.674, 17.518, 0.015, 0., 0., 0., 17.768, 17.737, 0.

The unevaluated limits take about the same amount of time.

So what remains to explain is, why does limit succeed on its own for k = 8 and k = 9, i.e. corresponding to dlist-values 1.166666667*Pi and 1.333333333*Pi?

In Maple 13 the first results of limit (S above) are the same as in Maple 14. But using evalf on S doesn't change S at all, i.e. evalf/limit gives up and does it quickly.

In Maple 14 seq(evalf(S[k]),k=1..nops(dlist)); actually evaluates the previously unevaluated limits and that takes up the time. Maybe it ought to have given up. 

@hirnyk 

This is in Maple 14. Comments about 13 at the end.

restart;
dlist := [0., .1666666667*Pi, .3333333333*Pi, .5000000000*Pi, .6666666667*Pi, .8333333333*Pi, 1.000000000*Pi, 1.166666667*Pi, 1.333333333*Pi, 1.500000000*Pi, 1.666666667*Pi, 1.833333333*Pi, 2.000000000*Pi]:
#convert~(dlist,rational) gives
#[0, (1/6)*Pi, (1/3)*Pi, (1/2)*Pi, (2/3)*Pi, (5/6)*Pi, Pi, (7/6)*Pi, (4/3)*Pi, (3/2)*Pi, (5/3)*Pi, (11/6)*Pi, 2*Pi]

#The following doesn't take much time:

S:=[seq(limit(csc(6*x),x=dlist[k]),k=1..nops(dlist))];
# There are 3 kinds of results, which are shown below.
Sn:=ListTools:-Enumerate(S):
#For these values of k the results are as they ought to be. Not surprising, since they
#correspond to integral multiples of Pi/2, so that no roundoff error is possible:
op~(1,select(has,Sn,Float(undefined)));
                       [1, 4, 7, 10, 13]
#For these values of k the limits are unevaluated:
op~(1,select(has,Sn,limit));
                      [2, 3, 5, 6, 11, 12]
#For the remaining two values of k the limits are evaluated to large absolute values:
op~(1,remove(has,Sn,{Float(undefined),limit}));
                             [8, 9]

#Now with these results in mind the following results are not surprising
seq(time(evalf(S[k])),k=1..nops(dlist));
 0.015, 17.472, 17.550, 0., 17.674, 17.518, 0.015, 0., 0., 0., 17.768, 17.737, 0.

The unevaluated limits take about the same amount of time.

So what remains to explain is, why does limit succeed on its own for k = 8 and k = 9, i.e. corresponding to dlist-values 1.166666667*Pi and 1.333333333*Pi?

In Maple 13 the first results of limit (S above) are the same as in Maple 14. But using evalf on S doesn't change S at all, i.e. evalf/limit gives up and does it quickly.

In Maple 14 seq(evalf(S[k]),k=1..nops(dlist)); actually evaluates the previously unevaluated limits and that takes up the time. Maybe it ought to have given up. 

Digging into where the procedure discont seems to take up the time, it appears that

`plot/discontplot`:-ModuleApply

calls the procedure discontlist local to `plot/discontplot`

which in turn calls another procedure alldisconts local to `plot/discontplot`.

Now `plot/discontplot`::alldisconts contains a procedure g which is mapped onto a list of problematic points called dlist.

So what seems to be happening is this:

restart;
dlist := [0., .1666666667*Pi, .3333333333*Pi, .5000000000*Pi, .6666666667*Pi, .8333333333*Pi, 1.000000000*Pi, 1.166666667*Pi, 1.333333333*Pi, 1.500000000*Pi, 1.666666667*Pi, 1.833333333*Pi, 2.000000000*Pi]:
g:=proc(z) local t,s;
   try
      t:=limit(csc(6*x),x=z);
      if type(evalf(t),'numeric') then
         try
            s:=eval(csc(6*x),x=z);
            if s=t then return NULL end if
         catch:  
         end try;
         return [z,t] else
         return NULL
      end if
   catch:
      return NULL
   end try
end proc:
time(map(g,dlist));
                            105.222

#It really is evalf on the unevaluted limits that takes up the time:

restart;
dlist := [0., .1666666667*Pi, .3333333333*Pi, .5000000000*Pi, .6666666667*Pi, .8333333333*Pi, 1.000000000*Pi, 1.166666667*Pi, 1.333333333*Pi, 1.500000000*Pi, 1.666666667*Pi, 1.833333333*Pi, 2.000000000*Pi]:
time(evalf(seq(limit(csc(6*x),x=dlist[k]),k=1..nops(dlist))));
                            105.285

In Maple 13:

restart;
dlist := [0., .1666666667*Pi, .3333333333*Pi, .5000000000*Pi, .6666666667*Pi, .8333333333*Pi, 1.000000000*Pi, 1.166666667*Pi, 1.333333333*Pi, 1.500000000*Pi, 1.666666667*Pi, 1.833333333*Pi, 2.000000000*Pi]:
time(evalf(seq(limit(csc(6*x),x=dlist[k]),k=1..nops(dlist))));
                                    0.748

Digging into where the procedure discont seems to take up the time, it appears that

`plot/discontplot`:-ModuleApply

calls the procedure discontlist local to `plot/discontplot`

which in turn calls another procedure alldisconts local to `plot/discontplot`.

Now `plot/discontplot`::alldisconts contains a procedure g which is mapped onto a list of problematic points called dlist.

So what seems to be happening is this:

restart;
dlist := [0., .1666666667*Pi, .3333333333*Pi, .5000000000*Pi, .6666666667*Pi, .8333333333*Pi, 1.000000000*Pi, 1.166666667*Pi, 1.333333333*Pi, 1.500000000*Pi, 1.666666667*Pi, 1.833333333*Pi, 2.000000000*Pi]:
g:=proc(z) local t,s;
   try
      t:=limit(csc(6*x),x=z);
      if type(evalf(t),'numeric') then
         try
            s:=eval(csc(6*x),x=z);
            if s=t then return NULL end if
         catch:  
         end try;
         return [z,t] else
         return NULL
      end if
   catch:
      return NULL
   end try
end proc:
time(map(g,dlist));
                            105.222

#It really is evalf on the unevaluted limits that takes up the time:

restart;
dlist := [0., .1666666667*Pi, .3333333333*Pi, .5000000000*Pi, .6666666667*Pi, .8333333333*Pi, 1.000000000*Pi, 1.166666667*Pi, 1.333333333*Pi, 1.500000000*Pi, 1.666666667*Pi, 1.833333333*Pi, 2.000000000*Pi]:
time(evalf(seq(limit(csc(6*x),x=dlist[k]),k=1..nops(dlist))));
                            105.285

In Maple 13:

restart;
dlist := [0., .1666666667*Pi, .3333333333*Pi, .5000000000*Pi, .6666666667*Pi, .8333333333*Pi, 1.000000000*Pi, 1.166666667*Pi, 1.333333333*Pi, 1.500000000*Pi, 1.666666667*Pi, 1.833333333*Pi, 2.000000000*Pi]:
time(evalf(seq(limit(csc(6*x),x=dlist[k]),k=1..nops(dlist))));
                                    0.748

@pagan If a function is not defined at a point, the question of its possible continuity or differentiability at that point cannot come up.

Notice also that the help page for isdifferentiable specifically refers to piecewise defined functions.
?isdifferentiable
f:=piecewise(x=3,7,1/(x-3)):
             
isdifferentiable(f,x,1);
                             false
isdifferentiable(1/(x-3),x,1);
                              true


@pagan If a function is not defined at a point, the question of its possible continuity or differentiability at that point cannot come up.

Notice also that the help page for isdifferentiable specifically refers to piecewise defined functions.
?isdifferentiable
f:=piecewise(x=3,7,1/(x-3)):
             
isdifferentiable(f,x,1);
                             false
isdifferentiable(1/(x-3),x,1);
                              true


@snackman 

Integrate eqn1 from s = 0 to s = 1.5:

int(15*(diff(theta(s), s, s)), s = 0 .. 1.5) = int(N3*cos(theta(s))-N1*sin(theta(s)), s = 0 .. 1.5);

Now the requirements eqn4, eqn5 and bcns imply that the right hand side is just N3.

Thus 15*(D(theta)(1.5)-D(theta)(0)) = N3.

But from the equation for diff(theta(s),s) we find (taking signs into account!):

D(theta)(0) = - sqrt( (2/15)*N1 + C)

D(theta)(1.5) = sqrt( (2/15)*N1 + C)

By combining these results, C is easily found.

@snackman 

Integrate eqn1 from s = 0 to s = 1.5:

int(15*(diff(theta(s), s, s)), s = 0 .. 1.5) = int(N3*cos(theta(s))-N1*sin(theta(s)), s = 0 .. 1.5);

Now the requirements eqn4, eqn5 and bcns imply that the right hand side is just N3.

Thus 15*(D(theta)(1.5)-D(theta)(0)) = N3.

But from the equation for diff(theta(s),s) we find (taking signs into account!):

D(theta)(0) = - sqrt( (2/15)*N1 + C)

D(theta)(1.5) = sqrt( (2/15)*N1 + C)

By combining these results, C is easily found.

A minor modification (admittedly making the technique more complicated)  makes it possible to avoid determining the points of intersection:

plottools[transform]((x,y) -> [x,y + x + 5])(plot(min(x^2-x-5,0),x=-3..4,filled=true));

A minor modification (admittedly making the technique more complicated)  makes it possible to avoid determining the points of intersection:

plottools[transform]((x,y) -> [x,y + x + 5])(plot(min(x^2-x-5,0),x=-3..4,filled=true));

Your integrand is not only oscillating, it is oscillating wildly.

Even on a much smaller interval you will have problems.

This works though:

evalf(Int(sin(x^19),x=0..2,method = _d01akc,maxintervals=10000));
                         0.08028631940

See ?evalf,Int

Preben Alsholm

I don't know what you are trying to do, but I would start like this:

restart;
C:=Vector(3):
A := evalf(Matrix([[2, -1, sqrt(2)], [3, 2, -3], [3, sqrt(2), -15/7]]));
B := evalf(<5+7*sqrt(2), -24, -12-3*sqrt(2)>);
N := 3: t := 1e-6:
P := evalf((1/2)*sqrt(t*N)):

You don't need the packages linalg and student, and they are superseded anyway by the LinearAlgebra and the Student packages.

Preben Alsholm

I don't know what you are trying to do, but I would start like this:

restart;
C:=Vector(3):
A := evalf(Matrix([[2, -1, sqrt(2)], [3, 2, -3], [3, sqrt(2), -15/7]]));
B := evalf(<5+7*sqrt(2), -24, -12-3*sqrt(2)>);
N := 3: t := 1e-6:
P := evalf((1/2)*sqrt(t*N)):

You don't need the packages linalg and student, and they are superseded anyway by the LinearAlgebra and the Student packages.

Preben Alsholm

The removal of evalf in my previous comment could be replaced by the version below, which seems to take care of the title problem.

The first change to `plots/animate` is the same as before.

`plots/animate`:=subs(subs=((x,y)->eval(y,x)),eval(`plots/animate`)):
`plots/animate`:=subs(evalf=(proc(x::uneval) local a;a:=[eval(x)];if _npassed>1 then evalf(a[1],_rest) else op(a) end if;end proc), eval(`plots/animate`)):
`plots/animate`:=subs(numeric=realcons,eval(`plots/animate`)):


Preben Alsholm

The removal of evalf in my previous comment could be replaced by the version below, which seems to take care of the title problem.

The first change to `plots/animate` is the same as before.

`plots/animate`:=subs(subs=((x,y)->eval(y,x)),eval(`plots/animate`)):
`plots/animate`:=subs(evalf=(proc(x::uneval) local a;a:=[eval(x)];if _npassed>1 then evalf(a[1],_rest) else op(a) end if;end proc), eval(`plots/animate`)):
`plots/animate`:=subs(numeric=realcons,eval(`plots/animate`)):


Preben Alsholm

First 222 223 224 225 226 227 228 Page 224 of 231