Robert Israel

6577 Reputation

21 Badges

18 years, 207 days
University of British Columbia
Associate Professor Emeritus
North York, Ontario, Canada

MaplePrimes Activity


These are answers submitted by Robert Israel

It didn't display because MaplePrimes treats the "lessthan" symbol as the start of an HTML tag. You have to use < ; (without the space) to get <. Assuming you had a linear inequality such as x + y < 3, you can use inequal in the plots package. Thus:
> with(plots):
  inequal(x + y < 3, x = 0 .. 4, y = 0 .. 4);
There are many possible options, for which you should look at the help page ?inequal.
Actually that should be c(t)=1/50, h(t)=11/10, k(t) = -55/7, by my calculations. This equilibrium is a saddle point with two negative eigenvalues and one positive eigenvalue. So there's a two-dimensional surface, the stable manifold, in [c,k,h]-space consisting of points attracted to this saddle point. You said you want k(0) = 55/6, but you have an extra degree of freedom: this should intersect the stable manifold in a curve (if it intersects at all). We might try running the solution backwards, starting at points close to the equilibrium near the stable manifold, and seeing if we can hit k = 55/6. For example, with initial condition c(0) = 1/50, h(0) = 11/10, k(0) = -55/7 + 0.01:
> sys:= {diff(c(t),t) = 2*(c(t))^2 + 21/50*c(t) -  
    1/100*(h(t))^2 - 1/5*c(t)*h(t) - (9/2000)*h(t) + 49/4000,
  diff(h(t),t) = h(t)*(c(t) - 1/5*(h(t)-1)),
  diff(k(t),t) = 8/100*k(t) - h(t) - 1/5*k(t)*h(t)};
  Sol0:= dsolve(sys union 
        {c(0)=1/50,k(0)=-55/7+0.01,h(0)=11/10},
     {c(t),k(t),h(t)},numeric,output=listprocedure);
  odeplot(Sol0,[t,k(t)],t=-100 .. 0, view=[-100..0,
    -55/7 .. 55/6]);
it looks like we hit k = 55/6 at about t = -53. fsolve confirms this:
> ksol:= subs(Sol0,k(t));
  t0:= fsolve(ksol(t)=55/6, t=-60 .. -40);
t0 := -53.14296014 So the initial condition should be the right sides of
> Sol0(t0);
[t(-53.14296014) = -53.14296014, c(t)(-53.14296014) = .200000000000000004e-1, h(t)(-53.14296014) = 1.09999999999999986, k(t)(-53.14296014) = 9.16666667290210490] Hmm, those c and h values look suspicious. In fact it turns out that the straight line c = 1/50, h = 11/10 is a trajectory of the system. Let's try a slightly different starting point.
> Sol1:= dsolve(sys union 
         {c(0)=1/50+0.0002,k(0)=-55/7+0.01,h(0)=11/10},
      {c(t),k(t),h(t)},numeric,output=listprocedure);
   odeplot(Sol1,[t,k(t)],t=-100 .. 0, view=[-100..0,
     -55/7 .. 55/6]);
  ksol:= subs(Sol1,k(t));
  t1:= fsolve(ksol(t)=55/6, t=-60 .. -40);
t1 := -53.37456172 And now the initial condition would be
> [c(0) = subs(Sol1,c(t))(t1),
   h(0) = subs(Sol1,h(t))(t1),
   k(0) = 55/6];
[c(0) = -.177757107426838104e-1, h(0) = .519002763775084696, k(0) = 55/6]
Yes, it does converge, but there's probably no closed form. You could get an asymptotic series for this, using the asymptotic series for 1/sqrt(a^2+w^2). The coefficients are expressed in terms of Whittaker M functions.
You didn't specify an initial condition for your differential equation: I'll suppose that's what the w is for. So perhaps you're looking at the initial value problem diff(y(t),t) = a + b*y(t) + c*y(t)^2 with y(0) = w. Now a suitable transformation of t and y will transform the differential equation into diff(y(t),t) = 1 + y(t)^2 (leaving out exceptional cases). Let's say you're interested in how the solution at t=1 depends on w. This initial value problem actually has a closed form solution y(t)=tan(t+arctan(w)) but let's pretend it doesn't. We can approximate y(t,w) as a power series in w:
> de:= diff(y(t),t) -(1 + y(t)^2);
  des:= {seq(coeff(eval(de, y(t)=add(Y[j](t)*w^j, j=0..20)),
      w,j), j=0..20)};
  ics:= {Y[0](0)=0, Y[1](0)=1, seq(Y[j](0)=0,j=2..20)};
  Sol:= dsolve(des union ics, {seq(Y[j](t),j=0..20)},numeric);
  f1:= subs(Sol(1), add(Y[j](t)*w^j, j=0..20));
Now assume there is one singularity closest to 0 which is a simple pole. Then there should be a good Pade approximant with a denominator of degree 1.
> fpade:= numapprox[pade](f1,w,[19,1]);
fpade := (1.557407728+.9999998218*w-.3636971552e-6*w^2-.4159409454e-6 *w^3-.5525785507e-6*w^4-.7233755573e-6*w^5-.1054922688e-5*w^6 -.1587407473e-5*w^7-.2411251858e-5*w^8-.3616877786e-5*w^9 -.5626254334e-5*w^10-.8238443847e-5*w^11-.1215672812e-4*w^12 -.1758204480e-4*w^13-.2511720685e-4*w^14-.3415940132e-4*w^15 -.4521097233e-4*w^16-.5626254334e-4*w^17-.6329536126e-4*w^18 -.5224379025e-4*w^19)/(1.000000000-1.557407863*w)
> pole:= solve(denom(%),w);
pole := .6420925589 This agrees well with the exact value of cot(1).
I don't recall seeing that feature before. Is it documented in any of the help pages? I couldn't find it.
It's under properties.
Maple can (sometimes) translate from Maple code to Matlab or C, but AFAIK it has no facility to translate from Matlab to anything else.
You have to tell plot that x is the independent variable. So:
> plot(1/7-5*(sin(x)),x = -10..10, -10 .. 10);
You would only omit the "x = " if you gave plot a function rather than an expression to plot. Thus
> plot(1/7 - 5*sin, -10 .. 10, -10 .. 10);
One thing you can do is leave dm undefined but define `expand/dm` (note the back-quotes):
`expand/dm`:= z -> diff(z,t)+v[c]*d_[c](z);
Then dm(rho) will stay as dm(rho). If you want the expanded form of this or any expression containing it, you use expand(...).
The function you use transform on should output a list of three coordinates (if you want to go to 3D). For example, if you want to map a 2D object into the plane y = 3 in 3D, you could use transform((x,z) -> [x, 3, z]).
This is a rather odd bug. It's not a distinction between horizontal and vertical, I think, rather that tickmarks on one axis are specified using strings. If you change ROMAN to ITALIC or BOLD, you'll see that this does affect the horizontal tickmarks as well as the vertical ones. The font size and font family don't affect tickmarks specified using strings. A work-around is to use names rather than strings. Thus
> dates :=  [0 = `21/8/00`, 140 = `8/1/01`, etc ]
ought to do it.
Obviously a bug. The MESH object in your plot3d is in fact MESH("unable to evaluate function `%1` in evalhf", sin). A simpler version of this bug:
> with(RealDomain):
  plot3d([(x,y)->x,(x,y)->y,(x,y) -> sin(x)],0..1,0..1);
Basically the problem is that Maple is trying to evaluate these procedures using evalhf, but RealDomain has its own versions of sin and many other functions, and this confuses evalhf. But I would have thought that plot3d would fall back on using software floats if presented with something for which evalhf didn't work.
If the curves are given parametrically, you can use spacecurve to plot the curves, and display to combine this with a plot of the surface.
Any sequence of plots can be animated. In your example:
> with(plots):
  animate(tubeplot,[[20*t*cos(t+s),20*t*sin(t+s),20*t],
     t=0..20*Pi, radius=-100*exp(.00001*t), tubepoints=20,
     numpoints=100, color=grey], s = 0 .. 2*Pi);
I don't think it's very much like a tornado, but maybe with a bit of imagination...
It's not that Maple is too smart, it's that it isn't smart enough. If you want a loop here, you have to make a loop.
P:= proc(eps,v,n)
  local k, Mr, newMr, Mi;
  Mr:= 1; Mi:= 0;
  for k from 0 to n do
    newMr:= -Mr*(eps + k) - Mi * v;
    Mi:= Mr*v - Mi*(eps + k);
    Mr:= newMr;
  end do;
  Mr
  end proc;

CodeGeneration[C](P, optimize);
First 130 131 132 133 134 135 136 Page 132 of 138