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

The difference between solve(a,x); and solve (a,x); in 2D Math Notation input is that the space is interpreted as multiplication, i.e., solve*(a,x);

To see this, insert the following command immediately after the offending solve command (the one with the space):

op(0,%);

I predict Maple will respond with `*` -- multiplication.

As others have already written, while 2D Math Notation might look nicer than traditional Maple Notation, it is very easy to get into situations like this one where you really have to know something about how Maple works to diagnose the problems that you encounter. (You can omit the * in the definition of a; you don't even need a space: a := b = 3x+5; is fine, in 2D Math Input.)

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.ed

I'm not sure how useful this is going to be, but it is instructive.

While Alex's use of dsolve and odeplot is very clean and simple, it does not yield a lot of insight into the problem. I thought  I would look at a series solution to the problem. What I found is interesting (at least to me).

First, this system can be reduced to a single ODE for c. First, integrate the first equation and substitute the expression for T back into the second ODE:

eval(dsolve( {ode1,ic1}, T(t) ),ic2);
                                         217     11003
                      T(t) = -125 c(t) - --- t + -----
                                         500      20 
eval( ode2, % );
       d                      2    /            2660              133\
      --- c(t) = -0.01725 c(t)  exp|- ------------------------- + ---|
       dt                          |              217     11003   15 |
                                   |  -125 c(t) - --- t + -----      |
                                   \              500      20        /
 

This first-order (nonlinear) equation for c can now be analyzed in a variety of ways. But, since we are ultimately interested in the maximum value of T, I will continue to work with the full system. The series solution to the system can be found as

dsolve( {ode1,ode2,ic1,ic2}, {T(t),c(t)}, series ):
evalf(%);
   /                                                   2                  3
  { T(t) = 300.1500000 + 8.229303194 t + 0.7522873015 t  + 0.06760870893 t
   \                                                                      

                       4                    5    / 6\          
     + 0.005422448484 t  + 0.0003310839444 t  + O\t /, c(t) = 2.

                                         2                    3
     - 0.06930642555 t - 0.006018298412 t  - 0.0005408696715 t

                         4                      5    / 6\\
     - 0.00004337958788 t  - 0.000002648671555 t  + O\t / }
                                                         /
 

The T solution can now be plotted (after it is converted to a polynomial):

plot( convert(eval(T(t),%),polynom), t=0..20, view=[DEFAULT,0..1000] );

No maximum! Let's increase the order of the series:

Order:=10:
dsolve( {ode1,ode2,ic1,ic2}, {T(t),c(t)}, series ):
evalf(%);
 /                                                   2                  3
{ T(t) = 300.1500000 + 8.229303194 t + 0.7522873015 t  + 0.06760870893 t
 \                                                                      

                     4                    5                      6
   + 0.005422448484 t  + 0.0003310839444 t  + 0.000004594887699 t

                        7                 -7  8                 -8  9    / 10\ 
   - 0.000002868409417 t  - 6.175501179 10   t  - 8.765791486 10   t  + O\t  /,

                                                2                    3
  c(t) = 2. - 0.06930642555 t - 0.006018298412 t  - 0.0005408696715 t

                       4                      5                 -8  6
   - 0.00004337958788 t  - 0.000002648671555 t  - 3.675910158 10   t

                   -8  7                 -9  8                 -10  9    / 10\
   + 2.294727534 10   t  + 4.940400944 10   t  + 7.012633193 10    t  + O\t  /

  \
   }
  /
plot( convert(eval(T(t),%),polynom), t=0..20, view=[DEFAULT,0..1000] );

Now we have a maximum, but it appears to be occuring quite a bit earlier than for the solution returned by dsolve,numeric. But, let's try increasing the order again.

Order:=15:
dsolve( {ode1,ode2,ic1,ic2}, {T(t),c(t)}, series ):
evalf(%);
  /                                                   2                  3
 { T(t) = 300.1500000 + 8.229303194 t + 0.7522873015 t  + 0.06760870893 t
  \                                                                      

                      4                    5                      6
    + 0.005422448484 t  + 0.0003310839444 t  + 0.000004594887699 t

                         7                 -7  8                 -8  9
    - 0.000002868409417 t  - 6.175501179 10   t  - 8.765791486 10   t

                    -9  10                 -10  11                 -11  12
    - 9.841490669 10   t   - 8.667194433 10    t   - 4.788979252 10    t 

                    -12  13                 -12  14    / 15\          
    + 1.726240218 10    t   + 1.008828893 10    t   + O\t  /, c(t) = 2.

                                        2                    3
    - 0.06930642555 t - 0.006018298412 t  - 0.0005408696715 t

                        4                      5                 -8  6
    - 0.00004337958788 t  - 0.000002648671555 t  - 3.675910158 10   t

                    -8  7                 -9  8                 -10  9
    + 2.294727534 10   t  + 4.940400944 10   t  + 7.012633193 10    t

                    -11  10                 -12  11                 -13  12
    + 7.873192532 10    t   + 6.933755547 10    t   + 3.831183402 10    t 

                    -14  13                 -15  14    / 15\\
    - 1.380992175 10    t   - 8.070631149 10    t   + O\t  / }
                                                            /
plot( convert(eval(T(t),%),polynom), t=0..20, view=[DEFAULT,0..1000] );

This is starting to show some of the all-too-common oscillations that can arise in a series expansion. The situation becomes even more interesting when one more term is added to the series solution:

Order:=16:
dsolve( {ode1,ode2,ic1,ic2}, {T(t),c(t)}, series ):
evalf(%);
 /                                                   2                  3
{ T(t) = 300.1500000 + 8.229303194 t + 0.7522873015 t  + 0.06760870893 t
 \                                                                      

                     4                    5                      6
   + 0.005422448484 t  + 0.0003310839444 t  + 0.000004594887699 t

                        7                 -7  8                 -8  9
   - 0.000002868409417 t  - 6.175501179 10   t  - 8.765791486 10   t

                   -9  10                 -10  11                 -11  12
   - 9.841490669 10   t   - 8.667194433 10    t   - 4.788979252 10    t 

                   -12  13                 -12  14                 -13  15  
   + 1.726240218 10    t   + 1.008828893 10    t   + 1.869215726 10    t   +

   / 16\                                                2                    3
  O\t  /, c(t) = 2. - 0.06930642555 t - 0.006018298412 t  - 0.0005408696715 t

                       4                      5                 -8  6
   - 0.00004337958788 t  - 0.000002648671555 t  - 3.675910158 10   t

                   -8  7                 -9  8                 -10  9
   + 2.294727534 10   t  + 4.940400944 10   t  + 7.012633193 10    t

                   -11  10                 -12  11                 -13  12
   + 7.873192532 10    t   + 6.933755547 10    t   + 3.831183402 10    t 

                   -14  13                 -15  14                 -15  15  
   - 1.380992175 10    t   - 8.070631149 10    t   - 1.495372579 10    t   +

   / 16\\
  O\t  / }
        /
plot( convert(eval(T(t),%),polynom), t=0..20, view=[DEFAULT,0..1000] );

What's the moral to all of this. Don't forget about series solutions, but also don't put too much blind faith in results obtained from a finite number of terms of the series solution.

Hoping this has not been a waste of time and space,

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.ed

I can see both sides of this. There are questions where you do want to insist that the "correct" answer has a particular form and others where any mathematically equivalent expression should be accepted.

As an extreme example, in the differentiation example above, how should MapleTA respond to a student response of:

diff( -cos(2*x - 4)/sin(2*x - 1, x );

Thrown to Toby's generic is( simplify( ... )=0 ); this would be correct.

I know there are simple ways within MapleTA to prevent this, but this is not always clear from the presentation of the problem. If question banks are going to be shared among users, there needs to be some way for instructors to know what Maple is going to do and to modify this as necessary. At present, the only real solution is to dive into the code - TA and Maple. I long for the day when this is no longer the case.

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.ed

 

I believe the place where Maple has trouble with your expressions is sqrt(q0*k0). You want this to be treated as sqrt(-q0)*sqrt(-k0) but Maple is going to look at it as sqrt(q0)*sqrt(k0). Then, with the assumptions on q0 and k0, Maple factors out an I (sqrt(-1)) from each of the two terms; when these are multiplied, the factor of -1 appears in the result.

Now, I do not understand why Maple returns the expressions with sqrt(-k0) and sqrt(-q0) in the following:

simplify(exp(-1/2*ln(-q0)+1/2*ln(-k0)));
                                      (1/2)
                                 (-k0)     
                                 ----------
                                      (1/2)
                                 (-q0)     

Based on my explanation above I would expect sqrt(k0)/sqrt(q0). Maybe Maple does not combine the two square roots because it cannot match the branch cuts. If this is overlooked, then it will create many other situations.

It looks to me as though your simplest solution might be to ask Maple to do symbolic simplification:

simplify(exp(-1/2*ln(-q0)+1/2*ln(-k0))*A,symbolic);
                                       2
                                     k0 
simplify(exp(-1/2*ln(-q0))*exp(1/2*ln(-k0))*A,symbolic);
                                       2
                                     k0 

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.ed

What is your function?

I assume it is:

f := tan( sin(x^2) ) - sin( tan(x^2) );
                            /   / 2\\      /   / 2\\
                         tan\sin\x // - sin\tan\x //

And your question is to find the first time a derivative of f, evaluated at zero, is not zero, right?

How much mathematics do you know? Are you supposed to do this by brute force? Or, can you use something else - like series?

I have the answer to the question, but I don't know which solution to suggest to you - or if I should even post a solution if this is a homework problem. If it is a homework problem, please show us what your attempts with Maple and we'll guide you from there.

Please confirm my understanding of the question, and provide some additional information. I'm sure you'll find some helpful guidance on this site.

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.ed

Acer,

That's part of the mathematics I was leaving for the original poster. This looks like an assignment. I'm happy to answer Maple questions, but not mathematics (at least not yet).

But, this is an interesting function, in a some square sense. (:-))

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.ed

This looks like it could be a homework assignment. I'll provide some guidance on the Maple, but not the overall mathematics.

Your function does not depend on N (n is the summation index and the result will depend only on x). To define a function, you would use:

f := x -> (definition of function, with variable x)

For a summation, you could use:

sum( (-1)^(n+1)/n * sin(n*x), n=1..infinity );

Now, you have to figure out how to modify this to include only the odd terms. When you combine all of this, you'll have your function.

Once you define this function, what are you going to do with it?

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.ed

Have you looked on the web for some ideas? I just did a Google search with keywords "maple cobweb" and found a list that contains a wide variety of resources. I would think that at least one of these would be of use to you. If not, please give a more complete explanation of your needs AND what you have tried. (Information about your level of experience with Maple and the mathematics that you know would also be 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.ed

What you will accept as a proof depends on what you already know.

As has already been pointed out, using l'Hopital's Rule is suspect because that presumes knowledge of the derivative of a^x, or exp(x). There are ways to get this using other facts, e.g., series and inverse function. But, this brings in more questions about what is proven and what is accepted fact for working problems.

Can you provide us with some additional background for your question?

I assume it's for an assignment you have been given. For what course? What do you know from previous courses? Can you provide any similar problems (or examples) that would illuminate us on the general level of your course?

Lastly,  your initial post asked for a "proof using Maple". What does this mean, exactly? Do you simply want to verify the calculations at each step? Or, might a "graphical proof" be acceptable?

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.ed

From what you have written, you have not used Maple at all in this problem. The questions you are asking are about the problem on your assignment. You are supposed to answer these questions on your own, to show your understanding of the material.

I have two comments about your work:

  1. Does it satisfy the original problem? This is one nice fact about differential equations - it's always possible to check your answer.
  2. Check the definitions of terms.

These will enable you to answer most of your questions.

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.ed

I think I see the problem. You really need to have nested if statements. Try this:

if mode=TM then
  if n=0 then Q:=QTM2
         else Q:=QTM1
  end if;
elif mode=TE then
  if   m=0 then Q:=QTE2
  elif n=0 then Q:=QTE3
           else Q:=QTE1
  end if;
end if

There are many other ways to do this, including some that might be slightly more efficient. But, this is the general idea.

I hope this helps,

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.ed

What about using the following:

stop_cond = [ (y(t)-2)*min(t-4,0), (y(t)-2)+min(t-4,0) ]

The first will be zero when either of the two terms is zero. Both will be zero only when both terms are zero. This is based on the fact that the unique solution to { x*y=0, x+y=0 } is { x=0, y=0 }. You could use any homogeneous nonsingulary linear system, but I thought the nonlinear system was actually a little clearer in how this works.

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.ed

Your system of ODEs has a singularity at t=0, so specifying ICs there can be problematic. Are you sure this is the problem you want to solve? You might be able to get away with changing the initial time to a "small" time, but I'll use t=1 just to check that Maple is capable of giving a solution to your system.

While I have not done a full analysis, this system has all the markings of a stiff system. Some coefficients are 1e-10, some are 1e-17, and others are order 1. This is likely to make this problem very difficult to solve.

As for the Maple, the only problem I see is the definition of the initial conditions (ics). You have only one IC, followed by separate equations. Change the semi-colon's to commas, and you will have a full set of ICs.

Making the syntax change and changing the initial time to t=1, Maple will give a numerical solution.

Watching the information from infolevel, I also see that the equations are complex valued? Is this correct?

Here is some output from the numerical solver:

 

eqsol(0.01);
[          
[t = 0.01, 

                                   513                           495    
  Phi(t) = 1.7561641855734908122 10    - 8.8862566109881525003 10    I, 

                                       510                           493    
  Theta0(t) = -2.9293833977001900868 10    + 1.4822789803922300535 10    I, 

                                      506                           489    
  Theta1(t) = 3.1814523334036783304 10    - 1.6098268067697731582 10    I, 

                                     511                           494    
  delta(t) = 2.6364822192224239381 10    - 1.3340698861061906950 10    I, 

                                  489                           506  ]
  v(t) = -4.8295006444616233419 10    - 9.5443969685953777495 10    I]

 

eqsol(0.1);
[                                          44                           26    
[t = 0.1, Phi(t) = 1.6086847331654627340 10   - 8.1400050533538737390 10   I, 

                                       42                           25    
  Theta0(t) = -2.7038561315335253698 10   + 1.3681613382950502379 10   I, 

                                      39                           22    
  Theta1(t) = 3.0043068299750770752 10   - 1.5201905179830742978 10   I, 

                                     43                           26    
  delta(t) = 2.4338163766357651140 10   - 1.2315202100393673506 10   I, 

                                  22                           39  ]
  v(t) = -4.5625739400441621539 10   - 9.0168777453811162121 10   I]

 

eqsol(0.5);
 [                                                                  -15    
 [t = 0.5, Phi(t) = 326.19430771200165313 - 1.6505489550294421852 10    I, 

                                                                -14    
   Theta0(t) = -28.419868131110325123 + 3.4209043642188798897 10    I, 

                                                                -16    
   Theta1(t) = 0.17714921898605595204 - 1.2390596890630632867 10    I, 

   delta(t) = 255.98106937583365131 + 0.0065000000007336814662 I, 

                                                           ]
   v(t) = 0.40000000003012568494 - 0.53612638482822179542 I]

Here's a plot of the real part of each of the 5 components of the solution. Note that I have scaled Theta0 and v to show more detail of the other components. (Call me lazy. I did not re-upload the plot after I changed the labels in the legend.)

odeplot( eqsol, [[t,Re(Phi(t))],[t,Re(Theta0(t))/5],[t,Re(Theta1(t))],[t,Re(delta(t))],[t,Re(v(t))/10]], t=0.5..10, legend=["Phi","Theta0/5","Theta1","delta","v/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.ed

Why do you want to convert this system of ODEs to cartesian coordinates?

As you have presented the system, the two equations are not connected. You can solve the theta equation: theta(t)=t+C1. The r equation is separable, and it's possible to obtain an explicit expression for r(t).

SYS := [ diff(r(t),t)=sin(Pi*r(t)), diff(theta(t),t)=1 ];
SOL := dsolve( SYS, [r(t),theta(t)] );
[                      
[                      
[                      
[{theta(t) = t + _C1}, 
[                      

   /             /   2 exp(Pi t + _C2 Pi)       exp(2 Pi t + 2 _C2 Pi) - 1\\ ]
   |       arctan|--------------------------, - --------------------------|| ]
  <              \1 + exp(2 Pi t + 2 _C2 Pi)    1 + exp(2 Pi t + 2 _C2 Pi)/ >]
   |r(t) = ----------------------------------------------------------------| ]
   \                                      Pi                               /

From here you can find explicit expressions for the Cartesian coordinates

CART := [ x(t)=r(t)*cos(theta(t)),  y(t)=r(t)*sin(theta(t)) ]:
eval( q1, map(op,SOL) );
[       
[       
[       
[x(t) = 
[       

        /   2 exp(Pi t + _C2 Pi)       exp(2 Pi t + 2 _C2 Pi) - 1\               
  arctan|--------------------------, - --------------------------| cos(t + _C1)  
        \1 + exp(2 Pi t + 2 _C2 Pi)    1 + exp(2 Pi t + 2 _C2 Pi)/               
  -----------------------------------------------------------------------------, 
                                       Pi                                        

  y(t) = 

        /   2 exp(Pi t + _C2 Pi)       exp(2 Pi t + 2 _C2 Pi) - 1\             
  arctan|--------------------------, - --------------------------| sin(t + _C1)
        \1 + exp(2 Pi t + 2 _C2 Pi)    1 + exp(2 Pi t + 2 _C2 Pi)/             
  -----------------------------------------------------------------------------
                                       Pi                                      

  ]
  ]
  ]
  ]
  ]

But, if you really want to put these ODEs into cartesian coordinates, here is how I would do it:

POLAR := [ r(t)=sqrt(x(t)^2+y(t)^2), theta(t)=arctan(y(t)/x(t)) ];
dCART := diff( CART, t );
simplify( eval( eval(dCART,SYS), POLAR ), symbolic );
[                 /                  (1/2)\                       (1/2)       
[                 |   /    2       2\     |        /    2       2\            
[ d           -sin\Pi \x(t)  + y(t) /     / x(t) + \x(t)  + y(t) /      y(t)  
[--- x(t) = - --------------------------------------------------------------, 
[ dt                                              (1/2)                       
[                                  /    2       2\                            
[                                  \x(t)  + y(t) /                            

                /                  (1/2)\                       (1/2)     ]
                |   /    2       2\     |        /    2       2\          ]
   d         sin\Pi \x(t)  + y(t) /     / y(t) + \x(t)  + y(t) /      x(t)]
  --- y(t) = -------------------------------------------------------------]
   dt                                           (1/2)                     ]
                                 /    2       2\                          ]
                                 \x(t)  + y(t) /                          ]

 

This looks like a real mess. And, when Maple tries to solve these ODEs, it encounters some difficulties.

dsolve( %, [x(t),y(t)] );
Error, (in dsolve) numeric exception: division by zero

If I have missed something, please do not hesitate to ask again.

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.ed

 

What version of MapleTA are you using? What version of Maple is it accessing? If I had to guess, I'd guess your MapleTA is not accessing Maple 11.

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.ed
First 26 27 28 29 30 31 32 Last Page 28 of 44