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

To save the values of all user-assigned variables, try
>  parse(sprintf("save %q,\"myfile.txt\";", 
         anames(user)),statement);
If you've used names beginning with an underscore or containing forward slashes, try "alluser" instead of "user". See the help page ?anames.
For a function whose body is a single expression
F := (x1, x2) -> expression;
I think codegen[GRADIENT] gives pretty much the same result as diff. The advantage comes in when the procedure has more than one step. However, it may be possible to further improve the code produced by GRADIENT using optimize, which will help by looking for common subexpressions. As far as I know, the "advantage" here is in terms of efficiency, rather than accuracy.
This is a Java version problem. See http://www.mapleprimes.com/forum/delete-key-problems. It should be fixable as outlined in http://docs.info.apple.com/article.html?artnum=303526.
I think you mean something like this.
> him := proc(i)
  local j, z, s, R;
  R:= NULL;
  for j in [k[1],k[2],k[3]] do
   for z from 1 to 3 do
     s:= solve(A[i][3,z] = T[z] - 1, j);
     if type(s, integer) then R:= R,j=s; break end if;
   end do
 end do;
 R
 end proc;
Some points to note: 1) You can build up a sequence by starting with NULL and adding one member at a time. For a very long sequence this wouldn't be very efficient (much better to use seq if possible), but here it's no problem. 2) Don't call "solve" twice with the same input; instead, save its output to a variable and use that. 3) Use "type" rather than "is" to check if something is literally an integer; "is" should be reserved for those cases where some expression might not be a literal integer but might be shown to represent one, e.g.
> is(a^2+1, integer) assuming a::integer;
You could try this:
c2:= subs(map(t -> (t=evalf(t)), 
  indets(convert(c,CompSeq), fraction)), eval(c));
Compiler:-Compile(c2);
As acer says, in this case you can use limit(Opl, hoekvers=2*Pi). There's a mathematical technicality here: it's not always true you can take the pointwise limit of the solution of a DE initial value problem depending on a parameter and get a solution of the limit of the problem. For example, consider { diff(y(t),t) = 2*n*t/(1+n^2*t^4), y(0) = 0 } which has solution y[n](t) = arctan(n*t^2). The limit of the problem as n -> infinity is { diff(y(t),t) = 0, y(0) = 0} whose solution y(t) = 0 is not the limit of y[n](t) as n -> infinity. But in your case the limit works. However, I'm puzzled by the statement "I've tried to work around this problem by substituting 2*Pi before solving the equation but then Maple gives me no output so nothing to plot."
> dsolve({eval(DV,hoekvers=2*Pi), beginvw}, x(t));
x(t) = 25*sin(2*Pi*t)/Pi^2+5*cos(2*Pi*t)-50*cos(2*Pi*t)/Pi*t
If your differential equation is second-order and linear, as these example have been, you know the general solution is of the form f(x) = f[p](x) + c[1]*f[1](x) + c[2]*f[2](x) where f[p] is one particular solution (with whatever boundary conditions you wish), f[1] and f[2] are two linearly independent solutions of the homogeneous equation. If f(0) = b is the other boundary condition, you might take f[p] with f[p](0) = b, f[1] with f[1](0) = 0, so you know c[2] will have to be 0. Then you can determine c[1] by int(f(x),x=0..1) = int(f[p](x),x=0..1) + c[1] *int(f[1](x),x=0..1) For example (to take a less trivial one), {diff(f(x),x$2) + x*f(x) = 1, f(0) = 2, int(f(x),x=0..1) = 3}
> fp:=subs(dsolve({(D@@2)(f)(x)+x*f(x)=1,f(0)=2,D(f)(0)=0},
        numeric,output=listprocedure),f(x));
  f1:=subs(dsolve({(D@@2)(f)(x)+x*f(x)=0,f(0)=0,D(f)(0)=1},
        numeric,output=listprocedure),f(x));
  Jp:= evalf(Int(fp,0..1));
  J1:= evalf(Int(f1,0..1));
  c1:= solve(Jp + c*J1 = 3, c);
  # and the solution is:
  f:= fp + c1*f1;
Would you consider this elegant?
L:= curry(`*`,phi) @ D;
After multiplying by the denominator, you get a cubic equation in a. Although Maple can solve that, the result is not pleasant; moreover, it's not obvious how many of the cubic's roots will be real: it happens to be three. It would be easier to solve for b (since the equation depends only on b^2). Or, since you say you want a plot, you could give the equation (after substituting numerical values for c, d and e) to plots[implicitplot] or algcurves[plot_real_curve]. For example:
> eq := a^2 = (b^2*c^2) + (a*d^2)/(a-e^2);
  p := eval(numer(lhs(eq)-rhs(eq)),{c=1,d=2,e=3});
  algcurves[plot_real_curve](p, a, b);
or
> plots[implicitplot](p,a=-4..12,b=-15..15,grid=[100,100]);
The same example with 500000 points seemed to run out of RAM and start into virtual memory, so it took several minutes to plot, with Bytes showing 354M (I have 1GB of RAM, and Firefox was also running).
Just for the sake of having one in 3D, again with 100000 points:
a:= 1.4; b:= 0.3;
  HH3 := proc(x0,y0,n,A)
  local x,y,X,Y,i;
  X:= x0; Y:= y0;
  for i from 1 to n do
    # calculate the mapping from [X,Y] to [x,y]
    x:= 1 + Y - a*X^2;
    y:= b*X;
    # make the result X and Y for the next iteration
    X:= x;
    Y:= y;
    # and save in the Array
    A[i,1]:= x;
    A[i,2]:= y;
    A[i,3]:= sin(i);
  end do;
  end proc;
 
 A:= Array(1..100000,1..3,datatype=float[8]);
 evalhf(HH3(0.63135448, 0.18940634, 100000, A)):
 plots[pointplot3d](A,  symbol=point,axes=box);
This time computing the Array took 0.313 seconds. The plot took 4.906 seconds of CPU time, but about 15 seconds in real time. This is in Classic 11.01 under Windows XP on a 3 GHz Pentium 4.
For efficiency, it's best to generate your points in an Array with datatype=float[8]. If your procedure is suitable, you might use evalhf to do this. Then you can pass the Array directly to pointplot3d. I don't have a good 3-d attractor handy, but here's one in 2-d: an attractor of the Henon map.
> a:= 1.4; b:= 0.3;
  HH := proc(x0,y0,n,A)
  local x,y,X,Y,i;
  X:= x0; Y:= y0;
  for i from 1 to n do
    # calculate the mapping from [X,Y] to [x,y]
    x:= 1 + Y - a*X^2;
    y:= b*X;
    # make the result X and Y for the next iteration
    X:= x;
    Y:= y;
    # and save in the Array
    A[i,1]:= x;
    A[i,2]:= y;
  end do;
  end proc;
 
 A:= Array(1..100000,1..2,datatype=float[8]);
 evalhf(HH(0.63135448, 0.18940634, 100000, A)):
 plots[pointplot](A,  symbol=point);
Calculating 100000 points took me 0.171 seconds, plotting them took 3.469 seconds.
By "find the offending rows" I assume you mean find a nonzero vector w such that (w^%T).Z1 = 0, where Z1 is your matrix. That means the entries of this vector are coefficients in a linear combination of rows of Z1 that is 0. The way to do this in Maple is
> NullSpace(Z1^%T):
In this case, this produces one vector w whose entries are all 0 except w[9] = -59/108 and w[10] = 1. Thus the 10th row of Z1 is 59/108 times the 9th row.
The main problem is that your g(x,y) works only if x and y are numbers, not symbolic variables; plot3d, like most Maple commands, evaluates its arguments first, so it will call g with symbolic arguments x and y. This is called premature evaluation. One work-around is to change g(x,y) to something that works with symbolic variables, but will evaluate correctly when numbers are substituted for x and y. In this form of plot3d you want the colour argument to be a list of three expressions. So, for example:
> plot3d(1, x = 0 .. 2, y = 0 .. 2, colour = 
[`if`(x < y, 1, 0) $ 3]);
Another way is to use procedures instead of expressions. In this case the first argument to plot3d should be a procedure (but 1 works fine there), the x and y ranges should be given without names, and the colour argument should be a list of three procedures. So:
> g1:= proc(x,y) if x < y then 1 else 0 fi end proc:
  plot3d(1, 0 .. 2, 0 .. 2, colour = [g1, g1, g1]);
But you say you want a 2d plot, not a 3d plot. Now densityplot from the plots package should be what you want. The following should work, for example:
> plots[densityplot](`if`(x < y, 1, 0), x = 0 .. 2, y = 0 .. 2);
Unfortunately it doesn't work as well in Maple 11 as it should, because Maple seems to be trying to "smooth out" the colours. This is a "feature" that I consider a bug.
You might try some of the commands in the gfun package. For example,
> L := [1, 5/2, 19/6, 65/24, 211/120, 133/144, 2059/5040, 1261/8064, 
19171/362880, 2321/145152];
Guess a generating function.
> gfun[guessgf](L,t);
[(-exp(2*t)+exp(3*t))/t, ogf] Write as a series
> convert(%[1],Sum);
Sum((3*3^_k1/(1+_k1)!-2*2^_k1/(1+_k1)!)*t^_k1,_k1 = 0 .. infinity) So term number n in your sequence (starting at 1 rather than 0) is:
> simplify(eval(op(1,%), {_k1=n-1, t=1}));
(3^n-2^n)/n!
First 131 132 133 134 135 136 137 Page 133 of 138