Preben Alsholm

13728 Reputation

22 Badges

20 years, 251 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@laetititia Heaviside(t) = 1 for t > 0 and  Heaviside(t) = 0 for t < 0. In Maple it is left undefined for t = 0, but the user can change that.
Every piecewise defined function can be expressed in terms of translated versions of Heaviside.
In your example:
p:=x->piecewise(x<10,2,x<11,10,2);
convert(p(x),Heaviside);
#Result: -8*Heaviside(-11+x)+8*Heaviside(-10+x)+2
#You can convert this back tp piecewise:
convert(%,piecewise);
#Here you will see some undefined values due to H(0) not being defined.


It looks like a command line plot, but how you got it in the Standard interface I don't know. I have seen this in the past very rarely, but never found out what caused it.
Have you tried restart? Or if that doesn't make the problem go away, try closing Maple and opening it again.

@Alejandro Jakubi Yes, a similar problem also raised by Marko Riedel is found here:

http://www.mapleprimes.com/questions/141286-Computing-A-Residue

Another choice is to change residue only:

Residue:=subs(series=MultiSeries:-series,eval(residue));

or (in this case)
Residue:=subs(normal=simplify@normal,eval(residue));

One more:

Why make the variable change in residue from say s = s0 to x = s-s0 as is done in the code?
Why not just like this:
residue:=proc (f, a::(name = anything)) local g, i, t, x, t1;
   options remember, `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`;
   x := op(1, a);
   if has(op(2, a), x) then
      error "invalid point"
   elif type(op(2, a), 'infinity') then
      g := -(eval(f, x = 1/x))/x^2
   else
      g := f #Changed
   end if;
   if hastype(g, 'RootOf') then g := convert(g, 'radical') end if;
   for i from 0 to 5 do
      try
         t := series(g, a, i^2+2) #Changed
         catch:
         t := FAIL
      end try;
      if t = 0 then return 0 end if;
      if t = FAIL or not type(t, 'laurent') then next end if;
      t1 := order(t);
      if t1 = infinity or -1 < t1 then return coeff(t, x-op(2,a), -1) end if #Change in coeff
   end do;
   'residue(args)'
end proc;

@Alejandro Jakubi I think it is a problem with series, not coeff.
Debugging residue we see that g below appears. g is correct, but the series t for g is wrong:

restart;
g := 1/(2^s*2^((6*I)*Pi/ln(2))-1); #From de
g1:=simplify(g);
t:=series(g,s=0,2); #Wrong
series(g1,s=0,2); #Correct
coeff(%,s,-1);
#MultiSeries:-series does a better job:
MultiSeries:-series(g,s=0,2);
coeff(%,s,-1);
simplify(%);

Converting piecewise to Heaviside helps on the first order equation in u(x) = z'(x).
Since a graph is desired we need initial conditions. Here I take z(0) = 0, z'(0) = 0.

restart;
Eq2 := diff(z(x), x, x) = p(x)*sqrt(1+diff(z(x), x)^2);
p:=x->piecewise(x<10,2,x<11,10,2);
#The following doesn't work with or without conversion
dsolve({convert(Eq2,Heaviside),z(0)=0,D(z)(0)=0}); #error
#The first order equation in u(x) = z'(x):
Eq1:=subs(diff(z(x),x)=u(x),Eq2);
dsolve({convert(Eq1,Heaviside),u(0)=0}); #z'(0) = 0
res:=int(rhs(%),x=0..x)+C;
eval(res,x=0);
sol:=eval(res,C=0); #Since z(0) = 0 is required
plot(ln(sol),x=0..20);
#This compares well with the corresponding numerical solution:
res:=dsolve({Eq2,z(0)=0,D(z)(0)=0},numeric);
plots:-odeplot(res,[x,ln(z(x))],0..20);

Note 1: Using simplify/piecewise as Carl does is nicer I think than convert/Heaviside.
Note 2: However, for some reason the Heaviside version handles my initial conditions z(0)=0,z'(0)=0 better than does the simplify/piecewise version:
Eq1a:=simplify(Eq1,piecewise);
dsolve({Eq1a,u(0)=0}); #Two solutions: One of them wrong
combine(%);
zz:=op(indets(%,`local`));
about(zz);
#This could be handled artificially like this:
eval(dsolve({Eq1a,u(0)=epsilon}),epsilon=0);

@Carl Love Thanks for the warning. Even the title (question) is long!

@trace For whatever reason I cannot see any text in the file resub.pdf.

@trace If the point considered is (x,y,z) = (-1, 0, 0), then what would m(y/z) be, and what would r = z/y be?

So the eigenvalue result you mention doesn't even seem to make sense.

Under the assumption that m has a positive lower bound, i.e. m(s) >= epsilon > 0 it makes some sense to say that (-1,0,0) is an equilibrium point for your system: Taking g(-1,0,0) to mean the limit of g(-1,y,z) as (y,z) -> (0,0), and a similar interpretation for h.
However, if you look at the jacobian at (-1,y,z) then its element
J(-1,y,z)[2,3]
is
-1/m(z/y)+z*(D(m))(z/y)/(m(z/y)^2*y)-2*y
i.e. with r = z/y:
-1/m(r) +r*m'(r)/m(r)^2 -2*y.
But r doesn't have a limit as (y,z) -> 0 so we are stuck even with as simple an example as m(s) = 1+s^2.


As Carl is also pointing out you could replace m in your definitions of f,g,h by m(z/y).
With m a general nonconcrete function solve gives you some RootOf-solutions involving m that seem impossible to do anything about without making m a concrete function, but it also gives you
{x = 1, y = 0, z = 0}, {x = -4, y = 5, z = 0}, {x = -1, y = 0, z = 0}, {x = 0, y = -1, z = 2}
But since in general because of m(z/y) you cannot accept y = 0 you are left with
{x = -4, y = 5, z = 0} and  {x = 0, y = -1, z = 2}.
Applying J to those points you see that you need to know something about m(0) and m(-2).
But the problem remains: What can you say about
RootOf(2*m(3/(2*(_Z-2)))*_Z-4*m(3/(2*(_Z-2)))+2*_Z-1) ?
This represents possible solutions to
2*m(3/(2*(Z-2)))*Z-4*m(3/(2*(Z-2)))+2*Z-1 = 0.
After the change of variable 3/(2*(Z-2)) = s
this means finding solutions to m(s) + s + 1 = 0.
This looks neater, but remember that you have to get back to Z.

@zippon If we define Z[N+1] = Z[1] the relation mentioned above can be written
conjugate(Z[i]) =Z[N-i+2] for i=1..N. (*)

The theorem then is: z[i] are all real iff (*) holds.

And similarly with z[N+1] = z[1] we have:

conjugate(z[i]) =z[N-i+2] for i=1..N iff all Z[i] are real.

 

@zippon If we define Z[N+1] = Z[1] the relation mentioned above can be written
conjugate(Z[i]) =Z[N-i+2] for i=1..N. (*)

The theorem then is: z[i] are all real iff (*) holds.

And similarly with z[N+1] = z[1] we have:

conjugate(z[i]) =z[N-i+2] for i=1..N iff all Z[i] are real.

 

Maple has to isolate the derivative dS/dx. The image seems to show that that would be no problem. So please give us the code either as plain text or as an uploaded worksheet.

@Markiyan Hirnyk Try

addressof~([x-> x^2,x-> x^2]);
#different!
addressof~([x^2,x^2]);
#same!

@Markiyan Hirnyk Try

addressof~([x-> x^2,x-> x^2]);
#different!
addressof~([x^2,x^2]);
#same!

I suppose the problem Maple has with this example is the error handling?

The integrand G

G:=diff(-ln(1-Fy), y)*fy;
# is clearly a smooth function defined on all of R and falls of very rapidly as |y| -> infinity
plot(G,y=-10..10);
asympt(G,y);
#I got the following slightly disturbing error from this command (Digits = 10):
evalf(Int(G,y=-10..10));
# error message:
Error, (in property/LinearProp/+) too many levels of recursion
#However the results of the following were OK:
evalf[15](Int(G, y = -5 .. 5));
evalf[15](Int(G, y = -10 .. 10));
evalf[15](Int(G, y = -20 .. 20));

First 172 173 174 175 176 177 178 Last Page 174 of 230