PatrickT

Dr. Patrick T

2153 Reputation

18 Badges

16 years, 96 days

MaplePrimes Activity


These are replies submitted by PatrickT

The problem will remain even in Maple 17. Now that your problem has been cleaned up, I suggest you do this:

define numerical values of the parameters (all the non-tau symbols), set them equal to something like values you expect them to have, and if you're not sure set them equal to one. Then see if fsolve can solve it. Also worth investigating is whether you can eliminate some equations by eliminating some of the tau variables, if you can reduce the system from 5 to 3 or even to 2, then you can plot it and see visually whether there are any solutions.

I wanted to do the same thing for a plot created with odeplot. The following seems to work:

restart;
### Create plot with CURVES structure (odeplot)

sys := diff(y(x), x) = z(x), diff(z(x), x) = y(x):
var := {y(x), z(x)}:
sol := dsolve({sys, y(0) = 0, z(0) = 1}, var, type = numeric, method = classical):

p := plots:-odeplot(sol, [x, y(x), z(x)]
  , x = -4 .. 4
  , 'style' = point
  , 'numpoints' = 50
  , 'symbol' = diamond
  , 'color' = orange
  , 'labels' = [`x`,`y`,"z"]
  , 'title' = `Amazing Plot`
  , 'axes' = normal
  , 'view' = [-4..4,-20..20,0..30]
  , 'orientation' = [ 45, 45 ]
) :

# plots:-display(%, orientation = [45,45] );

### Procedure Animate Can Export Individual Frames of an Animation
### Works for an odeplot structure

Animate2 := proc(p :: evaln
    , { [n,numframes,frames] :: posint := 4 }
    , { [k,keep] :: boolean := false } # words "export" and "save" are reserved
    , { [x,ext,extension] := png}
    )
    description "based on acer's code 130009-Animation-Of-3D-Plot-Rotation" ;
    local A, o, i, j ;
    o := select(type,[op(eval(p))], specfunc(anything,{ORIENTATION})):
    A := [ seq( PLOT3D( remove(type,[op(eval(p))], specfunc(anything,{ORIENTATION}))[]
                          , ORIENTATION( op([1,1],o)+360/n*(j-1) , op([1,2],o) ) )
               , j = 1 .. n )
         ] ;
    if keep = true then
        try for i from 1 to n do
            plotsetup('x', 'plotoutput'=cat(convert(p,string),convert(i,string),".",convert(x,string)),
                'plotoptions'=`color,portrait,noborder,transparent`) ;
            print( plots:-display( A[i]
                , 'axesfont' = [ TIMES, 10 ]
                , 'labelfont' = [ TIMES, ROMAN, 10]
            ) ) ;
        end do; catch: print("something went wrong"); plotsetup(default); end try;
    end if;
    return plots:-display(A, 'insequence' = true);
end proc:

Animate2(p, 'frames'=100); # fails if there is no orientation option, set up without psi angle in orientation, to be fixed
Animate2(p, 'frames'=4, 'keep'=true);
 

Edit: changed ORIENTATION( 360/n*(j-1) , op([1,2],o) ) to ORIENTATION( op([1,1],o)+360/n*(j-1) , op([1,2],o) ), which is what I intended, to use the initial angle for the first frame

@acer

I've adapted your code to export each individual frame. It sort of works, with minor quirks. It just about does what I need for now, provided I stick to png format, but for the record, let me copy the code here with short comments.

Adapting your code I asked a procedure to use the existing orientation options and rotate according to the first angle. (it currently doesn't work if no orientation has been specified) Then I asked the procedure to export the individual plots for each frame, with a naming scheme like so: nameofplot1, nameofplot2, etc.. The procedure takes the number of frames to be produced as an input, e.g. 'numframes'=4. Then, to export the plot, you specify: 'keep'=true (words "export" and "save" are reserved, so I used "keep"). To export in a certain format, you specify 'ext'=png or 'ext'=ps, with the default set to png.

Exporting as png works on the test plot, named p, and produces correctly rotated p1,p2,p3,p4. Exporting as ps fails in Maple 16 (it works in Maple 15 in the limited sense that it produces the low-quality postscript plots Maple 15 usually produces). A funny thing happens with gif, it produces p1,p2,p3 as expected, but then p4 is an animated gif of all the frames.

Right now I'm using png, so this is good enough. Ultimately it would be nice to be able to export as postscript. A funny thing happened with ps: the first time around I was able to export p1 and p2 in postsript with Maple 16, it took a long time but the exported graphs were very nice. Unfortunately in subsequent attempts it will produce only p1 with an ugly black background and then hangs. Not sure what's going on. I've generally been unable to export postscript plots with Maple 16, so this is not news, but the surprising thing is that I was able to export 2 of the 4 plots at some point during my experiments. Can you export all 4 postscript plots?

 
Edit: Corrected a typo in the code below, where I was incorrectly selecting angles. Also, I think I just noticed that the third angle of ORIENTATION (the psi value) is dropped at times, which I suspect is after a plot device crash.

restart;

### The test plot

p := plot3d(1+x+1*x^2-1*y+1*y^2+1*x*y
  , x = -5 .. 5
  , y = -5 .. 5
  , 'orientation' = [ 45, 45, 45 ]
  , 'axes' = normal
  , 'labels' = [x,y,z]
) :

### Procedure Animate

Animate := proc(p :: evaln
    , { [n,numframes,frames] :: posint := 4 }
    , { [k,keep] :: boolean := false }
    , { [x,ext,extension] := png}
    )
    description "based on acer's code 130009-Animation-Of-3D-Plot-Rotation" ;
    local A, o, i, j ;
    o := select(type,[op(eval(p))], specfunc(anything,{ORIENTATION})):
    A := [ seq( PLOT3D(GRID( op([1,1..2],eval(p)), op([1,3],eval(p)) )
                    , remove(type,[op(eval(p))], specfunc(anything,{GRID,ORIENTATION}))[]
                    , ORIENTATION( 360/n*(j-1) , op([1,2],o) , op([1,3],o) )
                    #, ORIENTATION( 360/n*(j-1) , op([1,2],o) )
              ) , j = 1 .. n )
         ] ;
    if keep = true then
        try for i from 1 to n do
            plotsetup('x', 'plotoutput'=cat(convert(p,string),convert(i,string),".",convert(x,string)),
                'plotoptions'=`color,portrait,noborder,transparent`) ;
            print( plots:-display( A[i]
                , 'axesfont' = [ TIMES, 10 ]
                , 'labelfont' = [ TIMES, ROMAN, 10]
            ) ) ;
        end do; catch: print("something went wrong"); plotsetup(default); end try;
    end if;
    return plots:-display(A, 'insequence' = true);
end proc:

Animate(p, 'frames'=4); # fails if there is no orientation option, to be fixed
Animate(p, 'frames'=4, 'keep'=true); # Maple 15,16: p1,p2,p3 work as expected
                                                 # but p4.gif keeps all the animations together, not sure why
Animate(p, 'frames'=4, 'keep'=true, 'ext'=png); # Maple 15,16: works as expected
Animate(p, 'frames'=4, 'keep'=true, 'ext'=ps); # Maple 16: produces the first plot then hangs for a while
                                                             # Maple 15: produces fast but ugly
Animate(p, 'frames'=4, 'keep'=true, 'ext'=eps); # Maple 16: produces the first plot, with black background

@Alejandro Jakubi 

well...

donde hay vida hay esperanza

or for a politically incorrect colloquial translation:

 

it's not over till the fat lady sings

Thanks so much Preben,

  1. I had missed a pretty basic thing, that once you use the curly brackets, the order of the arguments can change, with the result (if I understand correctly) that the procedure does not "remember" or "care" about the order of the arguments. So you can write h1('n'=2,'l'=[1,2,3]), but if you write h1(2,[1,2,3]), the procedure does not know which is 'n' and which is 'l'. Fair enough. Thanks for clarifying this Preben, I'm a little new to this.
  2. "If you don't want h1 to accept any other arguments than the ones explicitly given in the definition of h1, then end the argument sequence with a dollar sign." That's great, thanks.
  3. "The default doesn't have to have the declared type!" That's great to know, I hadn't even thought of it.
  4. The "seed" procedure, I'm not sure it's working properly, is it working for you? p2 is the same proc you wrote, slightly simplified:

    > p2:=proc({ [s,seed]::posint:= NULL} )
      if s<>NULL then randomize(s) end if:
      rand()
    end proc:

    > p2(); # does not randomize properly
    > p2(seed=45); # works
                              565524699004
                              873186335134

    Without a restart, I get a number for p2(); the first time I execute and a different number the second time, but from then on the same (second) number keeps coming up. With a restart at the top, I consistently get the same number for p2().

    The proc below works as expected in that it produces different numbers every time I execute p1(); (unless I restart).
    p1:=proc()
      rand();
    end proc:

    p1(); # works, and really you would expect it to!

Thanks so much Preben,

  1. I had missed a pretty basic thing, that once you use the curly brackets, the order of the arguments can change, with the result (if I understand correctly) that the procedure does not "remember" or "care" about the order of the arguments. So you can write h1('n'=2,'l'=[1,2,3]), but if you write h1(2,[1,2,3]), the procedure does not know which is 'n' and which is 'l'. Fair enough. Thanks for clarifying this Preben, I'm a little new to this.
  2. "If you don't want h1 to accept any other arguments than the ones explicitly given in the definition of h1, then end the argument sequence with a dollar sign." That's great, thanks.
  3. "The default doesn't have to have the declared type!" That's great to know, I hadn't even thought of it.
  4. The "seed" procedure, I'm not sure it's working properly, is it working for you? p2 is the same proc you wrote, slightly simplified:

    > p2:=proc({ [s,seed]::posint:= NULL} )
      if s<>NULL then randomize(s) end if:
      rand()
    end proc:

    > p2(); # does not randomize properly
    > p2(seed=45); # works
                              565524699004
                              873186335134

    Without a restart, I get a number for p2(); the first time I execute and a different number the second time, but from then on the same (second) number keeps coming up. With a restart at the top, I consistently get the same number for p2().

    The proc below works as expected in that it produces different numbers every time I execute p1(); (unless I restart).
    p1:=proc()
      rand();
    end proc:

    p1(); # works, and really you would expect it to!

@Axel Vogt 

Indeed that is quite unexpected, why if k is unassigned the procedure not return an error that it is not a posint? The proc seems to be treating g() and g(k), for any unassigned k, as synonyms...?

@Axel Vogt 

Indeed that is quite unexpected, why if k is unassigned the procedure not return an error that it is not a posint? The proc seems to be treating g() and g(k), for any unassigned k, as synonyms...?

Preben, pagan, thanks a lot, ':-n' has worked so far, so that's a great one to know.

On the second, related question above: it's a little disappointing to have to write a line of code for each of the procedure's arguments, because it adds about 5 lines of code to my procedure and because I also have to define 5 new local variables. But I can certainly live with that if it is the only way. Thanks to both of you for your timely help!

I have some follow-up questions, which I have branched off following acer's suggestion.

Preben, pagan, thanks a lot, ':-n' has worked so far, so that's a great one to know.

On the second, related question above: it's a little disappointing to have to write a line of code for each of the procedure's arguments, because it adds about 5 lines of code to my procedure and because I also have to define 5 new local variables. But I can certainly live with that if it is the only way. Thanks to both of you for your timely help!

I have some follow-up questions, which I have branched off following acer's suggestion.

How about setting interface to quiet?
interface(quiet=true):

How about setting interface to quiet?
interface(quiet=true):

I have managed to set up the procedure to condition upon which events triggered. [I couldn't get eventclear to work, but didn't need it in the end]. I set up the events such that the odd-numbered events indicate a failure and even-numbered events indicate a success.

       if abs(op(rhs(op(rhs(sol('eventfired'=[rhs(rhs(sol('eventstop'))[1])])))[1]))) < 1 then
                error; # forcing an exception
       elif type(rhs(rhs(sol('eventstop'))[1]), odd) then
                error; # forcing an exception
       else
                break ;
       end if ;

The first conditional statement figures out the time at which an event triggered and returns an exception if it's too short (less than 1 unit of time). The second conditional statement returns an exception when the event that triggered is odd-numbered (events 1,3, or 5). It is rather fiddly to work out the proper combinations of op, rhs, and [].

EDIT 1: Replaced ln(0); with error; to force an exception. Thanks pagan for this. The Maple programming guide has section on this, which I had missed on a first reading.

I have 40 trajectories to simulate. To begin with, only 26 converged. Now I've got 30. So there's still 10 to go. Not as good as I was hoping given the hard work. But there is hope. I'll report back if I have any success. Thanks very much for your help.

EDIT 2: Got it to work!! Some very silly bugs I had introduced here and there. I have uploaded the worksheet a couple of posts up from this one. Thanks guys!

I have managed to set up the procedure to condition upon which events triggered. [I couldn't get eventclear to work, but didn't need it in the end]. I set up the events such that the odd-numbered events indicate a failure and even-numbered events indicate a success.

       if abs(op(rhs(op(rhs(sol('eventfired'=[rhs(rhs(sol('eventstop'))[1])])))[1]))) < 1 then
                error; # forcing an exception
       elif type(rhs(rhs(sol('eventstop'))[1]), odd) then
                error; # forcing an exception
       else
                break ;
       end if ;

The first conditional statement figures out the time at which an event triggered and returns an exception if it's too short (less than 1 unit of time). The second conditional statement returns an exception when the event that triggered is odd-numbered (events 1,3, or 5). It is rather fiddly to work out the proper combinations of op, rhs, and [].

EDIT 1: Replaced ln(0); with error; to force an exception. Thanks pagan for this. The Maple programming guide has section on this, which I had missed on a first reading.

I have 40 trajectories to simulate. To begin with, only 26 converged. Now I've got 30. So there's still 10 to go. Not as good as I was hoping given the hard work. But there is hope. I'll report back if I have any success. Thanks very much for your help.

EDIT 2: Got it to work!! Some very silly bugs I had introduced here and there. I have uploaded the worksheet a couple of posts up from this one. Thanks guys!

@acer 

absolutely, it ran superbly with Maple 15 and I played around with it for a while, I didn't experience any crashes, and all of the components were very easy to understand. The initialization took merely a few seconds. None of what I tried took longer than about 4 seconds. Tooltips are not a top priority, it's very clear what every thing does. It just about fit on a single page too, which is really nice.

First 14 15 16 17 18 19 20 Last Page 16 of 93