Constructing a contour plot is difficult.

Numerically, it is difficult to determine if two points are on the same or different contours. I believe Maple's approach does not use an adaptive method. A grid is placed on the domain and points on this grid are assigned to different contours. There are a few things you can do to help Maple: give an explicit collection of contours to be plotted and to specify the size of the grid to be used. For example,

with(plots):
ODE := t*diff( x(t),t ) = x(t);
ISO := solve( ODE, diff(x(t),t) );
contourplot( ISO, t=-10..10, x=-10..10,
             contours=[-2,-1,0,1,2], grid=[51,51] );

Compare this with

contourplot( ISO, t=-10..10, x=-10..10,
             contours=[-2,-1,0,1,2], grid=[50,50] );

The drastic difference is caused by the fact that the points on both axes are grid points in the first case and not the second.

You might also be interested in the contourplot3d command that produces a 3D plot using the same syntax.

contourplot3d( ISO, t=-10..10, x=-10..10,
               contours=[-2,-1,0,1,2], grid=[51,51] );

or, to give the usual 2D perspective, use

contourplot3d( ISO, t=-10..10, x=-10..10,
               contours=[-2,-1,0,1,2], grid=[51,51],
               axes=normal, orientation=[-90,180] );

I hope this sheds some light on 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.edu/~meade/

Robert's suggestion does the job, but only for plotting individual points (which is likely to be the most common usage). Here are two variations using adaptive=false and either sample= or numpoints=. That these produce significantly different results is not expected.

integerplot2:= proc(f,r::name=range )
    local x, a,b;
    x := op(1,r); a,b:= op(op(2,r));
    plot(f, x = ceil(a) .. floor(b), numpoints=floor(b)-ceil(a)+1,
            adaptive=false, _passed[3..-1]);
 end proc:
integerplot3:= proc(f,r::name=range)
    local x, a,b;
    x := op(1,r); a,b:= op(op(2,r));
    plot(f, x = a .. b, sample=[$ceil(a)..floor(b)],
            adaptive=false, _passed[3..-1]);
 end proc:

These can be used as follows:

integerplot2( cos(Pi*x/2), x=0..100/Pi, color=blue, symbol=box );
integerplot2( cos(Pi*x/2), x=0..100/Pi, style=point, color=blue, symbol=box );
integerplot3( cos(Pi*x/2), x=0..100/Pi, color=blue, symbol=box );

I was surprised to find that the automatic sample in integerplot2 did NOT use exactly evenly spaced points.

My additional $0.50 (2x$0.25),

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/

This is a fairly common request. Over the years, the options to the plot commands have expanded to give users much greater control over their plots. Unfortunately, in my opinion, most users (including myself) have not kept up to date on all of the available options.

When i saw your post I decided to take a look at the on line help for the options to the plot command (?plot,options). From this long list, I found numpoints (old), adaptive, and sample. The descriptions of numpoints and sample include reference to the "minimum number of points" and "initial sampling" which do not sound promising. But, reading further, the description of sample mentions using adaptive=false to strictly adhere to the user-provided sample points.

Even though the description of numpoints does not make explicit reference to adaptive=false, I decided to give it a try.

plot( x^2, x=-5..5, numpoints=11, adaptive=false, style=point );

This appears to work as you have requested.

It would be a fairly simple task to create a procedure that converted the domain into one with integer endpoints and automatically generated the appropriate numpoints= option. (I have to go teach, so cannot do this now. I'll bet you or someone else will have this done by the time I get back to MaplePrimes.)

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/

While your expressions do look complicated, you should not give up so quickly on an explicit evaluation, at least for ETa.

ETa:= 9.285051068*Int(9.285051068*min(.5676676416*Pnm, 50.)*exp(-9.285051068*Pnm/Pna)/Pna, Pnm = 0. .. infinity);
                /  /infinity)
                | |                 
                | |                 
    9.285051068 | |                 
                |/                  
                \ 0.                

                                                /  9.285051068 Pnm\     \
      9.285051068 min(0.5676676416 Pnm, 50.) exp|- ---------------|     |
                                                \        Pna      /     |
      ------------------------------------------------------------- dPnm|
                                   Pna                                  |
                                                                        /

This integrand is, for each given Pna, piecewise linear times an exponential, which can be evaluated explicitly. The key here is to split the integral into two pieces based on the part of the min that is to be used. Maple can write the integrand in piecewise form:

convert( ETa, piecewise, Pnm );
                   /  /infinity                 /                    
                   | |                          |                    
                   | |                          |                    
       9.285051068 | |                 piecewise|Pnm <= 88.07970780, 
                   |/                           \                    
                   \ 0.                                              

                        /  9.285051068 Pnm\                         
         5.270823042 exp|- ---------------| Pnm                     
                        \        Pna      /                         
         --------------------------------------, 88.07970780 < Pnm, 
                          Pna                                       

                        /  9.285051068 Pnm\\     \
         464.2525534 exp|- ---------------||     |
                        \        Pna      /|     |
         ----------------------------------| dPnm|
                        Pna                /     |
                                                 /

but there is no easy way to get Maple to write this as the sum of two integrals with different integrands. So, we do this by hand:

ETa2 := 9.285051068*(Int(9.285051068*.5676676416*Pnm*exp(-9.285051068*Pnm/Pna)/Pna, Pnm = 0. .. 50/0.5676676416)
                    +Int(9.285051068*50.            *exp(-9.285051068*Pnm/Pna)/Pna, Pnm = 50/0.5676676416 .. Finfinity));
            /  /88.07970780                /  9.285051068 Pnm\         \
            | |             5.270823042 exp|- ---------------| Pnm     |
            | |                            \        Pna      /         |
9.285051068 | |             -------------------------------------- dPnm|
            |/                               Pna                       |
            \ 0.                                                       /

                 /  /infinity                       /  9.285051068 Pnm\     \
                 | |                 464.2525534 exp|- ---------------|     |
                 | |                                \        Pna      /     |
   + 9.285051068 | |                 ---------------------------------- dPnm|
                 |/                                 Pna                     |
                 \ 88.07970780                                              /

These integrals can be evaluated explicitly:

value( ETa2 );
                                     /  817.8245844\    
  0.5676676416 Pna - 0.5676676416 exp|- -----------| Pna
                                     \      Pna    /    

                      /  817.8245844\               /     /
     - 464.2525531 exp|- -----------| + 9.285051068 |limit|
                      \      Pna    /               \     \

           /  2321262767 Pnm\         /  76952181988817 \                \\
    -50 exp|- --------------| + 50 exp|- ---------------|, Pnm = infinity||
           \  250000000 Pna /         \  94093750000 Pna/                //

Well, Maple does need to be told that Pna is positive (right?)

value( ETa2 ) assuming Pna>0;
                    -12 /              10        /817.8245844\
      6.033000509 10    |9.409375000 10   Pna exp|-----------|
                        \                        \    Pna    /

                         10                     13\    /  817.8245844\
         - 9.409375000 10   Pna - 7.695218199 10  | exp|- -----------|
                                                  /    \      Pna    /

                        464.2525534              
         + --------------------------------------
                     (76952181988817/94093750000)
           /   / 1 \\                            
           |exp|---||                            
           \   \Pna//                   

And, now, it's not hard to plot this expressiov vs. Pna.

plot( ETa2, Pna=0..100 );

Is this close to what you are trying to do?

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/

Try this.

restart;
eq1 := tan(beta)=-2/alpha*beta;
                                          2 beta
                            tan(beta) = - ------
                                          alpha 
eq2 := beta = sqrt(lambda*alpha-(alpha/2)^2);
                                                     (1/2)
                          1 /                      2\     
                   beta = - \4 lambda alpha - alpha /     
                          2                               

The first equation can be solved for alpha, and this can be substituted into the first to obtain:

eq3 := eval( eq2, isolate(eq1,alpha) );
                                                       (1/2)
                         /                         2  \     
                       1 |  8 lambda beta    4 beta   |     
                beta = - |- ------------- - ----------|     
                       2 |    tan(beta)              2|     
                         \                  tan(beta) /     

Solving for lambda is now straightforward:

eq4 := solve( eq3, lambda );
                                  /         2    \
                             beta \tan(beta)  + 1/
                           - ---------------------
                                  2 tan(beta)     
eq5 := simplify( eq4, trig );
                                     beta         
                           - ---------------------
                             2 cos(beta) sin(beta)
plot( eq5, beta=-Pi/2..Pi/2, view=[DEFAULT,-10..0] );

---------------------------------------------------------------------
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/

Why do you believe the plot is of k vs. P_bar?

If, in the definition of curve, you remove the definition of P, the plot is different. Create a copy of curve, say curve2, but remove the definition of P_bar. Then, plot both curve and curve2. Here's the command I used:

plot([curve, curve2], 0. .. .5, view = [0 .. .1, 0 .. 9000],
  axes = boxed, gridlines = true, title = "P - k RELATION",
  labels = ["k", "P"]);

The shapes are similar, but the plotted points are clearly different.

 

Am I missing something?

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 believe this question has been discussed previously. Check out . My suggestion in that discussion was to use something like:

plots[display](
  [plot( x^2, x=0..1, filled=true, color=white ),
   plot(   x, x=0..1, filled=true, color=blue  )]
);

I thought this was generalized, but I do not see that post. So, I created the following procedure:

PlotFill := proc( f, g )
  local C, opts;
  C,opts := selectremove( has, [_passed[3..-1]], {color,colour} );
  if C=NULL then C:=[color=red] end if;
  plots[display]( [plot(min(f,g), opts[], filled=true, color=white),
                   plot(max(f,g), opts[], filled=true, C[])] )
end proc:

To use this, try

PlotFill( x,x^2, x=0..2, color=blue );

I know this works only when the functions are positive. When they are negative, the roles of max and min are reversed. I do not have time right now to add this feature. I'm sure someone else will do it before I get back to the computer later tonight. 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/

As a general rule of thumb, numerically, all matrices are nonsingular. To accurately determine is a matrix is singular, numerically, requires some special considerations. In general, this is difficult to convey to students in an introductory course in linear algebra. Consider the following two matrices.

with( LinearAlgebra ):
M := M := Matrix([[1e-7,2e-8],[1e-14,2e-15]]);
M2 := Matrix([[1e3,2e-8],[1e-14,2e-15]])

By default, Maple uses 10 significant digits for its floating-point calculations. Here are the rank, rref (not shown), and condition numbers for these two matrices:

Digits:=10;
                                     10
Rank( M );
                                      1
ReducedRowEchelonForm( M );

Rank( M2 );
                                      1
ReducedRowEchelonForm( M2 );

ConditionNumber( M ), ConditionNumber( M2 );
                                   23                17
                     5.878072505 10  , 5.000000000 10  

And now the same results with Digits=15:

Digits := 15;
                                     15
Rank( M );
                                      1
ReducedRowEchelonForm( M );

Rank( M2 );
                                      1
ReducedRowEchelonForm( M2 );

ConditionNumber( M ), ConditionNumber( M2 );
                                                       17
                   Float(infinity), 5.00000000060000 10  

The results with Digits=14 agree with the default. The condition number is a common measure of the closeness of a matrix to being singular. In general, we want to keep the condition number small but there is no firm limit to what small means. This is a challenging topic to understand let alone to teach. 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/

It is not reasonable to ask a function to accept multiple inputs. That would change domain, and hence the function would be different.

What does make sense is to apply the function to a set (or list) of arguments. This is what the map function does.

Because of the way Maple treat parentheses, I do not know how to do this with your function of 2 variables. I would convert this to a function a single variable (a 2-element list), as follows:

f := X->X[1]^2+X[2]^2;

Then, I would use map as follows:

map( f, { [0,0], [1,1] } );
                                   {0, 2}
map( f, { [1,0], [0,1] } );
                                     {1}
map( f, [ [1,0], [0,1] ] );
                                   [1, 1]

 

Notice that Maple's sets really are sets, they do not allow for multiple copies of the same element. If you want to see the output for each argument - in the order corresponding to the order of the inputs - use a list (as in my third example).

Read the online help for map (and map2) and see the examples for more information.

Hoping this is relevant for your purposes,

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/

Change sum to add and all will be fine, I think.

M:=Matrix(3,(i,j)->f(i,j));
                          
g:=proc(M,i,j);
op(2,M[i,j])*op(1,M[i,j])*(2-op(1,M[i,j]))
end proc:
 

eq:=op(2,M[1,2])=g(M,2,2)+g(M,3,2);
                                   2 = -6

eq1:=op(2,M[1,2])=add( g(M,i,2) , i=2..3 );
                                   2 = -6

The output is garbage because I do not know your f. But, I get the same output both ways.

Your misuse of sum is very common (is it on the list of Top Ten Maple Errors?) but the online help is pretty clear. sum is for sums (finite or infinite) that you expect to have a closed form expression to replace the sum (e.g., sum( k^2, k=1..n )). add simply adds up a collection of numbers (must be finite and the endpoints must be definite numbers).

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 approach I typically take to creating a table of function values is to build a Matrix.

Here's the way I think about this task, and then a procedure to automate the process:

restart;
f := x->x^2;
a,b,n := 0,2,5;             # endpoints of interval and number of points
dx := (b-a)/n;              # stepsize between sampled points
X := [seq(a+dx*i,i=0..n)];  # list of points to use to create the table

Now that the function and sample points are defined, it is time to construct the table. First, create the column headers

hdr := < < `x`          | `f(x)` >,
         < `----------` | `----------` > >;

Now, the body of the table can be constructed as a collection of rows, one for each sample point:

body := 

And, finally, combine the header and body to form the complete table:

< hdr, body >;

Of course, if you want the output as floating point numbers you can do something like the following:

body2 := 

Here is a procedure to automate this process:

T2 := proc( f::procedure, X::list )
  local body, hdr;
  hdr := < < `x`          | `f(x)` >,
         < `----------` | `----------` > >;
  body := 

And the previous two examples using the procedure:

T2( f, X );
T2( f, evalf(X) );

Note that if the Matrix is too large, Maple will not display the full matrix unless you do something like:

interface( rtablesize=50 );

where 50 is the (row or column) dimension of the largest Matrix that Maple will display.

This method generalizes to more columns in a natural way. I hope this will be 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.edu/~meade/

Your question is not completely clear to me, but I believe you are asking how to convert the output of solve (a set of equations) to an ordered pair. Your procedure currently looks like:

beta_k := proc(P)
 ...
 sol := fsolve({Eq17, Eq19}, {k = 0 .. .5, beta_star = Pi .. (3/2)*Pi}):
end proc;

with typical output

                 {beta_star = 4.491245861, k = 0.1974116613}
                 {beta_star = 4.491447306, k = 0.1883934427}
                 {beta_star = 4.491614448, k = 0.1805073158}
                 {k = 0.1735347756, beta_star = 4.491755359}

If you want to see the output as an ordered pair [k,beta_star] you need to add one line immediately before the end proc. Here is how that part of the code should look:

beta_k := proc(P)
 ...
 sol := fsolve({Eq17, Eq19}, {k = 0 .. .5, beta_star = Pi .. (3/2)*Pi}):
 eval( [k,beta_star], sol );
end proc;

With this change, the above output appears as:

                         [0.1974116613, 4.491245861]
                         [0.1883934427, 4.491447306]
                         [0.1805073158, 4.491614448]
                         [0.1735347756, 4.491755359]

Of course, if you want the values in the opposite order, the modification is obvious. What might not be quite so obvious is that you can use eval to compute complex expressions (e.g., k^2+beta_star^2) that can be expressed in terms of the values in the set (or list) in the second argument. 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/

There are many options to the Matrix command (see ?Matrix). It takes a lot of reading to see everything, but what you want is the shape= option (see ?symmetric).

A := Matrix(3,3,symbol=a, shape=symmetric );

---------------------------------------------------------------------
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/

Yep. That's what I had in mind when I started to work on this.

But, as I got into it, I lost my focus.

Acer's explanation is much better than what I would have written. My point about the package was that thereis no need to load the entire package globally.
All I did was a shortcut to using the long form of the name.You are correct that the real performance hit is the creation of the procedures to do the
sampling. This is not an area where I am very knowledgeable.

As I wrote in my post, I was initially pleased with my improvements, but these disappeared with additional testing.

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 would have done this with seq, as follows:

q:={x = 0, y = -1}, {x = 3/4, y = -7/16};
                                       /    3      -7\ 
                     {y = -1, x = 0}, { x = -, y = -- }
                                       \    4      16/ 
seq( eval([x,y],p), p=[q] );
                                       [3  -7]
                              [0, -1], [-, --]
                                       [4  16]

---------------------------------------------------------------------
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/