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

This is a problem that I believe Maple can handle with relative ease. I have uploaded my approach as a worksheet: View 178_TotalDiff.mw on MapleNet or Download 178_TotalDiff.mw Since some of this requires some explanation, I will also show the code here: Start by telling Maple that the p's are functions of e1, e2, and e3. This cuts down on input, and output.
restart;
PDEtools:-declare( seq( p(e1,e2,e3), m=1..3 ) );
                  p(e1, e2, e3) will now be displayed as p
alias( seq(p[m]=p[m](e1,e2,e3),m=1..3 ) );
Now, for the T operator. This is pretty straightforward. I added the collect on diff to help organize the output. (If you want the output in a different form, you might need to replace the -> form with a formal proc.)
T := f -> collect( add( i*e||i*diff(f,e||i), i=1..3 ), diff );
Here are two tests of this operator (the output is much nicer in 2D in the worksheet):
T(p[1]);
             / d                   \        / d                   \
          e1 |---- p[1](e1, e2, e3)| + 2 e2 |---- p[1](e1, e2, e3)|
             \ de1                 /        \ de2                 /

                    / d                   \
             + 3 e3 |---- p[1](e1, e2, e3)|
                    \ de3                 /
T(T(p[1]));
      / d                   \        / d                   \
   e1 |---- p[1](e1, e2, e3)| + 4 e2 |---- p[1](e1, e2, e3)|
      \ de1                 /        \ de2                 /

             / d                   \     2 / d   / d                   \\
      + 9 e3 |---- p[1](e1, e2, e3)| + e1  |---- |---- p[1](e1, e2, e3)||
             \ de3                 /       \ de1 \ de1                 //

                / d   / d                   \\
      + 4 e1 e2 |---- |---- p[1](e1, e2, e3)||
                \ de2 \ de1                 //

                / d   / d                   \\
      + 6 e1 e3 |---- |---- p[1](e1, e2, e3)||
                \ de3 \ de1                 //

            2 / d   / d                   \\
      + 4 e2  |---- |---- p[1](e1, e2, e3)||
              \ de2 \ de2                 //

                 / d   / d                   \\
      + 12 e2 e3 |---- |---- p[1](e1, e2, e3)||
                 \ de3 \ de2                 //

            2 / d   / d                   \\
      + 9 e3  |---- |---- p[1](e1, e2, e3)||
              \ de3 \ de3                 //
We turn now to the total derivative operator. I will not use D ad Maple already defines this as a builtin differentiation operator. The only real trick here is to explicitly replace each ei with ei(x) before taking the derivatives. Once again, collecting on diff helps to put the output in a form that highlights the linearity in the derivatives.
Dtot := f -> collect( diff(eval(f,[seq(e||i=e||i(x),i=1..3)]),x), diff );
To suppress the dependence on x in the output, we use:
alias( seq(e||i=e||i(x),i=1..3) );
This generates 9 warning messages, 1 for each of the three variables, ei, for the three functions p[m]. And, now, some tests, with the disclaimer that the output is much more pleasing in the worksheet:
Dtot(p[1]);
                             / d    \                          / d    \
      D[1](p[1])(e1, e2, e3) |--- e1| + D[2](p[1])(e1, e2, e3) |--- e2|
                             \ dx   /                          \ dx   /

                                  / d    \
         + D[3](p[1])(e1, e2, e3) |--- e3|
                                  \ dx   /

Dtot(T(p[1]));
(D[1](p[1])(e1, e2, e3) + 2 e2 D[1, 2](p[1])(e1, e2, e3)

                                                                    / d    \   
   + e1 D[1, 1](p[1])(e1, e2, e3) + 3 e3 D[1, 3](p[1])(e1, e2, e3)) |--- e1| + 
                                                                    \ dx   /   

  (2 e2 D[2, 2](p[1])(e1, e2, e3) + e1 D[1, 2](p[1])(e1, e2, e3)

                                                                / d    \      
   + 2 D[2](p[1])(e1, e2, e3) + 3 e3 D[2, 3](p[1])(e1, e2, e3)) |--- e2| + (2 
                                                                \ dx   /      

  e2 D[2, 3](p[1])(e1, e2, e3) + e1 D[1, 3](p[1])(e1, e2, e3)

                                                                / d    \
   + 3 D[3](p[1])(e1, e2, e3) + 3 e3 D[3, 3](p[1])(e1, e2, e3)) |--- e3|
                                                                \ dx   /

Did I come close to the functionality you requested? 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/
Robert's response is correct, but if you are a newcomer to Maple, you may not know how to enter "your equations". Here is an excerpt from a Maple 11 worksheet showing a few ways to approach your problem: The first step is to define the three equations (note the use of * for multiplicaiton and := for assignment):
eq1 := 3*x+2*y+z=2;
                              3 x + 2 y + z = 2
eq2 := 4*x+2*y+2*z=8;
                             4 x + 2 y + 2 z = 8
eq3 := x-y+z=4;
                                x - y + z = 4
Here is what Robert's solution would look like:
solve( {eq1,eq2,eq3} );
                           {y = 2, x = -4, z = 10}
It is exactly the same as you would see if you explicitly entered the unknowns as a set:
solve( {eq1,eq2,eq3}, {x,y,z} );
                           {y = 2, x = -4, z = 10}
If you specify the unknowns as a list (which has order), the output is slightly different:
solve( {eq1,eq2,eq3}, [x,y,z] );
                          [[x = -4, y = 2, z = 10]]
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/
Here is a complete solution to your question. The main approach is the same as suggested by Robert. I have guessed at the function that you could have intended. You can, of course, change this to what you really need. Note that I use the numpoints= option to ensure the final frames have sufficient resolution. Also, the frames are numbered. I like to use sprintf. The title, and the background, are assembled for each frame as it is created inside the loop. The final display simply displays these frames.
restart;
with(plots):
f:=piecewise(-Pi<x and x<0, -1, x>0 and x<Pi, 1 );
functionplot:=plot(f,x=-Pi..Pi):
N := 50;
a[0]:=int(f,x=-Pi..Pi)/2/Pi:
F[0]:=evalf(a[0]):
for i from 1 to N do
  a[i]:=evalf(int(f*cos(i*x),x=-Pi..Pi)/Pi):
  b[i]:=evalf(int(f*sin(i*x),x=-Pi..Pi)/Pi):
  F[i]:=F[i-1]+a[i]*cos(i*x)+b[i]*sin(i*x):
  figure[i]:=display( [plot(F[i],x=-Pi..Pi,color=blue,numpoints=500),functionplot], title=sprintf("Fourier Approximations to a Step Function\nFrame %a of %a", i, N) ):
end do:
display([seq(figure[i],i=1..N)],insequence=true);
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/
I will let Scott speak for himself, but the misuse of quotes is one of the standard things I check for when debugging Maple code. It really is a fairly common mistake for students to want to use the double-quote key as a shortcut for two single-quotes. An opinion from the trenches, 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/
Further to this discussion, am I correct in thinking that the integral of any continuous function between 0- and 0+ (or more generally a- to a+) is always equal to zero? You should be careful to indicate exactly where the function is continuous: at 0 or at a for the two statements in your post. The result is still zero even for some discontinuous functions. If the integrand has a simple jump discontinuity at a (different left and right limits), then there is still a continuous antiderivative and the definite integral from a- to a+ will be zero. Things get interesting when the integrand is something like a impulse function, such as a Dirac distribution. The Dirac distribution at x=a is the "function", f, with the property that f(x)=0 for all x<>a and int( f(x), x=-infinity..infinity ) = 1. Thus, int( f(x), x=-infinity..a- ) = 0, int( f(x), x=a+..infinity ) = 0, and int( f(x), x=a-..a+ ) = 1. For more information, take a look at the Wikipedia page on Dirac delta functions. At some point in this discussion you will have to take some time to learn about the different types of integral. Most people start with the Riemann integral, but that is only the tip of the proverbial mathematical iceberg. (Hey, it's cold here today.) For starters, read the Wikipedia page on Lebesgue integration.
---------------------------------------------------------------------
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/
The previous posts have all written Maple procedures that return your desired list of matrices. You need to i) copy the code for the procedure into an active Maple worksheet, ii) execute the code for each procedure, and iii) call the procedure, as shown below (output omitted).
mzp(3);
m2p(3);
myp(3);
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/
The code for isolve is open. Anyone can see the code using, e.g., showstat as shown below. From my experience, this is usually more than enough information for an interested user to create their own extension or improvement. When complete, submit your work to the Application Center for other users to download. Also, let Maplesoft Customer Support know about your improvements. If they see benefits in your work, you could find your improvements in a future release of Maple. I know this works, sometimes, from personal experience. My experience also suggests that the developers will likely make additional revisions beyond the ones you implemented. So, let's start seeing some improved isolve commands.
showstat( isolve );

isolve := proc(eqs, solsin)
local eqn, eqns, sol, solsinset, vars;
   1   if nargs < 1 or 2 < nargs or not type(eqs,{'relation', 'algebraic', 'set'}) or nargs = 2 and not type(solsin,{'set', 'name'}) then
   2     error "invalid arguments"
       end if;
   3   if type(eqs,set) then
   4     vars := indets(eqs);
   5     eqns := eqs
       else
   6     vars := indets(eqs);
   7     eqns := {eqs}
       end if;
   8   vars := map(proc (x) if type(x,name) then x end if end proc,vars);
   9   eqns := map(proc (x) if type(x,equation) then op(1,x)-op(2,x) else x end if end proc,eqns);
  10   if nargs = 2 then
  11     if type(solsin,set) then
  12       solsinset := solsin
         else
  13       solsinset := {solsin}
         end if
       else
  14     solsinset := {_Z1}
       end if;
  15   if type(eqns,('set')(('polynom')('anything',vars))) then
  16     eqns := remove(x -> evala(x = 0),eqns);
  17     if select(x -> type(x,'realcons'),eqns) <> {} then
  18       return {}
         end if;
  19     if `minus`(indets(eqns,'realcons'),indets(eqns,'rational')) <> {} then
  20       eqns := map(`@`(evala,Factor),eqns);
  21       eqns := map(proc (x) if type(x,`*`) then remove(y -> `minus`(indets(y,realcons),indets(y,rational)) <> {},x) else x end if end proc,eqns);
  22       if `minus`(indets(eqns,'realcons'),indets(eqns,'rational')) <> {} then
  23         return {}
           end if
         end if
       end if;
  24   sol := `isolve/isolve`(eqns,vars,solsinset);
  25   if sol = {} then
  26     userinfo(1,isolve,`Warning: unable to find solution over the integers`);
  27     return NULL
       end if;
  28   return op(sol)
end proc
showstat( `isolve/isolve` );
showstat( `isolve/quadratic` );
showstat( `isolve/reduce' );
showstat( `isolve/compose` );
...
---------------------------------------------------------------------
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/
Since 0+ and 0- are not actual numbers, you will need to go back to what this notation actually represents: a limit.
assume( s>0 );
f := t;
                                      t
F := int( f*exp(-s*t), t=a..infinity );
                             (1 + s a) exp(-s a)
                             -------------------
                                      2         
                                     s          
limit( F, a=0, right );
                                     1 
                                     --
                                      2
                                     s 
This gives the limit from the right (0+). Replace the last argument of the limit command with left to take the limit from the left (0-). Of course, in this case the same result can be obtained without limits. 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/
Here is one way to do what you are trying to do. The key to my approach is to delay evaluation of the integral until you have a specific value for n, or you really want the general case.
B := 4/Pi*Int(cos(3*x)*cos(n*x),x=0..Pi/2);
                       /  /1/2 Pi                     \
                       | |                            |
                     4 | |        cos(3 x) cos(n x) dx|
                       | |                            |
                       \/0                            /
                     ----------------------------------
                                     Pi                
This defines the Fourier coefficient as an integral but does not attempt to evaluate the integral. To get Maple to evaluate the integral, use value:
BB := value( B );
                                     /1     \
                               12 cos|- Pi n|
                                     \2     /
                     BB :=     --------------
                                   /      2\ 
                                Pi \-9 + n / 
But you cannot evaluate this result for n=3. You can take a limit to get the B3 coefficient.
eval( BB, n=3 );
Error, numeric exception: division by zero
limit( BB, n=3 );
                                      1
Another way to get B3 without a limit is to evaluate the unevaluated integral when n=3 and ask Maple to evaluate this result.
value( eval(B,n=3) );
                                      1
---------------------------------------------------------------------
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/
As you have described the problem, it sounds as though you have a situation where you have a set of initial conditions that allow you to solve the system from the initial time up to the time of the first "perturbation". When the first perturbation occurs a new system has to be solved. This system uses the solution from the previous time interval as its initial condition - with one (or some) values reset. Now, solve up to the second time when a perturbation occurs. Repeat the above for as many perturbations as necessary. If you want to put everything together in one "solution", I would recommend using the piecewise command. Here is how I put this together for the example you describe in your post:
restart;
eq_1:=diff(X(t), t)=-X(t)*Y(t);
eq_2:=diff(Y(t),t)=X(t)*Y(t);

IC1:={eq_1, eq_2,X(0)=10,Y(0)=1};
F1:=dsolve(IC1, {X(t),Y(t)},type=numeric, method=rosenbrock);
h1 := theta->eval([X(t),Y(t)],F1(theta));

IC2 := {eq_1,eq_2, X(10)=5, Y(10)=h1(10)[2]};
F2:=dsolve(IC2, {X(t),Y(t)},type=numeric, method=rosenbrock);
h2 := theta->eval([X(t),Y(t)],F2(theta));

H := proc(theta)
  if type(theta,name) then
    'procname'(args)
  else
    piecewise( theta<10, h1(theta), h2(theta) )
  end if;
end proc;

plot( [H(t)[1],H(t)[2]], t=0..15 );
I do not know what you want to do with the solutions. It's possible my implementation should be modified to fit your specific needs. But, I think this shows one way to approach this type of problem. All of this could be put in a proc is desired. 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/
There are many ways to do this. I'll give one of the most straightforward (but not necessariy fastest or most efficient).
F := (n::posint,x) ->
  if n=1 then 1
  elif n=2 then x
  else simplify( x*F(n-1,x)+F(n-2,x) )
  end if:
F(1,x);
                                      1
F(3,x);
                                    2    
                                   x  + 1
F(10,x);
                        9      7       5       3      
                       x  + 8 x  + 21 x  + 20 x  + 5 x
F(0,x);
Error, invalid input: F expects its 1st argument, n, to be of type posint, but received 0

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/
Look carefully at your post. You have
SetProperty('PlotInput2','value',plot(GetProperty(MathContainer0','value')),x=0..10);
and
SetProperty('PlotInput2','value',plot(parse(GetProperty('TextArea0','value')),x=0..10));
It appears as though the first is missing a parse and the x=1..10 is an argument to SetProperty, not to plot. The real problem is likely to be that the value field in a MathContainer is expressed in MathML. Instead of parse, you should probably be using MathML[Import]. (See ?MathExpressionComponent.) I don't know why the training materials seem to encourage putting both the get and the set in one line. It make much more sense, to me to implement this as:
F := MathML[Import]( GetProperty( 'MathContainer0', 'value' ) );
P := plot( F, x=1..10 ):
SetProperty( 'PlotInput2', 'value', P );
If there is a chance F might be an expression sequence, then you should probably put [] or {} around F in the plot command. Can you provide the exact worksheet that you are using? This would help us to better understand your situation. 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/
It sounds as though you can get started by implementing your definition of a binary operation (addition) and equality:
MyAdd   := (a,b)->a+b;
MyEqual := (a,b)->evalb(sort(a)=sort(b));
Note that sort is used to test if two elements have the same entries regardless of order. Now, to create the equivalence classes among pairs of elements of a set, A, when combined with a binary operation BinOper and testing equality with Equal can be implemented as:
EquivalenceClass := proc( A, BinOper, Equal )
  local b, B, s, S;
  S := combinat[choose](A,2);
  B := {};
  while S<>{} do
    for b in B do
      S := select( s->not Equal(BinOper(s[1],s[2]),BinOper(b[1],b[2])), S );
    end do;
    if S<>{} then B := B union {S[1]} end if;
  end do;
  return B;
end proc;
Here is how this would look in your specific case:
a := [ [1,0,1,0], [1,1,0,0], [0,0,1,1] ]:
A := { a[1], a[2], a[3] };
                 {[1, 1, 0, 0], [1, 0, 1, 0], [0, 0, 1, 1]}
B := EquivalenceClass( A, MyAdd, MyEqual );
        {{[1, 1, 0, 0], [1, 0, 1, 0]}, {[1, 1, 0, 0], [0, 0, 1, 1]}}
Hoping 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/
Thanks for the clarification. This should make it easier for people to know how to respond to your posts. As a professor who has TAs help with the Maple Lab part of the course I would want to know about the TAs inability to provide adequate instructional support. Does the professor know about these difficulties? As a student (i.e., paying customer), the university has a responsibility to provide support for its courses (i.e., its product). I tell this to my students and, as the Directory of Undergraduate Studies, receive complaints from students about this type of problem. These are not issues for you to be concerned about tonight. You have your assignment to complete, and the exams tomorrow. 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/
I apologize in advance if this is out of line, but I can't help but to wonder if the advise being given here is appropriate. I am teaching a Calculus I course, which happens to also be numbered MATH 141, and am giving my final exam tomorrow. This is not one of my students sice we do not have a Maple component to the final. But, these questions appear to be directly related to the final assessment for this course. The questions are not about how to use Maple to answer these questions; they are asking about the mathematics: continuity, differentiability, .... To the student:
  1. Have you asked your instructor, or other local support, for help with these questions?
  2. Do your instructions allow for outside help, or are you supposed to be doing these on your own?
  3. Good luck on your exams tomorrow!
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 33 34 35 36 37 38 39 Last Page 35 of 44