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 replies submitted by Doug Meade

Let's not jump to conclusions, too quickly. I have not done the analysis (more on this later), but I do want to issue one caution.

First, the red and blue curves are the nullclines of the system. These are the curves in phase space where the v' (red) and w' (blue) are zero. Intersections of v' and w' nullclines are the equilibrium solutions. Across each v nullcline (red), solutions have no horizontal movement, only up or down. And, across each w nullcline (blue), solutions have no vertical movement, only left and right. Because of the scaling in Robert's plot, the horizontal movement is not visible. It should be possible to see this with an appropriate rescaling of the problem. Does it make sense to put this in a non-dimensional form?

As I admitted in my opening, I have not explored this to see what is really happening in this situation. This was not for a lack of trying. I copied Robert's code, but the DEplot terminated due to an "integer too large in context". Replacing sys100 with evalf(sys100) allowed Maple to work longer, but with fatal results (connection to kernel lost). I do not have more time right now to investigate this further.

I look forward to what others have to say about this system.

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 was not comfortable with my assumption that varphi=t. Your message confirms this. It would be most helpful to see the full set of three ODEs. It's quite possible Maple can be used to create a nice visualization of this system of ODEs.

In the meantime, I want to show you how to create some animations of the slope field and equilibrium solutions for your original equation. We start as before:

 

restart;
with( plots ):
with( DEtools ):
with( plottools ):
param := [ B=1, a=0.3 ];
x_:=diff(u(t),t)=(B/a-B+B*u(t)-varphi)*u(t);
F := rhs(x_);
equilib_soln := solve( F=0, u(t) );

Creating separate animations of the equilibrium solution and slope field is now a simple modification of the code in my previous post:

 

aeqP := animate( plot, [ [eval(equilib_soln,param)], t=-2..5,
                 thickness=2, color=blue], varphi=-2..5, frames=16 ):
aeqP;


asfP := animate( DEplot, [eval(x_,param), u(t), t=-2..5, u=-5..5,
                 dirgrid=[40,40]], varphi=-2..5, frames=16 ):
asfP;

To superimpose the equilibrium solutions on top of the slope field took a little more work than I thought it would. (The problem appears to be the premature evaluation of DEplot before it has a value of varphi.) Here is what I thought would work, and the warning and error that it generates):

 

ppP := animate( display, [ plot([eval(equilib_soln,param)], t=-2..5, thickness=2, color=blue),
                           DEplot(eval(x_,param), u(t), t=-2..5, u=-5..5, dirgrid=[40,40]) ],
                varphi=-2..5, frames=16 );
Warning, unable to evaluate 1 of the 2 functions to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct
Error, (in DEtools/DEplot/CheckDE) extra unknowns found: varphi

Here is the code I wrote to get around this problem. The main trick here is not let DEplot be called before varphi has a numerical value.

 

ppP := proc( varphi )
  local param, x_, F, equilib_soln, eqP, sfP;
  uses plots;
  if not is(evalf(varphi),numeric) then
    return 'procname'(_passed);
  end if;
  param := [ B=1, a=0.3 ];
  x_:=diff(u(t),t)=(B/a-B+B*u(t)-varphi)*u(t);
  F := rhs(x_);
  equilib_soln := solve( F=0, u(t) );
  eqP := plot( [eval(equilib_soln,param)], t=-2..5, thickness=2, color=blue ):
  sfP := DEplot( eval(x_,param), u(t), t=-2..5, u=-5..5, dirgrid=[40,40] ):
  display( eqP, sfP );
end proc:
 
animate( ppP, [varphi], varphi=-2..5, frames=16 );


Unfortunately, MaplePrimes does not support animated gifs. The direct link to the animated GIF file is http://www.mapleprimes.com/files/178_aP3.gif

These pictures clearly indicate that this differential equation is similar to a logistic equation. There are some sign differences that make it different, and the stability of the nonzero equilibrium changes at zero, but maybe varphi is always positive.

While I doubt I need to see the full derivation of the equations, it would be interesting to see the full set of 3 ODEs with relevant values for the parameters.

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 was not comfortable with my assumption that varphi=t. Your message confirms this. It would be most helpful to see the full set of three ODEs. It's quite possible Maple can be used to create a nice visualization of this system of ODEs.

In the meantime, I want to show you how to create some animations of the slope field and equilibrium solutions for your original equation. We start as before:

 

restart;
with( plots ):
with( DEtools ):
with( plottools ):
param := [ B=1, a=0.3 ];
x_:=diff(u(t),t)=(B/a-B+B*u(t)-varphi)*u(t);
F := rhs(x_);
equilib_soln := solve( F=0, u(t) );

Creating separate animations of the equilibrium solution and slope field is now a simple modification of the code in my previous post:

 

aeqP := animate( plot, [ [eval(equilib_soln,param)], t=-2..5,
                 thickness=2, color=blue], varphi=-2..5, frames=16 ):
aeqP;


asfP := animate( DEplot, [eval(x_,param), u(t), t=-2..5, u=-5..5,
                 dirgrid=[40,40]], varphi=-2..5, frames=16 ):
asfP;

To superimpose the equilibrium solutions on top of the slope field took a little more work than I thought it would. (The problem appears to be the premature evaluation of DEplot before it has a value of varphi.) Here is what I thought would work, and the warning and error that it generates):

 

ppP := animate( display, [ plot([eval(equilib_soln,param)], t=-2..5, thickness=2, color=blue),
                           DEplot(eval(x_,param), u(t), t=-2..5, u=-5..5, dirgrid=[40,40]) ],
                varphi=-2..5, frames=16 );
Warning, unable to evaluate 1 of the 2 functions to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct
Error, (in DEtools/DEplot/CheckDE) extra unknowns found: varphi

Here is the code I wrote to get around this problem. The main trick here is not let DEplot be called before varphi has a numerical value.

 

ppP := proc( varphi )
  local param, x_, F, equilib_soln, eqP, sfP;
  uses plots;
  if not is(evalf(varphi),numeric) then
    return 'procname'(_passed);
  end if;
  param := [ B=1, a=0.3 ];
  x_:=diff(u(t),t)=(B/a-B+B*u(t)-varphi)*u(t);
  F := rhs(x_);
  equilib_soln := solve( F=0, u(t) );
  eqP := plot( [eval(equilib_soln,param)], t=-2..5, thickness=2, color=blue ):
  sfP := DEplot( eval(x_,param), u(t), t=-2..5, u=-5..5, dirgrid=[40,40] ):
  display( eqP, sfP );
end proc:
 
animate( ppP, [varphi], varphi=-2..5, frames=16 );


Unfortunately, MaplePrimes does not support animated gifs. The direct link to the animated GIF file is http://www.mapleprimes.com/files/178_aP3.gif

These pictures clearly indicate that this differential equation is similar to a logistic equation. There are some sign differences that make it different, and the stability of the nonzero equilibrium changes at zero, but maybe varphi is always positive.

While I doubt I need to see the full derivation of the equations, it would be interesting to see the full set of 3 ODEs with relevant values for the parameters.

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 there might be some confusion with terminology.

When I think of isocline, I think of a curve where the slopes are constant. Plotting isoclines would be plotting level curves of the slope (in (t,k) space). This is doable. I'll bet someone will have it done before I get back from my afternoon meetings.

Keep telling us exactly what you want and there's a good chance we can 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.ed

I think there might be some confusion with terminology.

When I think of isocline, I think of a curve where the slopes are constant. Plotting isoclines would be plotting level curves of the slope (in (t,k) space). This is doable. I'll bet someone will have it done before I get back from my afternoon meetings.

Keep telling us exactly what you want and there's a good chance we can 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.ed

I agree with Robert. It is unusual to talk about a plot of the phase plane for a single ODE. You can ask for the slope field, but this does not show k_dot vs k. The only plot of k_dot vs. k is the one in your original post:

y:=A*k^a:
S:=s*y:
BE:=(n+g+d)*k:
k_dot:=S-BE;
s:=0.7: A:=1: a:=0.3: n:=0.06: g:=0.04: d:=0.03:
plot( k_dot, k=0..20, labels=["k","k_dot"] );

The slope field for this ODE can be produced with DEplot as follows:

with( DEtools ):
ODE := diff( k(t),t ) = eval( 100*k_dot, k=k(t) ):

Note that I rescaled the RHS by a factor of 100. This is equivalent to changing the time scale by a factor of 100. Without this adjustment, all of the slopes appear to be very close to zero.

Since this ODE is autonomous, the slopes are the same at every instant in time. So, maybe you want to see only one column of slopes. Here is how I did this

DEplot( ODE, k(t), t=0..1, k=0..20, arrows=line,
        dirgrid=[2,20],
        view=[-0.1..0.1,DEFAULT] );

I hope something in here is useful 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.ed

I agree with Robert. It is unusual to talk about a plot of the phase plane for a single ODE. You can ask for the slope field, but this does not show k_dot vs k. The only plot of k_dot vs. k is the one in your original post:

y:=A*k^a:
S:=s*y:
BE:=(n+g+d)*k:
k_dot:=S-BE;
s:=0.7: A:=1: a:=0.3: n:=0.06: g:=0.04: d:=0.03:
plot( k_dot, k=0..20, labels=["k","k_dot"] );

The slope field for this ODE can be produced with DEplot as follows:

with( DEtools ):
ODE := diff( k(t),t ) = eval( 100*k_dot, k=k(t) ):

Note that I rescaled the RHS by a factor of 100. This is equivalent to changing the time scale by a factor of 100. Without this adjustment, all of the slopes appear to be very close to zero.

Since this ODE is autonomous, the slopes are the same at every instant in time. So, maybe you want to see only one column of slopes. Here is how I did this

DEplot( ODE, k(t), t=0..1, k=0..20, arrows=line,
        dirgrid=[2,20],
        view=[-0.1..0.1,DEFAULT] );

I hope something in here is useful 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.ed

Maple is not going to take the place of thinking. There is no magic command that will give you the functions that you need in your Questions 1 and 2. You have to do some mathematical thinking. Only after you have done this can Maple be used as a tool to show that your ideas are sound, or to point out errors in your thinking.

Personally, I think your first stop should be your instructor.

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

Maple is not going to take the place of thinking. There is no magic command that will give you the functions that you need in your Questions 1 and 2. You have to do some mathematical thinking. Only after you have done this can Maple be used as a tool to show that your ideas are sound, or to point out errors in your thinking.

Personally, I think your first stop should be your instructor.

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 do you know about Maple? Any reasonable instructor will have shown you some similar examples or given some instruction on how to use Maple. How is Maple being used in your course? Have you asked your instructor for help?

There is no infinity key. In Maple, you just type out infinity, or click on the infinity symbol in a pallette. (What version of Maple are you using?)

Why do you believe the exponential funciton is appropriate for Question 1? In Maple you would enter the exponential function as exp(x). What are the critical points of the exponential function?

What mathematical property does the function in Question 2 have to satisfy? (It's not sufficient to simply say "I don't know". This is not a treasure hunt, it's a learning exercise.)

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 do you know about Maple? Any reasonable instructor will have shown you some similar examples or given some instruction on how to use Maple. How is Maple being used in your course? Have you asked your instructor for help?

There is no infinity key. In Maple, you just type out infinity, or click on the infinity symbol in a pallette. (What version of Maple are you using?)

Why do you believe the exponential funciton is appropriate for Question 1? In Maple you would enter the exponential function as exp(x). What are the critical points of the exponential function?

What mathematical property does the function in Question 2 have to satisfy? (It's not sufficient to simply say "I don't know". This is not a treasure hunt, it's a learning exercise.)

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

Software change requests (SCRs) include bugs.

You can submit a Maple SCR through MaplePrimes. Just click on the link in the list under Create post on the left hand side of the MaplePrimes window.

How does one put a symbol like ∈ in the subject?

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

Software change requests (SCRs) include bugs.

You can submit a Maple SCR through MaplePrimes. Just click on the link in the list under Create post on the left hand side of the MaplePrimes window.

How does one put a symbol like ∈ in the subject?

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 remember working this problem. It is a little tricky, until you catch on to the essentials of this approach.

The important part of this is that you are looking for the curves that are orthogonal to each of the circles centered at (c,0) with radius c.

The key is to find an expression for dy/dx that does not depend on the parameter c. If you simply differentiate y^2+(x-c)^2=c^2 you will find dy/dx, but only in terms involving c:

 

eq1 := y^2+(x-c)^2=c^2;

                              2          2    2
                             y  + (x - c)  = c 

implicitdiff( eq1, y, x );

                                     x - c
                                   - -----
                                       y  

But, note that expanding the original description of the circles leaves only one occurrence of c, which can then be solved for.

 

q1 := isolate( eq1, c );

                                       2    2
                                     -y  - x 
                               c = - --------
                                       2 x   

Now, use implicit differentiation to find an expression for the slopes of these curves:

 

q2 := implicitdiff( q1, y, x );

                                     2    2
                                    x  - y 
                                  - -------
                                     2 x y 

The orthogonal trajectories have a slope that is the negative reciprocal of this:

 

eq2 := -1/-(1/2)*(x^2-y^2)/(x*y);

                                    2 x y 
                                   -------
                                    2    2
                                   x  - y 

So, we are looking for the solutions to

 

de2 := diff( y(x), x ) = eval( -1/q2, y=y(x) );

                             d          2 x y(x) 
                            --- y(x) = ----------
                             dx         2       2
                                       x  - y(x) 

Here are some of the steps in finding its solutions:

 

q3 := dsolve( de2, y(x), implicit );
                              / 2       2\                  
                   /y(x)\     |x  + y(x) |                  
                 ln|----| - ln|----------| - ln(x) - _C1 = 0
                   \ x  /     |     2    |                  
                              \    x     /  
                
q4 := simplify( q3, symbolic );
                                  / 2       2\          
                     ln(y(x)) - ln\x  + y(x) / - _C1 = 0

q5 := combine( map(exp,eval(q4,_C1=-ln(2*C))) );
                                2 y(x) C     
                               ---------- = 1
                                2       2    
                               x  + y(x)     

q6 := q5*denom(lhs(q5));
                                        2       2
                            2 y(x) C = x  + y(x) 

q7 := student[completesquare](rhs(q6)-lhs(q6), y(x)) = 0;
                                    2    2    2    
                          (y(x) - C)  - C  + x  = 0

q7+C^2;
                                      2    2    2
                            (y(x) - C)  + x  = C 

There is one more curve to add to this family of curves. The constant function y=0 obviously solves the differential equation for the orthogonal trajectories (de2).

So, the orthogonal trajectories are the circles centered at (0,c) with radius c and the horizontal line y=0 (x-axis).

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 remember working this problem. It is a little tricky, until you catch on to the essentials of this approach.

The important part of this is that you are looking for the curves that are orthogonal to each of the circles centered at (c,0) with radius c.

The key is to find an expression for dy/dx that does not depend on the parameter c. If you simply differentiate y^2+(x-c)^2=c^2 you will find dy/dx, but only in terms involving c:

 

eq1 := y^2+(x-c)^2=c^2;

                              2          2    2
                             y  + (x - c)  = c 

implicitdiff( eq1, y, x );

                                     x - c
                                   - -----
                                       y  

But, note that expanding the original description of the circles leaves only one occurrence of c, which can then be solved for.

 

q1 := isolate( eq1, c );

                                       2    2
                                     -y  - x 
                               c = - --------
                                       2 x   

Now, use implicit differentiation to find an expression for the slopes of these curves:

 

q2 := implicitdiff( q1, y, x );

                                     2    2
                                    x  - y 
                                  - -------
                                     2 x y 

The orthogonal trajectories have a slope that is the negative reciprocal of this:

 

eq2 := -1/-(1/2)*(x^2-y^2)/(x*y);

                                    2 x y 
                                   -------
                                    2    2
                                   x  - y 

So, we are looking for the solutions to

 

de2 := diff( y(x), x ) = eval( -1/q2, y=y(x) );

                             d          2 x y(x) 
                            --- y(x) = ----------
                             dx         2       2
                                       x  - y(x) 

Here are some of the steps in finding its solutions:

 

q3 := dsolve( de2, y(x), implicit );
                              / 2       2\                  
                   /y(x)\     |x  + y(x) |                  
                 ln|----| - ln|----------| - ln(x) - _C1 = 0
                   \ x  /     |     2    |                  
                              \    x     /  
                
q4 := simplify( q3, symbolic );
                                  / 2       2\          
                     ln(y(x)) - ln\x  + y(x) / - _C1 = 0

q5 := combine( map(exp,eval(q4,_C1=-ln(2*C))) );
                                2 y(x) C     
                               ---------- = 1
                                2       2    
                               x  + y(x)     

q6 := q5*denom(lhs(q5));
                                        2       2
                            2 y(x) C = x  + y(x) 

q7 := student[completesquare](rhs(q6)-lhs(q6), y(x)) = 0;
                                    2    2    2    
                          (y(x) - C)  - C  + x  = 0

q7+C^2;
                                      2    2    2
                            (y(x) - C)  + x  = C 

There is one more curve to add to this family of curves. The constant function y=0 obviously solves the differential equation for the orthogonal trajectories (de2).

So, the orthogonal trajectories are the circles centered at (0,c) with radius c and the horizontal line y=0 (x-axis).

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 43 44 45 46 47 48 49 Last Page 45 of 76