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

Have you looked at Maple's event handler for differential equations? See ?dsolve,numeric,Events

If your equations change at times known in advance, then piecewise might be best (as Markiyan suggests). But, if your transitions occur only after the solution reaches a specific state, the event handler might be the answer you seek. The above-cited help page is pretty extensive, with several examples at the end.

Let us know if you need more help or if these suggestions are not helpful. Posting a worksheet would allow us to see exactly what you are trying to do, and might allow us to give more helpful advise.

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

The evaln approach is a good answer when you want to have the functions defined for later use.

What I don't understand is why you need the Equ definitions. Other than simply displaying the definition of the expression, I don't see how Equ[i] will be used.

If you are going to be doing this a lot, I'd suggest defining a special Maple procedure to make these displays:

showfn := a::uneval -> a = eval(a):
showfn(u[1](t));

1
u[1](t) = - L cos(t w + phi[1])
2
This does not work well inside loops due to the way the uneval feature works:

for i to 4 do showfn(u[i](t)) end do;
print(`output redirected...`); # input placeholder
1
u[i](t) = - L cos(t w + phi[1])
2
             1 
u[i](t) = - L cos(t w + phi[2])
2
1
u[i](t) = - L cos(t w + phi[3])
2
1
u[i](t) = - L cos(t w + phi[4])
2

If you are interested, I'm sure we can come up with a better implementation of showfn.

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

While the elementwise operators are nice, they sometimes are not the most natural idea. Here's another approach in which we construct a list the same size as the list of fractions where each element's value is e:

restart;
b := [ 3, 8/3, 11/4, 19/7, 87/32, 106/39, 193/71]:
c := [ exp(1) $ nops(b) ]:
d := evalf( c-b );
[-0.281718172, 0.051615161, -0.031718172, 0.003996114, -0.000468172, 0.000333110, -0.000028031]

I hope this is helpful to you, and to others,

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

I can't give you an answer directly for your system of DEs because I don't know the values for all of your parameters.

You also don't say what version of Maple you are using. I'm going to show you how I would do this, using the "parameters=" feature in dsolve,numeric. If you are using a version that does not support this, let me know and I'll show you the old way of doing this as well.

To start, here's the simple ODE I'm using for my example:

restart;
with( plots ):
ode := diff(y(x),x)=y(x)/x^p;
ic := y(1)=2;

To have Maple produce a numerical solution, with p as a paramter:

sol := dsolve( {ode,ic}, numeric, parameters=[p] );

To get a plot for a specific value of p:

sol(parameters=[p=2]);
odeplot( sol, 0..2 );

To create a list of plots for different values of p:

P := NULL:
for pp from 1 to 10 do
sol( parameters=[p=pp] ):
P := P, odeplot( sol, 1..2, title=sprintf("Solution for p=%d",pp) );
end do:

To see a specific plot:

P[2];

Or, to see the list as an animation:

display( P, insequence );

I hope this helps you get closer to what you are trying to accomplish.

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

The previous posts are correct, but I did not see that they explained WHY you can't use gamma as a variable name in Maple.

gamma is Euler's constant

 

If you really need to use gamma in your work, and aren't using a version that support "local gamma" I would encourage you to look at the unprotect command to break the built-in link between the name gamma and Euler's constant. But, before you do use unprotect, read the help page for unprotect (?unprotect) and understand the example that redefines Pi and then computes arctan(1).

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

Why do you say these are the same expression?

exp(I*2*t)-exp(-I*2*t) = 2*I*sin(2*t)

and cos(t)-I*sin(t) has a non-zero real part (and a different imaginary part).

Maple's plot command does not know what to do with a complex-valued function. The result is usually empty, but the viewing window could be different - because the functions are different.

As you show at the end of your post, the complexplot command is able to create a plot of a complex-valued expression, The example you included in your post plots the (parametric) curve in the complex plane.

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

Maple interprets 1..130 as an interval, as in specifying the viewing window for a plot or the interval of integration for a definite integral. I'm going to guess that what you are trying to get is a list of points. For this you need to use the sequence operator ($).

> $(1 .. 10);

1, 2, 3, 4, 5, 6, 7, 8, 9, 10

If you really want 131 points, Maple will not show the full matrix - just its size and basic information about the data. To get Maple to show a larger matrix you need to increase the rtablesize setting:

> interface( rtablesize=150 ):

The default setting for rtablesize is 10.

Here are a few examples if you want to use points with steps other than 1:

> ($0..10)/10;  # rational

> ($0..10)/10.; # floating-point

> 3/2+~($0..10) # half-integers, starting with 3/2

> a+~b*~($0..10); # using the elementwise + and * operators

> a+k*b $ k=0..10; # same as the previous, using a different form of the sequence operator

Please let us know if this is not what you are trying to do, or have more questions. (Just don't ask us to do your HW 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

In addition to the previous comments about pi and Pi and about possible problems with implied multiplication, I wonder why you chose to convert everything to floating-point numbers? If you keep the coefficients are integer and rational numbers, there's a lot better chance that Maple will be able to simplify these expressions. After you've gotten the exact equations, if you then need Maple to find an approximate solution you can use fsolve instead of solve.

Looking forward to seeing if all of these comments are of any help to you. (If you post the actual .mw file it will be easier for us to make more specific and useful comments.)

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

The specific question asked by the OP is about finding the difference between elements in a list. Here is how I would do this for a previously defined list L:

dL := L[2..-1] - L[1..-2];

The first term on the RHS, L[2..-1], is the list with the first element omitted; the second term on the RHS, L[1..-2], is the list with the last term omitted. Note that the use of negative indices allows this to work without regard for the number of elements in L. (If L has only one element, this returns the empty list. If L is the empty list, it returns an error about an invalid subscript selector. Both are reasonable responses, IMHO.)

As has been previously mentioned, there are lots of problems with the code provided by the OP. I don't see that PosEigenvs contains only the positive eigenvectors and I don't understand the need for seqPosEigenvs (where each entry in this list is itself a list).

Since you are finding the eigenvalues of a symmetric matrix, you know all eigenvalues are real-valued. I would force Maple to remove the "+0.*I", and then extract the positive eigenvalues as follows:

PosEigenvs := select[flatten]( q -> (q>0), evalc(Re(Eigenvs)) );

The evalc(Re(Eigenvs)) ensures that the eigenvalues are actually real-valued (so the subsequent comparison won't have any problems). The first argument to select is the boolean function that returns true when the input is positive. Without the [flatten] option the result will be a list the same size as the original list, with removed entries replaced by NULL; adding [flatten] tells Maple to remove these NULL elements.

I hope these give you some ideas about how to fix your code to get what you need. Please feel free to post more questions. We are here to help.

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

From the online help for "next":

for n from 1 to 3 do
  if n=2 then next; end if;
  print(n);
end do;
   1
   3

This is, of course, equivalent to the first response. I prefer to use next because it keeps the code a little cleaner and easier to maintain - I don't have to worry about some commands falling on the wrong side of the then, else, or end if.

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

It is not clear to me why you want to "convert print(h, AA)" - or what this means.

The way I would produce the two lists you describe is as follows:

h := [$0.1 .. 50.0]/10;
AA := map( x->(1+x)/x^2, h );

The first command uses the $ operator to create the list of h values. Note that the stepsize has to be 1, so I multiplied your initial and final values by 10 and then divide the entire list by 10.

The second command applies your function to each element in the list h.

If you should want to combine these lists into a list of ordered pairs, one efficient way to do this is:

hAA := zip( (h,a)->[h,a], h, AA );

To convert this into a 50x2 matrix:

convert( hAA, Matrix );

I hope some of this is helpful to you.

Doug

 

There is no need to get overly fancy with the syntax. I would personally use the implicitplot as suggested above. But, if I did want to graph an explicit function (or two), I would suggest the following:

c := x^2+(y-3)^2=25:
cplot := solve( c, y ):
plot( [cplot], x=-5..5, scaling=constrained );

The result is the same as you would get with Carl's suggestion, but I find it easier to understand what is going on with the square brackets when they are used as above. An alternative is to put the square brackets in the definition of cplot, not in plot - or to even avoid the use of cplot altogether:

plot( [solve( c, y )], x=-5..5, scaling=constrained );

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

It sounds as though you are asking us to do your homework. I am not about to comply with your request. I will suggest that you spend some time looking for information about Maple's numerical solution to differential equations. There is lots of information about this. I believe you'll quickly find something that gives you what you need - once you update the command for your specific problem.

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

The choose command in the combinat package will help you with this:

restart;
with(combinat):

S := [{2, 8}, {2, 4, 6, 7, 8}, {4, 6, 8}, {1, 4, 8}, {1, 4, 8, 9}];
[{2, 8}, {2, 4, 6, 7, 8}, {4, 6, 8}, {1, 4, 8}, {1, 4, 8, 9}]

choose( S, 3 );
[[{2, 8}, {1, 4, 8}, {4, 6, 8}],

[{2, 8}, {1, 4, 8}, {1, 4, 8, 9}],

[{2, 8}, {1, 4, 8}, {2, 4, 6, 7, 8}],

[{2, 8}, {4, 6, 8}, {1, 4, 8, 9}],

[{2, 8}, {4, 6, 8}, {2, 4, 6, 7, 8}],

[{2, 8}, {1, 4, 8, 9}, {2, 4, 6, 7, 8}],

[{1, 4, 8}, {4, 6, 8}, {1, 4, 8, 9}],

[{1, 4, 8}, {4, 6, 8}, {2, 4, 6, 7, 8}],

[{1, 4, 8}, {1, 4, 8, 9}, {2, 4, 6, 7, 8}],

[{4, 6, 8}, {1, 4, 8, 9}, {2, 4, 6, 7, 8}]]
nops( % );
10

choose( S );
[[], [{2, 8}], [{1, 4, 8}], [{4, 6, 8}], [{1, 4, 8, 9}],

[{2, 4, 6, 7, 8}], [{2, 8}, {1, 4, 8}], [{2, 8}, {4, 6, 8}],

[{2, 8}, {1, 4, 8, 9}], [{2, 8}, {2, 4, 6, 7, 8}],

[{1, 4, 8}, {4, 6, 8}], [{1, 4, 8}, {1, 4, 8, 9}],

[{1, 4, 8}, {2, 4, 6, 7, 8}], [{4, 6, 8}, {1, 4, 8, 9}],

[{4, 6, 8}, {2, 4, 6, 7, 8}], [{1, 4, 8, 9}, {2, 4, 6, 7, 8}],

[{2, 8}, {1, 4, 8}, {4, 6, 8}],

[{2, 8}, {1, 4, 8}, {1, 4, 8, 9}],

[{2, 8}, {1, 4, 8}, {2, 4, 6, 7, 8}],

[{2, 8}, {4, 6, 8}, {1, 4, 8, 9}],

[{2, 8}, {4, 6, 8}, {2, 4, 6, 7, 8}],

[{2, 8}, {1, 4, 8, 9}, {2, 4, 6, 7, 8}],

[{1, 4, 8}, {4, 6, 8}, {1, 4, 8, 9}],

[{1, 4, 8}, {4, 6, 8}, {2, 4, 6, 7, 8}],

[{1, 4, 8}, {1, 4, 8, 9}, {2, 4, 6, 7, 8}],

[{4, 6, 8}, {1, 4, 8, 9}, {2, 4, 6, 7, 8}],

[{2, 8}, {1, 4, 8}, {4, 6, 8}, {1, 4, 8, 9}],

[{2, 8}, {1, 4, 8}, {4, 6, 8}, {2, 4, 6, 7, 8}],

[{2, 8}, {1, 4, 8}, {1, 4, 8, 9}, {2, 4, 6, 7, 8}],

[{2, 8}, {4, 6, 8}, {1, 4, 8, 9}, {2, 4, 6, 7, 8}],

[{1, 4, 8}, {4, 6, 8}, {1, 4, 8, 9}, {2, 4, 6, 7, 8}],

[{2, 8}, {1, 4, 8}, {4, 6, 8}, {1, 4, 8, 9}, {2, 4, 6, 7, 8}]]
nops(%);
32

 

This finds the 10 three element subsets of this set of 5 sets. It also shows that there are a total of 32 different subsets with between 0 and 5 elements.

I hope this is still of some use to 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

Satya,

 

This was fun. The key to my approach is the coeffs command. Not only does it return all coefficients of a multinomial polynomial but it can also return the term for each coefficient. I zip the two lists together and then select all pairs whose first term (coefficient) matches the given value. My final implementation is presented in the following procedure - with two test cases.

restart;

Terms := proc( f, vars, coef )
local c, t, T;
c := coeffs( f, vars, 't' );
T := zip( (x,y)->[x,y], [c], [t] );
return map2( op, 2, select( C->evalb( C[1]=coef ), T ) );
end proc:

f := 3*x^3 + 5*x^4 + 5*x^44;
3 4 44
3 x + 5 x + 5 x
Terms( f, x, 5 );
[ 4 44]
[x , x ]
Terms( f, x, 4 );
[]
Terms( f, x, 3 );
[ 3]
[x ]

f2 := 3*x^3*y + 4*x^2 + 5*y^33;
3 2 33
3 x y + 4 x + 5 y
Terms( f2, [x,y], 5 );
[ 33]
[y ]
Terms( f2, [x,y], 4 );
[ 2]
[x ]
Terms( f2, [x,y], 3 );
[ 3 ]
[x y]

Does this do what you request?

I hope this is helpful. I encourage you to consult the Maple help documents for coeffs, zip, select, and map2 (in particular) if you need help understanding what each command does.

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
3 4 5 6 7 8 9 Last Page 5 of 44