Preben Alsholm

13471 Reputation

22 Badges

20 years, 254 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@Alger When dsolve/numeric talks about a probable singularity one should (as you surely do) notice the adjective "probable".

What exactly is a singularity in this numerical context? What triggers the error (or warning) message?
Your system is very complicated, so it isn't easy to see. It would be nice to find a simpler example with the same problem.

@Alger You get a similar complaint when keeping C at 0.000050 but changing the range from 0..40 to 0..50:

Warning, cannot evaluate the solution further right of 47.288331, probably a singularity.

 

@arashghgood You ask in the title " Is there a solution?".

Clearly without additional requirements there are the very trivial solutions f(x,y,t) = F1(t)*y+F2(t) , (where F1 and F2 are arbitrary).

But I have no idea about how to solve this type of problem numerically.

You could set infolevel[int]:=5:

restart;
integrand1 := (x^3*exp(arcsin(x)))/sqrt(1 - x^2);
integrand2 := (x^3*exp(1)^arcsin(x))/sqrt(1 - x^2);
###
infolevel[int]:=5:
res1:=int(integrand1,x);

res2:=int(integrand2,x)

The first 4 steps are the same, but then they. Obviously the choice of approach is governed by the integrand and they are different.

exp(arcsin(x)) is of type function and exp(1)^arcsin(x) is of type `^` .

@acer Your solution is obviously simpler, but I was about to give this answer:
Wrap the last block of assignment code in a procedure call with SF declared global:

proc() global SF; 
  map(x->assign(SF[x],cat(`SF/`,x)),
  ['('Par')', '('add_basis')', '('char2sf')', '('conjugate')',
   '('dominate')', '('dual_basis')', '('evalsf')', '('hooks')',
   '('itensor')', '('jt_matrix')', '('nextPar')', '('omega')',
   '('plethysm')', '('scalar')', '('sf2char')', '('skew')', '('stdeg')',
   '('subPar')', '('theta')', '('toe')', '('toh')', '('top')',
   '('tos')', '('varset')', '('zee')'])
end proc();


 

@Anthrazit Have you tried my Round procedure?

If so, any comments?

Using an old procedure of mine (chop) I made up the following procedure Round which uses chop.
This combination is very likely unnecessarily clumsy and also not perfect otherwise.

Since, however, my computer is going to a repair shop for booting problems this morning, I dare give the code here with only a few examples.

restart;
chop:=proc (x::anything, n::nonnegint:=0) local op1, op2, k; 
   if not hastype(x,float) then return x end if; 
   if x::float then
     if n=0 then return round(x) end if; 
     op1, op2 := op(x); 
     k := length(op1); 
     evalf[n](round(op1*10^(-k+n))*10^(op2-n+k)) 
   else 
     map(procname,x,n) 
   end if 
end proc:
#####
chop(89.567,0);
chop(sin(1.)+cos(3));
chop(89.567,3);
#####
Round:=proc(x::anything,d::nonnegint:=0) local tr,fr,L,res;
  if not hastype(x,{realcons,complexcons}) then return x end if;
  if x::float or x::integer then 
     tr,fr:=(trunc,frac)(x);
     res:=chop(fr,d)+tr
  elif x::realcons then
     res:=procname(evalf(x),d)
  elif x::complexcons then
     res:=procname(Re(evalf(x)),d)+I*procname(Im(evalf(x)),d)
  else
     res:=evalindets(x,{realcons,complexcons},procname,d)
  end if;
  res
end proc:
#####
Round(4.56789,0);
sqrt(3.);
Round(sqrt(3),4);
Round(4567);
Round(sqrt(2*I),3);
Round([89,sqrt(2)],3);
Round(sin(Pi*x+I*sin(1)),0);
Round(Matrix([[78*Pi,sin(x*exp(1))],[exp(1)*I,erf(4)]]),4);


 

What should acer's Round(123456,2) return in your opinion?

@mmcdara There are in fact no imaginary parts. Those that do turn up after using evalf are artifacts of the floating point computation.

Try evalc(Im(EPS)) and also evalc(Re(EPS)) and you will see that Maple knows that Kummer functions in use are actually real. For perfect zeros in the first case (Im) it is important that the initial values are cleared of the small imaginary "artifacts".

An alternative is to avoid floats as in:

e4,e4_1:=op(eval([EPS1,diff(EPS1,t)],t=4));

and the other ones in the same way. It works, but the drawback is a builtup of huge symbolic expressions.

dsolve does indeed turn floats into rationals:

ode:=diff(y(x),x)=1.23*y(x);
dsolve(ode);

That has always been done in Maple at least as far as back to Maple 8 (the oldest on this computer).
This is definitely done purposely. I vaguely remember a discussion raised by somebody about that on MaplePrimes and that Edgardo Cheb-Terrab defended that handling.

@Carl Love Indeed weird with that number 10.

I do see that `limit/MrvRight` (called several times) in line 11 has the statement:

 if 10 <= nops(`limit/IndetsRange`(t)) then
        return _Range
 end if;

Whether that has anything to do with the importance of 10 I have no clue so far.

If you replace limit by MultiSeries:-limit you get 1.

restart;
I0:=(x - 11)^(n)*x/(x^(n + 1));                                              
I1:=(x - 11)^(n)/(x^(n));                                                    
is(I1=I0);  
Limit(I0,x=infinity)=MultiSeries:-limit(I0,x=infinity) assuming n::posint;
Limit(I1,x=infinity)=MultiSeries:-limit(I1,x=infinity) assuming n::posint;    
### Carl's example:
MultiSeries:-limit((x-a)^n*x/x^(n+1), x= infinity) assuming n::posint, a>10;
MultiSeries:-limit((x-a)^n*x/x^(n+1), x= infinity) assuming n::posint, a<10;

 

It has been noticed and mentioned by some (including very knowlegeable people) that the documentation for dsolve/numeric/events is very difficult to follow unless you in fact know all the terms used in advance. There ought to be many more examples to illustrate all the options mentioned.
As it is now you are basically on your own.
To get any help from this forum you need to be very concrete.

For 2015 I get 8.5 MB and for 2017 I get 7.7 MB using Windows 10.

With Maple 2019.2 I get a disastrous result: The program freezes. I have to kill it. The file was created, but had 0 bytes.
I tried 3 times with Maple 2019.2 with the same disastrous result.

I was using the following code in all 3 cases except that the file name test2019 changed.

with(plots);
plotsetup(postscript, plotoutput = `C:/Users/Bruger/test2019.eps`, plotoptions = `color,portrait,height=350,width=350`);
plot3d(x^2+y^2+1, x = -1 .. 1, y = -1 .. 1, axes = boxed, grid = [200, 200]);

 

Although my patience is not sufficiently large for me to see the message, I think that it is a warning, but not an error. Thus you should in principle have no problem.
I have seen it in other contexts.

@Kitonum The freely available package DirectSearch has the procedure GlobalOptima:
 

DirectSearch:-GlobalOptima(2*x^2 + 2*y^2 + y, {2*x + y <= 6, y^2 - x <= 2},maximize);

It returns: [46.1249999999993, [x = 4.24999999999998, y = -2.49999999999996], 173]

First 25 26 27 28 29 30 31 Last Page 27 of 225