Doug Meade

 

Doug

---------------------------------------------------------------------
Douglas B. Meade <><
Math, USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu
Phone: (803) 777-6183 URL: http://www.math.sc.edu

MaplePrimes Activity


These are answers submitted by Doug Meade

Here is one way to get a matrix returned from your maplet. This uses invisible fields to handle the communication. I tend to do this since I seem to have trouble when I try to pass arguments.

with(Maplets:-Elements):
makeMaplet := proc(m::posint,n::posint)
   Maplet([
     ["Enter the matrix, please"],
     seq([seq(TextBox[entry||i||j](width = 5),j = 1..n)],i = 1..m),
     [TextBox['TBm']('visible'='false', 'width'=1, 'value'=m ),
     TextBox['TBn']('visible'='false', 'width'=1, 'value'=n ),
     Button("Return the matrix",'onclick'='A1'),
     TextBox['TBM']('visible'='false', 'width'=1, 'value'="" ),
     TextBox['TBd']('visible'='false', 'width'=1, 'value'="" )]
     ],
     Action['A1']( Evaluate('function' = "GetMatrix"),
       Shutdown(['TBM']) )
  )
end proc:

GetMatrix := proc()
   local m, n, M;
   uses Maplets:-Tools;
   m := Get( 'TBm'::posint, corrections=true );
   n := Get( 'TBn'::posint, corrections=true );
   M := Matrix( m,n, (i,j)->Get(entry||i||j) );
   Set( 'TBM' = M );
end proc:
AA:=eval( parse( Maplets[Display](makeMaplet(3,5))[] ) );

I would be interested in other solutions that might be a little more elegant than this brute force solution. (The last TextBox is inserted to keep the Button centered.)

Is this helpful?

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

 

If you can give us more details about your problem and what you have already done, there is a better chance someone will be able to suggest a way to generate the plot that you are trying to create.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

Nice questions. Maple has some nice tools for working with PDEs. A good place to learn about these is to look at the online help for the PDEtools package (?PDEtools).

restart;
with( PDEtools ):

To give some specific answers to your questions, I will start by using alias to avoid having to type all arguments to a function:

alias( q=q(x,y,t), h=h(x,y,t) );
                                    q, h

pde := diff(q, t) +(g*DD-u^2)*(diff(h, x))-2*u*(diff(q, x))=0 ;

         / d   \   /        2\ / d   \       / d  \
         |--- q| + \g DD - u / |--- h| - 2 u |--- q| = 0
         \ dt  /               \ dx  /       \ dx /

This simplified the input, but did not help much with the output. To clean up the output, use the declare command (in the PDEtools package):

declare( q(x,y,t), h(x,y,t) );
                    q(x, y, t) will now be displayed as q
                    h(x, y, t) will now be displayed as h

pde := diff(q, t) +(g*DD-u^2)*(diff(h, x))-2*u*(diff(q, x))=0 ;
                   /        2\    
              q  + \g DD - u /  h  - 2 u q = 0
               t                 x        x

At this point I need to explain why I have replaced D with DD. Maple defines D to be the differentiation operator. We will see this appear in the result of substituting another equation into the PDE:

eqn := q=Q(x-a(t)*y,t);
                        q(x, y, t) = Q(x - a(t) y, t)

eval( pde, eqn );
   -D (Q)(x-a(t)y,t) a  y + D (Q)(x-a(t)y,t) 
     1                t      2  

       /        2\    
     + \g DD - u /  h  - 2 u D (Q)(x-a(t)y,t) = 0
                     x        1

I recommend using eval when substituting one expression into another. But, if things get sufficiently complicated, you should also be aware of algsubs.

algsubs( a+b=d, 1+a+b+c );
                                  1 + d + c

algsubs( a*b=c, 2*a*b^2-a*b*d );
                                -d c + 2 b c

algsubs( a*b/c=d, 2*a*b^2/c );
                                    2 d b

There are lots of other tools you could use, but these should be enough to get you started. Let us know if you have more questions.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><                  Math, USC, Columbia, SC 29208
E-mail: mailto:meade@math.sc.edu       Phone:  (803) 777-6183
URL:    http://www.math.sc.edu/~meade/

I agree with Scott's general suggestions, but see no reason to introduce the individual plots at all. Here is the way I would create a single plot:

plot( [ x*log(x^k) $ k=1..3 ], x=0..5 );

One advantage of this is each curve is a different color.

Here is how I would create the animation:

plots:-animate( plot, [x*log(x^k), x=0..5], k=1..3, frames=3 );

Note the use of frames=3 to get only integer values of k in the animation.

If you do want the log base 10, replace log with log[10]. And, lastly, note that in Scott's plots there is nothing plotted for x<0 for even k. This is as expected. The original post indicated that x is small and positive.

I hope this is useful.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

Someone might have a general procedure for distributing multiplication over the terms in a series, but here is a lower-level look at what is needed to make this work.

eq:=y*(sum(a[n]*y^(n-2)*n*(n-1), n = 0 .. infinity)) = (sum(a[n]*y^(n-1)*n*(n-1), n = 0 .. infinity));
      /infinity                       \   infinity                       
      | -----                         |    -----                         
      |  \                            |     \                            
      |   )          (n - 2)          |      )          (n - 1)          
    y |  /     a[n] y        n (n - 1)| =   /     a[n] y        n (n - 1)
      | -----                         |    -----                         
      \ n = 0                         /    n = 0                         

is(eq);    
                                    false

My approach is to separate the sides of the equation and to work on the left side.

L,R := [lhs,rhs](eq)[];
      /infinity                       \  infinity                       
      | -----                         |   -----                         
      |  \                            |    \                            
      |   )          (n - 2)          |     )          (n - 1)          
    y |  /     a[n] y        n (n - 1)|,   /     a[n] y        n (n - 1)
      | -----                         |   -----                         
      \ n = 0                         /   n = 0                         

The next command effectively moves the y inside the sum. This is done by replacing the first multiplicand (y) with 1 (subsop(1=1,L)) and to multiply the addends of the sum by the first multiplicand (subsop([2,1]=op(1,L)*op([2,1],L),L). These can be combined in the single command:

L2 := subsop(1=1,[2,1]=op(1,L)*op([2,1],L),L);
                      infinity                         
                       -----                           
                        \                              
                         )            (n - 2)          
                        /     y a[n] y        n (n - 1)
                       -----                           
                       n = 0                           

eq2 := L2=R;
     infinity                            infinity                       
      -----                               -----                         
       \                                   \                            
        )            (n - 2)                )          (n - 1)          
       /     y a[n] y        n (n - 1) =   /     a[n] y        n (n - 1)
      -----                               -----                         
      n = 0                               n = 0                         

is(eq2);
                                    true

To see the new expression simplified:

simplify( eq2 );
      infinity                          infinity                       
       -----                             -----                         
        \                                 \                            
         )          (n - 1)                )          (n - 1)          
        /     a[n] y        n (n - 1) =   /     a[n] y        n (n - 1)
       -----                             -----                         
       n = 0                             n = 0                         

It would be nice if the SumTools package had commands like the IntegrationTools package for working with a sum as an object.

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

Thanks, Scott, for filling in the details. I just did not have time to do the search for this at the time I wrote my response.

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

This is generally done with the DEplot command in the DEtools package. The only trick required here is to convert the second-order ODE into a system of two first-order ODEs:

with( DEtools ):
ODE := diff( x(t), t$2 ) + x(t) = 0;
                         / d  / d      \\           
                         |--- |--- x(t)|| + x(t) = 0
                         \ dt \ dt     //           
SYS := [ diff( y(t), t ) + x(t) = 0,
         diff( x(t), t )        = y(t) ];
                  [/ d      \              d             ]
                  [|--- y(t)| + x(t) = 0, --- x(t) = y(t)]
                  [\ dt     /              dt            ]
DEplot( SYS, [x(t),y(t)], t=0..1, x=-3..3, y=-3..3,
        scene=[x(t),y(t)] );

 

If you want to include some solution curves, here is how you can modify the DEplot command:

DEplot( SYS, [x(t),y(t)], t=0..2*Pi, x=-3..3, y=-3..3,
        [[x(0)=1,y(0)=0],[x(0)=2,y(0)=0]],
        scene=[x(t),y(t)], linecolor=[gold,blue],
        animatecurves=true );

 

Note that the animatecurves=true creates an animation (which cannot be displayed on MaplePrimes).

I hope this is helpful for you.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

I think you will be happy with the ApproximateInt command in the Student[Multivariate] package. The following example is from the online help for ApproximateInt:

 

I hope this is useful,

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

I do not "see" the elliptical intersection in this plot. Try this modification, with two more options and a different orientation:

y1 := (x,z) -> x^2 + z^2:
y2 := (x,z) -> 1 - z^2:
plot3d(
	[y1(x,z),y2(x,z)],x = -1.2..1.2,z = -1.2..1.2,
	orientation = [30,15],labels = [x,z,y],axes = BOXED,
        style=surface, color=[pink,cyan]
);

See ?plot3d,option for a complete set of options that you could use.

I hope this is helpful,

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

This request is almost the same as the one I answered earlier today about phase portraits (see http://www.mapleprimes.com/forum/phaseportrait#comment-14266)

Here's the setup:

restart;
with( DEtools ):
with( plots ):
sys := D(u)(t) = -a*u(t),
       D(v)(t) = u(t)^2/(1+u(t)^2)-v(t);
                                                   2         
                                               u(t)          
                D(u)(t) = -a u(t), D(v)(t) = --------- - v(t)
                                                     2       
                                             1 + u(t)        
ic := [ [u(0)=-1,v(0)=-1],
        [u(0)=-1,v(0)= 1],
        [u(0)= 1,v(0)=-1],
        [u(0)= 1,v(0)= 1] ];
   [[u(0) = -1, v(0) = -1], [u(0) = -1, v(0) = 1], [u(0) = 1, v(0) = -1], 

     [u(0) = 1, v(0) = 1]]

It was not necessary, but I did make up some initial conditions.

To create the direction field and solution curves for a=0.2 use the command:

DEplot( [eval(sys,a=0.2)], [u(t),v(t)], t=-10..10, u=-2..2, v=-2..2,
        ic, linecolor=gold);

If you do not want to see any solution curves, omit the last two arguments in DEplot. To omit the arrows, add arrows=none. See ?DEplot for a complete description of the acceptable options.

To animate this for different values of a, the animate command can be used as follows:

animate( DEplot,
  [ [sys], [u(t),v(t)], t=-10..10, u=-2..2, v=-2..2,
    ic, linecolor=gold ],
  a=-0.5..0.5, frames=50, trace=10 );

This file is an animate GIF, but MaplePrimes displays only the first frame - bummer!

The optional argument frames= is the number of frames (default is frames=25) and trace= sets the number of snapshots to display as the animation progresses (default is trace=0). See ?animate for a full description of all optional arguments.

Using trace with the changing direction fields creates an interesting effect. Try it!

I hope this has been helpful for you.

Doug

---------------------------------------------------------------------
Douglas B. Meade <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

One of the nicer additions to the plots package a few years ago was the animate command. What's so nice about this is that it allows a user to take a plot command for a static plot and easily create an animation based on a parameter in the problem. Here is how animate could be used in a situation like the one you describe.

First, we enter the ODE and some initial conditions:

restart:
with(DEtools):
with( plots ):
eqn := (D@@3)(y)(x)=-a*y(x)-x^2;
                                                    2
                        @@(D, 3)(y)(x) = -a y(x) - x 

ic := [ [y(0)= 0,D(y)(0)=0,(D@@2)(y)(0)=1],
        [y(0)= 1,D(y)(0)=0,(D@@2)(y)(0)=1],
        [y(0)=-1,D(y)(0)=0,(D@@2)(y)(0)=1] ]:

Here is a single phase portrait (for a=1)

ic := [ [y(0)=0,D(y)(0)=0,(D@@2)(y)(0)=1],
        [y(0)=1,D(y)(0)=0,(D@@2)(y)(0)=1],
        [y(0)=-1,D(y)(0)=0,(D@@2)(y)(0)=1] ];
               [[y(0) = 0, D(y)(0) = 0, @@(D, 2)(y)(0) = 1], 

                 [y(0) = 1, D(y)(0) = 0, @@(D, 2)(y)(0) = 1], 

                 [y(0) = -1, D(y)(0) = 0, @@(D, 2)(y)(0) = 1]]
phaseportrait( eqn1, y(x), x=-1..2.5, ic, title=`Asymptotic solution`,
               linecolor=[gold,yellow,wheat] );

And, here is all that is needed to create an animation in terms of the parameter a:

animate( phaseportrait,
  [ eqn, y(x), x=-1..2.5, ic, title=`Asymptotic solution`,
    linecolor=[gold,yellow,wheat] ],
  a=-1..1, frames=50, trace=5 );

NOTE: This is an animated GIF, but I doubt the animation displays in MaplePrimes.

The frames= option controls the number of frames in the animation (i.e., the smoothness of the animation); the default is 25 frames. The trace= option specifies the number of snapshots to leave during the animation; the default is 0. For more details about the options for animate, see the online help at ?animate.

Lastly, I will comment on the arrows question. The help for phaseportrait clearly states that arrows are possible only for first and second order equations. They do not make sense, in a 2D plot, for higher order equations. I believe the same restriction applies for DEplot.

I hope this is useful,

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

I think the real question boils down to this: What is meant by "The computations in all steps must be done in Maple"? This is a question for you to ask your professor.

There are some things that are easier to do by hand. Some with a calculator, and others with Maple.

Did you use Maple to find the eigenvectors? Is this your question? (I doubt it.)

There is no pre-packaged Maple command that I know for directly answering your questions. The suggestions you have gotten on this thread are about all we can give.

Good luck!

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

Joe's comments are all correct, except that he should have asked: is g just piecewise? and not piecewise continuous.

Joe's sample code does not assume continuity at x=c. Here are some additional suggestions for working with his suggestion:

A := int( f-g, x=0..1 );
A2 := map( collect, convert( A, piecewise, c ), c );

Only one of the parts of this piecewise defined function is relevant here, so let's look only at that one part:

A3 := collect( convert( A, piecewise, c ), c ) assuming c>0, c<1;
              /  1      1   \  2                  1        1   
              |- - m1 + - m2| c  + (-b1 + b2) c + - - b2 - - m2
              \  2      2   /                     4        2   

Now, differentiating yields

dA := diff( A3, c );
                          /  1      1   \            
                        2 |- - m1 + - m2| c - b1 + b2
                          \  2      2   /            
solve( dA=0, {c} );
                              /      -b1 + b2\ 
                             { c = - -------- }
                              \      m2 - m1 / 

This is completely equivalent to the assumption that g is continuous at x=c. Could that be the point of this exercise? Of course, we are still working under the assumption that f≥g for all 0≤x≤1.

Does this help you?

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

Have you looked at the arrow= option in the DEplot command?

From the help page for DEtools,DEplot:

arrows = arrowtype 
  arrowtype consists of 'small', 'smalltwo', 'medium', 'mediumfill', 'large', 'curve', 'comet', 'line', or 'none'
  The arrowtype indicates the type of arrow to be used in a direction field. The default, when a direction field is possible, is 'small'.
  Specifying an arrowtype of 'none' suppresses the calculation and display of any direction fields. 

As I do some testing with these options I do not see that they will give what you want. I tried to use a larger mesh of points (dirgrid=[50,50]) but the routine is smart enough to draw shorter arrows, ensuring that they do not interfere with one another.

So, I would recommend that there be another option added to DEplot that gives the user control over the lengths of the arrows.

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

Stephen,

Why do you need the extra map? Have you tried using:

$med=maple("solve(CurveFitting[Spline]($XY, x, degree=1)=$cf5/2,x)");

But, I believe the real problem is that some of the data in XY, or cf5, is not numeric. I get the solutions may be lost message if I call solve with symbolic values in XY.

Also, I notice that your XY ends with cf4 but you later refer to cf5. If all data is numeric, then I have no trouble.

Good luck

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/
First 31 32 33 34 35 36 37 Last Page 33 of 44