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

I think I have something you will like.

First, the "isoclines" you refer to are the equilibrium solutions of the differential equation.

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) );

 

These solutions can be plotted, and superimposed on top of the direction field as follows. (You do not tell us much about varphi. For this I have assumed that varphi is another name for the independent variable (t). If varphi is an actual parameter in the problem, then a different approach will need to be found. Please give some more information, if necessary.)


eqP := plot( [eval(equilib_soln,param)], varphi=-2..5, thickness=2, color=blue ):
sfP := DEplot( eval(x_,[varphi=t,param[]]), u(t), t=-2..5, u=-5..5, dirgrid=[40,40] ):
P := display( eqP, sfP ):
P;

This plot has the axes reversed from your request. This can be rectified by using the transform command from the plottools package.

TR := transform( (x,y)->[y,x] ):
display( TR( display( eqP,sfP ) ), labels=[u,varphi] );

 

This plot has the general properties you requested.

If this is not what you are seeking, might you be interested in a bifurcation diagram with varphi as the parameter? What values are reasonable for varphi?

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

You can use the assume facility. For complete information, see the online help. Here is how I would approach your problem.

assume( a::AndProp(NonZero,constant) );
about( a );
Originally a, renamed a~:
  is assumed to be: AndProp(constant,Non(0))

is(a,constant);
                                    true
is(a=0);
                                    false
is(a^2>0);
                                    FAIL

The last test failed because Maple does not assume a is real-valued. Here is how that could be added:


assume( a::AndProp(NonZero,constant,real) );
is(a^2>0);
                                    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.ed

While I can see how Maple can be used in this assignment, it is not where I would start.

Questions 1 and 2 ask you to find a function with certain properties. These functions are not just going to appear out of a Maple command. You have to have some ideas. What functions are you thinking about using? What would have to happen if a critical point is NOT an extreme point?

Questions 3 and 4 are somewhat easier, at least from the point of view that they involve given functions. Have you plotted the functions? What is n in the interval in Question 4?

f := 2^x * sin(x);
plot( f, x=0..10 );

g := 1/x + x^2;
plot( g, x=0..2, y=0..10 );
plot( g, x=0..3, y=0..10 );

Now that you have a picture of these functions, what properties of the function characterize the smallest and largest values of these functions?

If, after you answer these questions, you still need some help with the Maple, please let us know what you have done and what still needs to be done and I think you'll find some people eager to hellp you with this.

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

You will probably get the following two answers to your question:

  1. Use a worksheet (with 1D Maple input)
  2. Prepare your code in a text file and read the file into Maple.

If you want to have everything in one file, and you won't need to re-use it in other instances, I recommend the first option. For larger projects, or for code that will be used repeatedly, the extra complexity of multiple files can be preferrable.

This post is a self-fulfilling prophecy. You now have two options. If you give us some details about what you are trying to do we can probably give you more precise 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.ed

Please complete your question.

Most Maple users do not do any real programming with Maple. Those that do, usually start by learning to use Maple interactively. As they see a task that needs to be repeated frequently then they may do some programming.

What are you trying to do? What type of mathematics are you needing to use? Do you want symbolic results or are you trying to create graphics?

As we learn a little more about your specific needs we will be better able to direct you to appropriate resources and information.

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 do not know if this is possible or not. I have never called a separate maplet from within another maplet. It's not that difficult to define a window and to use it in multiple maplets.

Maybe something like the following:

 

MakeWindow := proc( i::posint, N::posint, image::string )
  local M, nxt, prv, W;
  uses Maplets[Elements];
  if i>N then error sprintf("invalid arguments, received i=%a and N=%a", i, N) end if;
  nxt := (i mod N) + 1;
  prv := `if`(i=1,N,i-1);
  W := map2( nprintf, "W%a", [i,nxt,prv] );
  M := Window[W[1]]('title' = sprintf("Maplet #%a",i),
         BoxLayout(
           BoxColumn(
             BoxRow("Fill this space with the contents of your choosing."),
             BoxRow(Label(Image(image))),
             BoxRow(
               Button("Next", Action(CloseWindow(W[1]),RunWindow(W[2]) ) ),
               Button("Prev", Action(CloseWindow(W[1]),RunWindow(W[3]) ) ),
               Button("Exit", Shutdown())
             )
           )
         )
       )
end proc:
maplet := Maplet('onstartup' = RunWindow('W1'),
  MakeWindow(1,3, "C:\\image1.jpg"),
  MakeWindow(2,3, "C:\\image2.jpg"),
  MakeWindow(3,3, "C:\\image3.jpg")
):

maplet1 := Maplet('onstartup' = RunWindow('W1'),
  MakeWindow(1,3, "C:\\image1.jpg")
):

maplet2 := Maplet('onstartup' = RunWindow('W1'),
  MakeWindow(1,3, "C:\\image2.jpg")
):

Maplets[Display](maplet1);
Maplets[Display](maplet1);


This is not what you asked for, directly, but it provides the same functionality, right?

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 knew I was tired yesterday, but I'm more tired today, but I do have to admit that I now do see the denominator of the exponent.

I'll bet the only people who would recognize that denominator as ell are those who have made this mistake in the past.

I apologize for my temporary blindness, but am glad you have been able to answer your question.

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

In your post you give different DEs in your text and in your Maple code. Also, your code is not correct (the dsolve is missing in your definition of differentialsolve. When I insert it I do get the reported response about non-convergence of Newton's Method.

Looking closer I see that you are using two different ODEs. You first write about a linear ODE ( a y'' - y' = 0 ) but your code is for a nonlinear ODE (a y'' - y y' = 0 ). Which one are you really interested in?

Maple has no trouble with the linear BVP:

singularperturbation:={0.01*diff(y(t),[t$2]) - diff(y(t),t) = 0, y(-1)=1,y(1)=-1};
differentialsolve:=dsolve(singularperturbation,y(t),type=numeric,method=bvp[midrich],abserr=1e-2);
plots[odeplot](differentialsolve,[t,y(t)],-1..1,color=green,style=line,symbol=circle,symbolsize=10);

Let's look a little closer at the nonlinear equation. This one can, in principle, be solved by hand. The key is to note that ( a  y' - y^2/2 )' = a y'' - y y' = 0 so that a y' = y^2 / 2 + K. This is a  Bernoulli equation that is not too difficult to solve. The key is then to find values of the constants that satisfy the boundary conditions. This is where there are problems, which I will now illustrate:

ODE2 := theta*diff(y(t),t) - y(t)^2/2 = K;
                             / d      \   1     2    
                       theta |--- y(t)| - - y(t)  = K
                             \ dt     /   2          
SOL2 := dsolve( {ODE2,y(-1)=alpha}, y(t) );
           //                                       /       (1/2)\\       \ 
           ||   (1/2)    (1/2)          (1/2)       |alpha 2     ||  (1/2)| 
           ||t K      + K      + theta 2      arctan|------------|| 2     | 
           ||                                       |     (1/2)  ||       | 
           |\                                       \  2 K       //       | 
 y(t) = tan|--------------------------------------------------------------| 
           \                           2 theta                            / 

    (1/2)  (1/2)
   K      2     
C2 := eval( SOL2, [t=1,y(t)=beta] );
          //                              /       (1/2)\\       \              
          ||   (1/2)          (1/2)       |alpha 2     ||  (1/2)|              
          ||2 K      + theta 2      arctan|------------|| 2     |              
          ||                              |     (1/2)  ||       |              
          |\                              \  2 K       //       |  (1/2)  (1/2)
beta = tan|-----------------------------------------------------| K      2     
          \                       2 theta                       /              
solve( C2, K );
          /                                                      (1/2)     
  1       |                               2      2   /  2      2\          
  - RootOf\tan(_Z) beta alpha + tan(_Z) _Z  theta  - \_Z  theta /      beta
  2                                                                        

                   (1/2)      \         
       /  2      2\           |        2
     + \_Z  theta /      alpha/^2 theta 

 

This is the end of the analytic solution approach. Let's give values to the parameters and see what we can find numerically and graphically.

param := [
 theta=0.01,
 alpha=1,
 beta=-1,
 NULL
]:
C3 := eval( C2, param );
        /            /                             /  (1/2) \\       \        
        |            |   (1/2)         (1/2)       | 2      ||  (1/2)|  (1/2) 
-1 = tan|50.00000000 |2 K      + 0.01 2      arctan|--------|| 2     | K      
        |            |                             |   (1/2)||       |        
        \            \                             \2 K     //       /        

   (1/2)
  2     
plot( [rhs(C3),lhs(C3)], K=-0.1..0.1, view=[DEFAULT,-2..2] );


 

So, there are lots of solutions to this equation. This is probably why dsolve had so much trouble finding a  numerical solution.

Let's try to find the smallest (positive) value of K that works:

fsolve( C3, K );
                                0.9818497427
KK := fsolve( C3, K=0..0.001 );
                               0.0005138206214

This gives us the solution as:


SOL3 := eval( SOL2, param ):
SOL4 := eval( rhs(SOL3), K=KK );
/ /
0.02266761173 tan\50.00000000 \0.02266761173 t + 0.02266761173

(1/2) / (1/2)\\ (1/2)\ (1/2)
+ 0.01 2 arctan\22.05790385 2 // 2 / 2

Let's check how this solutions satisfies the boundary condtions:

evalf( eval( SOL4, t=-1 ) );
                                0.9999999976
evalf( eval( SOL4, t=1 ) );
                                -1.000000034

plot( SOL4, t=-1..1, view=[DEFAULT,-2..2] );

 This gives a lot more information about the problem. I hope some of it is relevant, or is easily adapted to 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.ed

This has the definite feel of a homework problem. Many of us are willing to help, but not to do your work for you.

How are you supposed to access Newton's method? Through Maple's built-in command, an implementation you have been given, or an implementation that you are supposed to write on your own?

How much Maple do you know? How fancy a solution is expected? Show us how you started the problem and I'm confident you will get some feedback helping YOU to solve this problem.

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 is used in the nprintf command. nprintf is like C's printf, with the results a Maple "name", which is what is needed to refer to the windows in the maplet. Maple also has printf, sprintf (to produce a string) and fprintf (to send the result to a file). See ?printf for help on each of these.

The best way to see what this is doing is to test it on its own:

i := 1;
N := 3;
                                      1
                                      3
nxt := (i mod N) + 1;
prv := `if`(i=1,N,i-1);
                                      2
                                      3
W := map2( nprintf, "W%a", [i,nxt,prv] );
                                [W1, W2, W3]
i := 3;
                                      3
nxt := (i mod N) + 1;
prv := `if`(i=1,N,i-1);
                                      1
                                      2
W := map2( nprintf, "W%a", [i,nxt,prv] );
                                [W3, W1, W2]


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

 

Are you saying you want to connect the endpoints? There are a number of ways to do this. I won't offer a solution to this because I know Robert Israel will do this better and faster than I can - probably by grabbing the first and last points out of the plot structure.

This is a very different question from what you first requested. The plot is no longer the graph of a function. What if the function enters and leaves your range several times as the angle runs through its domain?

Do you want the endpoints connected by a straight line or by a radial curve?

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 can be done. How it is done depends on what you know and exactly what you find acceptable. I'll start with the code from your post:

restart;
with( plots ):
r := a -> 4*sin(a);
plot( r(a), a=0..Pi, coords=polar );

My approach is to manually truncate the function when the radius is outside a given interval, here [2,4]. I add the scaling and view options to improve the appearance of the plot.

plot( `if`(r(a)>2 and r(a)<4,r(a),undefined), a=0..Pi, coords=polar,
      scaling=constrained, view=[-2..2,0..4] );

It would be nice if you could specify view=[r0..r1,a0..a1] and have this interpretted as a polar window. While you can't (at present), I think the code I will present next shows how this could be implemented.

trim := proc(f,r::range)
  unapply( `if`(f(x)>lhs(r) and f(x)<rhs(r),f(x),undefined), x )
end proc:

r2 := trim(r,2..4);
plot( r2(a), a=0..Pi, coords=polar,
      scaling=constrained, view=[-2..2,0..4] );

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

 

Does this do what you want?

restart;
with( Maplets[Elements] ):
MakeWindow := proc( i::posint, N::posint, image::string )
  local M, nxt, prv, W;
  uses Maplets[Elements];
  if i>N then error sprintf("invalid arguments, received i=%a and N=%a", i, N) end if;
  nxt := (i mod N) + 1;
  prv := `if`(i=1,N,i-1);
  W := map2( nprintf, "W%a", [i,nxt,prv] );
  M := Window[W[1]]('title' = sprintf("Maplet #%a",i),
         BoxLayout(
           BoxColumn(
             BoxRow("Fill this space with the contents of your choosing."),
             BoxRow(Label(Image(image))),
             BoxRow(
               Button("Next", Action(CloseWindow(W[1]),RunWindow(W[2]) ) ),
               Button("Prev", Action(CloseWindow(W[1]),RunWindow(W[3]) ) ),
               Button("Exit", Shutdown())
             )
           )
         )
       )
end proc:
maplet := Maplet('onstartup' = RunWindow('W1'),
  MakeWindow(1,3, "C:\\image1.jpg"),
  MakeWindow(2,3, "C:\\image2.jpg"),
  MakeWindow(3,3, "C:\\image3.jpg")
):
Maplets[Display](maplet);

If the images do not exist, Maple does not complain and leaves an empty label.  You will need to modify the paths to the images to fit your needs.

The only oddity in this maplet is that the title bar turns grey for all but the initial window. Maybe I have missed something, but this looks like a minor bug to me.

I have uploaded the worksheet to MaplePrimes:
View 178_MultiWindowMaplet.mw on MapleNet
or Download 178_MultiWindowMaplet.mw
View file details

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

I don't think you can get the animation to display in WordPad.

But, if you export the animation as a GIF file, you will see the animation.

Unfortunately, MaplePrimes cannot display the animation. But, click here to see the animation from my previous post.

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 do not see any problems with this. This creates an animation of data that you generate.

Is your question about the animation or about the computations of the data that you plot?

It appears to me that you are using two different models (logistic and exponential) in the loop. You compute values for each model but, since you use the same name for both, are plotting only the last one computed (exponential).

There is a problem in the logistic model -  your formula has both k and K (both should be K).

Would you like to have the solutions for both models plotted in the same frame? If so, here is one way to do this

 

restart: with(plots):
N := Array( 1..2, 0..31 ):
plot_list1 := Array( 0..31 ):
plot_list2 := Array( 0..31 ):
N[1,0]:= 10;
N[2,0]:= 10;  K:= 1000 ; # initial condition and `carrying capacity`
title2:=`N(t) vs t when r= `: # title blanks for plots
sample_r:= [0.4, 0.7, 1.0, 1.3, 1.6, 1.9];
P := Array( 1..nops(sample_r) ):
for i from 1 to nops(sample_r)  do
 r := sample_r[i];
 plot_list1[0]:= [0, N[1,0]]:
 plot_list2[0]:= [0, N[2,0]]:

 for t from 0 to 30 do
  N[1,t+1]:= N[1,t] + ( r / K ) * N[1,t] * ( K - N[1,t] ):  # logistic model
  N[2,t+1]:= N[2,t] * exp(r * (1 - N[2,t] / K)):
  plot_list1[t+1]:= [t+1, N[1,t+1]]:
  plot_list2[t+1]:= [t+1, N[2,t+1]]:
 od;

 P[i] := plot( [convert(plot_list1,list),convert(plot_list2,list)],
               style=line, color=[plum,cyan],
               title=sprintf("%s%a",title2,r)):
od:
display( P );

display( convert(P,list), insequence=true );

One of the big changes in my code is the replacement of all dynamically generated lists with explicitly declared arrays. To keep track of both models you could use a second plot_list (as I have done) or an two-dimensional array (as I do for keeping track of the populations). Note that plot and display are happier when they receive a list (compared to an array).

If I were doing this problem from the beginning, I would approach the problem quite differently (using seq, etc.) but I have tried to use your solution as the basis for this response.

I hope you find this 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.ed
First 29 30 31 32 33 34 35 Last Page 31 of 44