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

Now we can make some progress. The problems are all in the piecewise definition of FOC. As we will see, this expression needs to be simplified. Then, a plot is almost possible. After simplification, however, it becomes almost reasonable to find an explicit formula for Q_opt. This does lead to a satisfactory plot. Let's look at this step-by-step. Along the way, we will see that there are some troubles with undefined and conversions between piecewise and Heaviside. All of these can be avoided with not-too-complicated manipulations. Your original definition of FOC is rather complicated.
restart;
FOC := s*Q+1/4*s^2-PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])^2*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])+Q^2+PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])^2+2*Q*PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])-2*Q*PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])-Q^2*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])-s*PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])+pa*PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])+1/4*pa^2*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])-s*(1-PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise]))+pa*(1-PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise]))+pa*(-1/2*pa-Q)*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])-p:
The easiest way to simplify a piecewise-defined function is often to convert it to a Heaviside representation and then back to piecewise.
convert( FOC, Heaviside );
                        1  2                            2   1   2
     s Q - s - p + pa + - s  + Heaviside(pa + 2 Q - 2) Q  - - pa 
                        4                                   4    

        + Heaviside(pa + 2 Q - 2) Q pa - 2 Heaviside(pa + 2 Q - 2) Q

                                       1                           2
        - Heaviside(pa + 2 Q - 2) pa + - Heaviside(pa + 2 Q - 2) pa 
                                       4                            

        + Heaviside(pa + 2 Q - 2) - Q pa + undefined Dirac(pa + 2 Q - 2)

        + undefined Dirac(pa + 2 Q - 2) Q + undefined
convert( %, piecewise, Q );
           /        1     
  piecewise|Q = 1 - - pa, 
           \        2     

      /    1   \                1  2   1   2   /    1   \                 
    s |1 - - pa| - s - p + pa + - s  - - pa  - |1 - - pa| pa + undefined, 
      \    2   /                4      4       \    2   /                 

                       1  2   1   2                   \
    s Q - s - p + pa + - s  - - pa  - Q pa + undefined|
                       4      4                       /
In this case, it appears that the undefined are not being properly handled. While this result is not entirely useful, neither is it useless. We see that there the definition depends on the relation between Q and 1-pa/2. Let's look at the three cases separately:
FOClt := simplify( FOC ) assuming Q<1-pa/2;
                         1  2   1   2                    
                   s Q + - s  - - pa  - Q pa - s + pa - p
                         4      4                        
FOCeq := simplify( FOC ) assuming Q=1-pa/2;
           1  2    2                         2                         
     s Q + - s  - Q  + undefined Q - Q pa + Q  undefined + undefined pa
           4                                                           

            2                              
        + pa  undefined - s - p + undefined
FOCgt := simplify( FOC ) assuming Q>1-pa/2;
                            1  2        2              
                      s Q + - s  + 1 + Q  - 2 Q - s - p
                            4                          
This looks manageable. Putting this back together as FOC, we have:
FOC := piecewise( Q<1-pa/2, FOClt, Q=1-pa/2, undefined, Q>1-pa/2, FOCgt );
         /        1           1  2   1   2                              1     
piecewise|Q < 1 - - pa, s Q + - s  - - pa  - Q pa - s + pa - p, Q = 1 - - pa, 
         \        2           4      4                                  2     

                 1               1  2        2              \
  undefined, 1 - - pa < Q, s Q + - s  + 1 + Q  - 2 Q - s - p|
                 2               4                          /
Now, to define Q_opt in a form suitable for plotting. We begin with the definition of Q_opt given previously.
Q_opt := proc( s_val, p_val, pa_val )
  local q;
  if not is( {args}, set(numeric) ) then
    return 'procname'('args')
  end if;
  q := fsolve( eval( FOC, [p=p_val,s=s_val,pa=pa_val]), Q );
#  if is(q,numeric) then q else NULL end if;
end proc:
The line that has been commented out detects instances when fsolve does not return a numeric value. For example:
Q_opt( 0.5, 0, 0.6 );
                                0.7500000000
Q_opt( .5075205526, 0, 0.6 );
           /         /                                                  
     fsolve\piecewise\Q < 0.7000000000, -0.0924794474 Q + 0.0668737252, 

       Q = 0.7000000000, undefined, 0.7000000000 < Q, 

                                        2\   \
       -1.492479447 Q + 0.5568737254 + Q /, Q/
How did I come up with this example? When I first looked at this, I wanted to see how many times, and with what arguments, Q_opt was called. For this I uncommented the following command:
#debug( Q_opt );
I saw that Q_opt was being called many times, with each call be answered quickly. Because of the way Maple uses (by default) an adaptive selection of points to construct the plot, it was finding lots of these points where Q_opt was not returning a numeric value. This is why you were not seeing any response from your attempts to plot Q_opt(s,0,0.6). We can reduce the amount of work Maple does by lowering the value of numpoints (or turning off adaptive point selection with adaptive=false):
P1 := plot( Q_opt(s,0,0.6), s=0.5..0.59, axes=BOXED, numpoints=10, thickness=3 ):
P1;
When you look at this plot (in the uploaded worksheet), you will see that it is mostly linear, with a few gaps in regions where Q_opt does not return a numeric value. As noted earlier, the definition of FOC is not so complicated. Let's try to solve each piece of the definition for Q in the hope that we can construct an explicit definition for Q_opt.
Qlt := solve( FOC=0, Q ) assuming Q<1-pa/2;
                                 2     2             
                         -4 s + s  - pa  + 4 pa - 4 p
                       - ----------------------------
                                 4 (-pa + s)         
Qgt := solve( FOC=0, Q ) assuming Q>1-pa/2;
                     1          (1/2)    1          (1/2)
                   - - s + 1 + p     , - - s + 1 - p     
                     2                   2               
Note that the second part returned two values for Q. I will choose to use the first one; you might have a better reason to make a different choice. (In your example, p=0, so this choice is irrelevant. Here is my explicit definition of Q_opt:
myQ_opt := proc( s, p, pa )
  local Qlt, Qgt;
  if not is({args},set(numeric)) then
    return 'procname'('args');
  end if;
  Qlt := -(1/4)*(-4*s+s^2-pa^2+4*pa-4*p)/(-pa+s);
  Qgt := -(1/2)*s+1+sqrt(p);#, -(1/2)*s+1-sqrt(p)
  if   Qlt<1-pa/2 then Qlt
  elif Qgt>1-pa/2 then Qgt
  else undefined
  end if;
end proc:
Let's see how this responds to the troublesome value of s used earlier:
myQ_opt(.5075205526,0,0.6);
                                0.7462397237
Looks good! To conclude, here is the plot of the simplified Q_opt. It's not necessary to use numpoints here.
P2 := plot( myQ_opt(s,0,0.6), s=0.5..0.59, axes=BOXED, color=blue ):
P2;
When the two plots are overlayed, it's clear we have been successful.
plots[display](P1,P2);
To conclude, it's now possible to plot the explicitly defined Q_opt over a larger interval:
plot( myQ_opt(s,0,0.6), s=0..2, axes=BOXED, color=blue );
Everything is in the worksheet that I have uploaded: View 178_Q_opt.mw on MapleNet or Download 178_Q_opt.mw
View file details 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/
Now we can make some progress. The problems are all in the piecewise definition of FOC. As we will see, this expression needs to be simplified. Then, a plot is almost possible. After simplification, however, it becomes almost reasonable to find an explicit formula for Q_opt. This does lead to a satisfactory plot. Let's look at this step-by-step. Along the way, we will see that there are some troubles with undefined and conversions between piecewise and Heaviside. All of these can be avoided with not-too-complicated manipulations. Your original definition of FOC is rather complicated.
restart;
FOC := s*Q+1/4*s^2-PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])^2*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])+Q^2+PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])^2+2*Q*PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])-2*Q*PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])-Q^2*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])-s*PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])+pa*PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise])*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])+1/4*pa^2*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])-s*(1-PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise]))+pa*(1-PIECEWISE([1, 1 < Q+1/2*pa], [Q+1/2*pa, otherwise]))+pa*(-1/2*pa-Q)*PIECEWISE([1, Q < 1-1/2*pa], [undefined, Q = 1-1/2*pa], [0, 1-1/2*pa < Q])-p:
The easiest way to simplify a piecewise-defined function is often to convert it to a Heaviside representation and then back to piecewise.
convert( FOC, Heaviside );
                        1  2                            2   1   2
     s Q - s - p + pa + - s  + Heaviside(pa + 2 Q - 2) Q  - - pa 
                        4                                   4    

        + Heaviside(pa + 2 Q - 2) Q pa - 2 Heaviside(pa + 2 Q - 2) Q

                                       1                           2
        - Heaviside(pa + 2 Q - 2) pa + - Heaviside(pa + 2 Q - 2) pa 
                                       4                            

        + Heaviside(pa + 2 Q - 2) - Q pa + undefined Dirac(pa + 2 Q - 2)

        + undefined Dirac(pa + 2 Q - 2) Q + undefined
convert( %, piecewise, Q );
           /        1     
  piecewise|Q = 1 - - pa, 
           \        2     

      /    1   \                1  2   1   2   /    1   \                 
    s |1 - - pa| - s - p + pa + - s  - - pa  - |1 - - pa| pa + undefined, 
      \    2   /                4      4       \    2   /                 

                       1  2   1   2                   \
    s Q - s - p + pa + - s  - - pa  - Q pa + undefined|
                       4      4                       /
In this case, it appears that the undefined are not being properly handled. While this result is not entirely useful, neither is it useless. We see that there the definition depends on the relation between Q and 1-pa/2. Let's look at the three cases separately:
FOClt := simplify( FOC ) assuming Q<1-pa/2;
                         1  2   1   2                    
                   s Q + - s  - - pa  - Q pa - s + pa - p
                         4      4                        
FOCeq := simplify( FOC ) assuming Q=1-pa/2;
           1  2    2                         2                         
     s Q + - s  - Q  + undefined Q - Q pa + Q  undefined + undefined pa
           4                                                           

            2                              
        + pa  undefined - s - p + undefined
FOCgt := simplify( FOC ) assuming Q>1-pa/2;
                            1  2        2              
                      s Q + - s  + 1 + Q  - 2 Q - s - p
                            4                          
This looks manageable. Putting this back together as FOC, we have:
FOC := piecewise( Q<1-pa/2, FOClt, Q=1-pa/2, undefined, Q>1-pa/2, FOCgt );
         /        1           1  2   1   2                              1     
piecewise|Q < 1 - - pa, s Q + - s  - - pa  - Q pa - s + pa - p, Q = 1 - - pa, 
         \        2           4      4                                  2     

                 1               1  2        2              \
  undefined, 1 - - pa < Q, s Q + - s  + 1 + Q  - 2 Q - s - p|
                 2               4                          /
Now, to define Q_opt in a form suitable for plotting. We begin with the definition of Q_opt given previously.
Q_opt := proc( s_val, p_val, pa_val )
  local q;
  if not is( {args}, set(numeric) ) then
    return 'procname'('args')
  end if;
  q := fsolve( eval( FOC, [p=p_val,s=s_val,pa=pa_val]), Q );
#  if is(q,numeric) then q else NULL end if;
end proc:
The line that has been commented out detects instances when fsolve does not return a numeric value. For example:
Q_opt( 0.5, 0, 0.6 );
                                0.7500000000
Q_opt( .5075205526, 0, 0.6 );
           /         /                                                  
     fsolve\piecewise\Q < 0.7000000000, -0.0924794474 Q + 0.0668737252, 

       Q = 0.7000000000, undefined, 0.7000000000 < Q, 

                                        2\   \
       -1.492479447 Q + 0.5568737254 + Q /, Q/
How did I come up with this example? When I first looked at this, I wanted to see how many times, and with what arguments, Q_opt was called. For this I uncommented the following command:
#debug( Q_opt );
I saw that Q_opt was being called many times, with each call be answered quickly. Because of the way Maple uses (by default) an adaptive selection of points to construct the plot, it was finding lots of these points where Q_opt was not returning a numeric value. This is why you were not seeing any response from your attempts to plot Q_opt(s,0,0.6). We can reduce the amount of work Maple does by lowering the value of numpoints (or turning off adaptive point selection with adaptive=false):
P1 := plot( Q_opt(s,0,0.6), s=0.5..0.59, axes=BOXED, numpoints=10, thickness=3 ):
P1;
When you look at this plot (in the uploaded worksheet), you will see that it is mostly linear, with a few gaps in regions where Q_opt does not return a numeric value. As noted earlier, the definition of FOC is not so complicated. Let's try to solve each piece of the definition for Q in the hope that we can construct an explicit definition for Q_opt.
Qlt := solve( FOC=0, Q ) assuming Q<1-pa/2;
                                 2     2             
                         -4 s + s  - pa  + 4 pa - 4 p
                       - ----------------------------
                                 4 (-pa + s)         
Qgt := solve( FOC=0, Q ) assuming Q>1-pa/2;
                     1          (1/2)    1          (1/2)
                   - - s + 1 + p     , - - s + 1 - p     
                     2                   2               
Note that the second part returned two values for Q. I will choose to use the first one; you might have a better reason to make a different choice. (In your example, p=0, so this choice is irrelevant. Here is my explicit definition of Q_opt:
myQ_opt := proc( s, p, pa )
  local Qlt, Qgt;
  if not is({args},set(numeric)) then
    return 'procname'('args');
  end if;
  Qlt := -(1/4)*(-4*s+s^2-pa^2+4*pa-4*p)/(-pa+s);
  Qgt := -(1/2)*s+1+sqrt(p);#, -(1/2)*s+1-sqrt(p)
  if   Qlt<1-pa/2 then Qlt
  elif Qgt>1-pa/2 then Qgt
  else undefined
  end if;
end proc:
Let's see how this responds to the troublesome value of s used earlier:
myQ_opt(.5075205526,0,0.6);
                                0.7462397237
Looks good! To conclude, here is the plot of the simplified Q_opt. It's not necessary to use numpoints here.
P2 := plot( myQ_opt(s,0,0.6), s=0.5..0.59, axes=BOXED, color=blue ):
P2;
When the two plots are overlayed, it's clear we have been successful.
plots[display](P1,P2);
To conclude, it's now possible to plot the explicitly defined Q_opt over a larger interval:
plot( myQ_opt(s,0,0.6), s=0..2, axes=BOXED, color=blue );
Everything is in the worksheet that I have uploaded: View 178_Q_opt.mw on MapleNet or Download 178_Q_opt.mw
View file details 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/
Shawn, If, as in your example, the equations are linear then you can use the theory of linear algebra to get an answer to your question. Here's one way that this could be done:
with( LinearAlgebra ):
EQS := {
  x   + y   = a,
  x   - y   = b,
  2*x + 2*y = c
};
A,b := GenerateMatrix( EQS, [x,y] );
M := < A | b >;
GaussianElimination( M );
Now, if each row has a pivot entry, then the system is consistent for every right-hand side. If there is a row that does not have a pivot entry, then the RHS must be zero for the system to be consistent. Setting these elements of the RHS to zero gives the constraints that must be satisfied in order for the system to be consistent. Depending on the specific form of your equations, you can get fancier. Does this help you 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.edu/~meade/
It's possible that plot is trying to sample too many points to construct the plot. If so, this can be controlled via additional options to plot. To help me diagnose the specific problems, could you share the specific definition of FOC? (Upload the worksheet if that would be easier.) Do you see the plot (a straight line) when you use the specific definition for FOC that I assumed in 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.edu/~meade/
It's possible that plot is trying to sample too many points to construct the plot. If so, this can be controlled via additional options to plot. To help me diagnose the specific problems, could you share the specific definition of FOC? (Upload the worksheet if that would be easier.) Do you see the plot (a straight line) when you use the specific definition for FOC that I assumed in 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.edu/~meade/
It's morning, and nobody has followed up on the suggestion I made late last night. So, as promised, here it comes. To make this a more complete demonstration I made up an equation for FOC. The Q_opt procedure is very similar to the one originally posted. I prefer to use eval rather than subs, and all three values can be given at once (as a list or set). The real difference is the initial check that all three arguments are numeric. If an argument is not numeric, the return value is the procedure call unevaluated. Note that it is now much more straightforward to plot this function. No more unapply or rcurry. (Of course, technically, you are not plotting the procedure.)
> FOC := p+s=pa+Q;
                               p + s = pa + Q
> Q_opt := proc( s_val, p_val, pa_val )
>   if not is( {args}, set(numeric) ) then
>     return 'procname'('args')
>   end if;
>   fsolve( eval( FOC, [p=p_val,s=s_val,pa=pa_val]), Q );
> end proc:
> Q_opt(a,1,1);
                                Q_opt(a,1,1)
> plot( Q_opt(s,0,0.6), s=0.5..0.59, axes=BOXED );
While this plot is not very exciting for my definition of FOC, I hope this is helpful. Doug
---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/
It's morning, and nobody has followed up on the suggestion I made late last night. So, as promised, here it comes. To make this a more complete demonstration I made up an equation for FOC. The Q_opt procedure is very similar to the one originally posted. I prefer to use eval rather than subs, and all three values can be given at once (as a list or set). The real difference is the initial check that all three arguments are numeric. If an argument is not numeric, the return value is the procedure call unevaluated. Note that it is now much more straightforward to plot this function. No more unapply or rcurry. (Of course, technically, you are not plotting the procedure.)
> FOC := p+s=pa+Q;
                               p + s = pa + Q
> Q_opt := proc( s_val, p_val, pa_val )
>   if not is( {args}, set(numeric) ) then
>     return 'procname'('args')
>   end if;
>   fsolve( eval( FOC, [p=p_val,s=s_val,pa=pa_val]), Q );
> end proc:
> Q_opt(a,1,1);
                                Q_opt(a,1,1)
> plot( Q_opt(s,0,0.6), s=0.5..0.59, axes=BOXED );
While this plot is not very exciting for my definition of FOC, I hope this is helpful. Doug
---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/
I believe I have found the source of the problem. Maple 11 does try to combine the two RealRanges, but fails because signum(0.) returns as a float (0.) and not an integer (0). The code for RealRange_union looks for overlapping, or adjacent, intervals. It does this by looking at the signum of the difference of the endpoints of the ranges. In most cases, signum returns 1 or -1 (integers), but signum(0.) returns 0. (floating point). Here is the source for RealRange_union, with some comments added:
showstat( RealRange_union );

RealRange_union := proc()
local i, j, r, r1, r2, t1, t2, lo, hi, si;
   1   if nargs < 2 or not has([args],RealRange) then
   2     return args
       elif 2 < nargs then
for more than 2 ranges, work pairwise (recursion)
   3     for i to nargs do
   4       for j from i+1 to nargs do
   5         r := traperror(procname(args[i],args[j]));
   6         if r <> lasterror and nops([r]) = 1 then
   7           return procname(op(subsop(i = r,j = NULL,[args])))
             end if
           end do
         end do;
   8     return args
       else
extract endpoints of range 1
   9     r1 := traperror(RealRange_ends(args[1]));
  10     if r1 = lasterror then
  11       return args
         end if;
extract endpoints of range 2
  12     r2 := traperror(RealRange_ends(args[2]));
  13     if r2 = lasterror then
  14       return args
         end if;
  15     if type(r1[1],Open(algebraic)) then
  16       t1 := op(1,r1[1])
         else
  17       t1 := r1[1]
         end if;
check order of left endpoints of each range
  18     if type(r2[1],Open(algebraic)) then
  19       t2 := op(1,r2[1])
         else
  20       t2 := r2[1]
         end if;
  21     if has(t1,infinity) and has(t2,infinity) then
  22       si := 1
         else
  23       si := signum(t2-t1)
         end if;
  24     if not type(si,integer) then
  25       return args
         end if;
  26     if si = -1 or si = 0 and type(r1[1],Open(algebraic)) then
reverse order of the ranges
  27       t1 := r1;
  28       r1 := r2;
  29       r2 := t1
         end if;
At this time it is known that the left endpoint of range 1 <= right endpoint of range 2. Now, look for overlap by comparing the right endpoint of range 1 and the left endpoint of range 2.
  30     if type(r1[2],Open(algebraic)) then
  31       t1 := op(1,r1[2])
         else
  32       t1 := r1[2]
         end if;
  33     if type(r2[1],Open(algebraic)) then
  34       t2 := op(1,r2[1])
         else
  35       t2 := r2[1]
         end if;
  36     si := signum(t2-t1);
Reduction is NOT possible when: i) t2>t1 or ii) both ranges are open @ t1=t2 In all other cases, simply return args unchanged
  37     if not type(si,integer) then
  38       return args
         end if;
  39     if si = 1 or si = 0 and type(r1[2],Open(algebraic)) and type(r2[1],Open(algebraic)) then
  40       return args
         end if;
Finally, begin construction of the merged range. The left endpoint is easy (because of original ordering of the ranges).
  41     lo := r1[1];
The right endpoint takes more effort. (But, isn't there a better way to do this?)
  42     if type(r2[2],Open(algebraic)) then
  43       t2 := op(1,r2[2])
         else
  44       t2 := r2[2]
         end if;
  45     if has(t1,infinity) and has(t2,infinity) then
  46       si := 1
         else
  47       si := signum(t2-t1)
         end if;
  48     if not type(si,integer) then
  49       return args
         end if;
  50     if si = 1 or si = 0 and type(r1[2],Open(algebraic)) then
  51       hi := r2[2]
         else
  52       hi := r1[2]
         end if;
Final processing to handle infinity, then return the single combined range
  53     if has(lo,infinity) then
  54       lo := -infinity
         end if;
  55     if has(hi,infinity) then
  56       hi := infinity
         end if;
  57     RealRange(lo,hi)
       end if
end proc
This procedure appears to be needlessly complicated, and is certainly quite old. (This code is exactly the same as the code found in Maple 10.) I do not completely understand the rationale for using signum and not (simple) inequalities, but I trust there is (or was) a reason for this. Maybe something with the handling of infinity or symbolic arguments? My suggestion for fixing this routine is to replace the use of
if not type(si,integer) then
in lines 24, 37, and 48 with
if not type(convert(si,rational),integer) then
While I do not have the time to completely rewrite this procedure, I can see that it would benefit from
  • replacement of traperror/lasterror to try ... catch ... end try (lines 9-11 and 12-14)
  • use of multiple assignments to interchange two quantities (lines 27-29)
  • and typematch and patmatch
Of course, a full solution to this problem will involve a comprehensive look at RealRange_intersect and all of the related routines. I hope this diagnosis of the problem is helpful to someone who can implement a complete solution. 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/
Jacques, I beg to differ. The problem I reported about solving with inequalities and floats arose precisely because Maple 11 behaves different from Maple 10. So, this is a case where a change has been made and the change is NOT backward compatible. I am starting to dig into the code for solve and RealRanges. I believe there will be more to report, soon. 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/
Yes, I admit my original post did refer to both sets and RealRanges. When I suggested that the result of solve should be the union of two sets I was talking about two mathematical sets - which Maple represents, in this case, as RealRanges. Thank you for pointing out this ambiguity in my posting. Another question that I had, and did not include in my original post, is: how does Maple convert from the set
{ x < 2, 1 < x }
to the range
RealRange(Open(1),Open(2))
It is at this point that it would seem most likely to expect Maple to also look for combinations that can be simplified. As Robert points out, it appears this could be done with OrProp:
OrProp(RealRange(0,Open(1)), RealRange(1,2));
                               RealRange(0, 2)
Note that this slight revision of Robert's example clearly shows Maple has the ability to combine adjacent open and closed intervals as we would do by hand. Moreover, this is not affected by floating point numbers:
OrProp(RealRange(0.,Open(1.)), RealRange(1.,2.));
                              RealRange(0., 2.)
So, it looks as though my conjecture about a floating point gap was incorrect. (What a relief!) So, now it appears there was a change between Maple 10 and 11 that omits this simplification of the results from solve.
---------------------------------------------------------------------
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 Roman was posting his latest comments about functional programming and making Maple open source, I added a blog entry that is relevant. See Solving Inequalities with Piecewise and Floats. This shows a case where there has been a subtle change from Maple 10 to 11 that could have additional ramifications. It could be that the Maple 11 implementation is more mathematically correct but this has the consequence of making the results less useful. While there is not a formal mechanism for adding user-supplied code to Maple, this is not uncommon within Maple. The old Maple Share Library is an excellent example. Users could post code for others to access. Code that was particularly useful (measured in part by number of downloads) did make its way into the actual code - after thorough evaluation and revision by Maplesoft staff. I, for one, would like to see the resurrection of a 21st century version of the Maple Share Library. 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/
Randy, It's been a while since I have worked in MapleTA but my experience there - and in regular Maple - suggest that you need to use simplify (and sometimes normal or another command of this type) to force Maple to put all terms in the difference between the user's response and the "correct" answer. Let me illustrate:
F1 := sqrt(x):
F2 := x^(1/2):
F3 := x^(0.5):
S := {F1,F2,F3};
                              / 0.5   (1/2)\ 
                        S := { x   , x      }
                              \            / 
simplify( S );
                                  / (1/2)\ 
                                 { x      }
                                  \      / 
evalb( F1-F2=0 );
evalb( F1-F3=0 );
evalb( simplify( F1-F3 ) = 0 );
                                    true
                                    false
                                    true
The output from the definition of S indicates that Maple recognizes F1 and F2 as the same expression. After application of simplify, the set contains only a single element - so all three expressions are recognized as being equivalent. Note that while it is not an issue in this example, it is generally better to check if a difference is zero than to ask if two expressions are equal to one another. I hope this is helpful, Doug
---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/
Randy, It's been a while since I have worked in MapleTA but my experience there - and in regular Maple - suggest that you need to use simplify (and sometimes normal or another command of this type) to force Maple to put all terms in the difference between the user's response and the "correct" answer. Let me illustrate:
F1 := sqrt(x):
F2 := x^(1/2):
F3 := x^(0.5):
S := {F1,F2,F3};
                              / 0.5   (1/2)\ 
                        S := { x   , x      }
                              \            / 
simplify( S );
                                  / (1/2)\ 
                                 { x      }
                                  \      / 
evalb( F1-F2=0 );
evalb( F1-F3=0 );
evalb( simplify( F1-F3 ) = 0 );
                                    true
                                    false
                                    true
The output from the definition of S indicates that Maple recognizes F1 and F2 as the same expression. After application of simplify, the set contains only a single element - so all three expressions are recognized as being equivalent. Note that while it is not an issue in this example, it is generally better to check if a difference is zero than to ask if two expressions are equal to one another. I hope this is helpful, Doug
---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/
Jacques, Thanks for the suggestions for modernizing the code I proposed for convert/sectan. The code I wrote was, as I admitted, an adaption of the code Maple is currently using for convert/sincos. What does this say about the efficiency of the current code used in Maple. Is there a performance reason to prefer one style over another? Should we expect all of the Maple code to be (re-)written in the modern style? and updated as the language is updated? What are the benefits of the modern style, other than appearance and ease-of-HUMAN-understanding? 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/
mint is nice, for what it does, but it's does not help in the situations I described. I agree that it is nice to get mint's information about declared and unused variables and undeclared variables. But, I find mint of almost no use for SYNTAX checking. For managing large projects I, too, like the $include functionality. I wish it could be used within the user interface as well. Doug
---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/
First 64 65 66 67 68 69 70 Last Page 66 of 76