In my previous blog entry, I neglected to explain how I discovered this problem with iscont. Recall:

f := 2^x+x^4:

In my work on my Maplets for Calculus project, I used a general command similar to the following to create a plot within a maplet:

plot( sqrt(f), x=1..10, discont=true );
Error, (in plot) invalid arguments

What could be invalid here? None of the following cause any trouble:

plot( f, x=1..10, discont=true );
plot( f, x=1..10 );
plot( sqrt(f), x=1..10 );

After more time than I want to admit, it dawned on me that the discont=true option is causing the trouble. As I dug deeper into this, I discovered the iscont oddity reported in my previous post. I overcame this situation with the use of try ... catch as follows:

try
  __p4:=plot( sqrt(f), x=1..10, discont=true, color=red );
catch:
  __p4:=plot( sqrt(f), x=1..10, color=red );
end try;

Why does this example cause so many problems? We now look at how discont handles this problem:

discont(f,x);
                                     {}
discont(sqrt(f),x);
               /            /       /1   1  \        (1/2)\  
               |  4 LambertW|_NN18, |- - - I| ln(2) 2     |  
              <             \       \8   8  /             /  
               |- -----------------------------------------, 
               \                    ln(2)                    

                            /       /  1   1  \        (1/2)\  
                  4 LambertW|_NN17, |- - - - I| ln(2) 2     |  
                            \       \  8   8  /             /  
                - -------------------------------------------, 
                                     ln(2)                     

                            /       /1   1  \        (1/2)\  
                  4 LambertW|_NN19, |- + - I| ln(2) 2     |  
                            \       \8   8  /             /  
                - -----------------------------------------, 
                                    ln(2)                    

                            /       /  1   1  \        (1/2)\\ 
                  4 LambertW|_NN20, |- - + - I| ln(2) 2     || 
                            \       \  8   8  /             / >
                - -------------------------------------------| 
                                     ln(2)                   / 
evalf(%);
      {-5.770780164 LambertW(_NN18, 0.1225322679 - 0.1225322679 I), 

        -5.770780164 LambertW(_NN19, 0.1225322679 + 0.1225322679 I), 

        -5.770780164 LambertW(_NN20, -0.1225322679 + 0.1225322679 I), 

        -5.770780164 LambertW(_NN17, -0.1225322679 - 0.1225322679 I)}

The uses of discont=true and iscont (with a real interval) trapped me into thinking only about the real-valued cases. It never dawned on me (until seeing the output from discont) that this result was being complicated by the more general complex-valued case. Knowing what I now know, is there a clean way to overcome this problem (without my try ... catch kludge)? I do see that fdiscont can be used in this situation, but I am not sure about using this in general (and it's not an option within the plot command).

fdiscont(f,x=1..infinity);
                                     []
fdiscont(sqrt(f),x=1..infinity);
                                     []

Hopefully, some of this insight will help someone to either i) understand the problem in a way that enables them to provide a better solution to this problem or ii) suggest an efficient way to handle such situations. 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 looks as though you are trying to say that h is a piecewise defined function. If so, then you will need to use the piecewise command. For example,

h := x->piecewise( x>0 and x<1/2, 2*x, 0 );

There are other issues with the little snippet of code that I can see in your original post. First, sum should be replaced with add. (While it works with sum, sum is for use with formal sums that you expect will have a closed form or other simplification.) You can now plot the FUNCTION f as follows:

plot( f, 0..1 );

or the EXPRESSION f(x) as follows:

plot( f(x), x=0..1 );

Since this function is discontinuous, you may want to include the option discont=true in either of the above. If you want to include the base function, h, in the plot you might want to do something like

plot( [f,h], 0..1, discont=true, style=[line,point] );

If you happen to want to look at the value of f(x), you will see the sum of 21 piecewise functions. We can force Maple to simplify this to a single piecewise-defined function with the following:

f(x);
convert( f(x), piecewise );

Guessing at your ultimate interests, this form might be very enlightening or instructive. 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 looks as though you are trying to say that h is a piecewise defined function. If so, then you will need to use the piecewise command. For example,

h := x->piecewise( x>0 and x<1/2, 2*x, 0 );

There are other issues with the little snippet of code that I can see in your original post. First, sum should be replaced with add. (While it works with sum, sum is for use with formal sums that you expect will have a closed form or other simplification.) You can now plot the FUNCTION f as follows:

plot( f, 0..1 );

or the EXPRESSION f(x) as follows:

plot( f(x), x=0..1 );

Since this function is discontinuous, you may want to include the option discont=true in either of the above. If you want to include the base function, h, in the plot you might want to do something like

plot( [f,h], 0..1, discont=true, style=[line,point] );

If you happen to want to look at the value of f(x), you will see the sum of 21 piecewise functions. We can force Maple to simplify this to a single piecewise-defined function with the following:

f(x);
convert( f(x), piecewise );

Guessing at your ultimate interests, this form might be very enlightening or instructive. 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 looks to be another instance of an < in the text of a post. To be able to see < in the body of a post, you will need to enter it as &lt; (with the semi-colon) or use an HTML tag to force the desired appearance.

---------------------------------------------------------------------
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 looks to be another instance of an < in the text of a post. To be able to see < in the body of a post, you will need to enter it as &lt; (with the semi-colon) or use an HTML tag to force the desired appearance.

---------------------------------------------------------------------
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 you have already answered your original question, I want to respond to the issue raised by the first suggestion. If you look carefully, you will see that plot was plotting [sin,D(sin)] these are the sine and cosine FUNCTIONS. If you changed this to [sin(x),D(sin)(x)], I bet you'd get the same response as for a polynomial. You would need to either define the polynomial as a function with -> or create a function from an expression with unapply.

F := x -> x^4-3*x^2;
g := x^4-3*x^2;
G := unapply( g, x );

Then, in the plot you would use [F,D(F)] and [G,D(G)], respectively. Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183

While you have already answered your original question, I want to respond to the issue raised by the first suggestion. If you look carefully, you will see that plot was plotting [sin,D(sin)] these are the sine and cosine FUNCTIONS. If you changed this to [sin(x),D(sin)(x)], I bet you'd get the same response as for a polynomial. You would need to either define the polynomial as a function with -> or create a function from an expression with unapply.

F := x -> x^4-3*x^2;
g := x^4-3*x^2;
G := unapply( g, x );

Then, in the plot you would use [F,D(F)] and [G,D(G)], respectively. Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183

Re-read the following from my previous post: You start with an initial condition (u[i,0], i=0..nxpoints). This allows you to compute the solution at the first timestep (u[i,1], i=1..nxpoints-1) and thenyou could use the boundary conditions to solve for u[0,1] and u[nxpoints,1]. Now, repeat for later timesteps. You cannot apply the boundary conditions before you do any computations. They have to be applied one timestep at a time. Solve in the interior, then compute the values at the boundary that satisfy the discrete boundary condition. Then, move on to the next time step. One last comment. Your initial condition appears to be u(x,0)=sin(x). So, u[x](x,0)=cos(x) and so u[x](0,0)=cos(0)=1 and u[x](1,0)=cos(1). So, your initial condition does not satisfy your stated boundary conditions. As such, your problem is inconsistent. I would encourage you to take some time to learn a little more about the PDE. This will make the coding much easier because you will see how each step meshes with the mathematics and will be able to look at your computed solution and know if it is reasonable. I do not write this to be critical or condescending. It's the advise I would give anyone. The problems you are trying to solve are not trivial and trying to code a solution without knowing what is reasonable and what problems have to be avoided is asking for trouble. 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/

Re-read the following from my previous post: You start with an initial condition (u[i,0], i=0..nxpoints). This allows you to compute the solution at the first timestep (u[i,1], i=1..nxpoints-1) and thenyou could use the boundary conditions to solve for u[0,1] and u[nxpoints,1]. Now, repeat for later timesteps. You cannot apply the boundary conditions before you do any computations. They have to be applied one timestep at a time. Solve in the interior, then compute the values at the boundary that satisfy the discrete boundary condition. Then, move on to the next time step. One last comment. Your initial condition appears to be u(x,0)=sin(x). So, u[x](x,0)=cos(x) and so u[x](0,0)=cos(0)=1 and u[x](1,0)=cos(1). So, your initial condition does not satisfy your stated boundary conditions. As such, your problem is inconsistent. I would encourage you to take some time to learn a little more about the PDE. This will make the coding much easier because you will see how each step meshes with the mathematics and will be able to look at your computed solution and know if it is reasonable. I do not write this to be critical or condescending. It's the advise I would give anyone. The problems you are trying to solve are not trivial and trying to code a solution without knowing what is reasonable and what problems have to be avoided is asking for trouble. 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 should have been more explicit, N=nxpoints, corresponding to x=1 (a boundary). The boundary condition should be u[x](0,t)=4 and u[x](1,t)=8 for all time t>0. It does not make any sense to apply these for N=0..nxpoints.

(u[N,j]-u[N-1,j])/h=4

This is a discrete approximatino to the x-derivative of u at x=1. In the limit, as h->0, the left-hand side converges to u[x](1,t). (With your revision, the right-hand side should now be 8.) These conditions are to be applied at each time step, i.e., for each j. Look at the equations that you use in the interior. These should be discrete approximations to the PDE (heat equation). Something like u[t](x,t)=u[xx](x,t) would be discretized as (u[i,j+1]-u[i,j])/k = (u(i-1,j)-u(i+1,j))/h^2 where k is the size of the timestep and h is the size of the spatial step. You start with an initial condition (u[i,0], i=0..nxpoints). This allows you to compute the solution at the first timestep (u[i,1], i=1..nxpoints-1) and thenyou could use the boundary conditions to solve for u[0,1] and u[nxpoints,1]. Now, repeat for later timesteps. Different approximations for any of these derivatives lead to different approximate solutions. If you approximate the time derivative in the opposite direction, (u[i,j]-u[i,j-1])/k, then you cannot solve directly for u[i,j]. Instead, you have to solve the full set of linear equations simultaneously on ech timestep. All of this can be found in almost any introductory textbook on numerical solutions to PDEs. 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 should have been more explicit, N=nxpoints, corresponding to x=1 (a boundary). The boundary condition should be u[x](0,t)=4 and u[x](1,t)=8 for all time t>0. It does not make any sense to apply these for N=0..nxpoints.

(u[N,j]-u[N-1,j])/h=4

This is a discrete approximatino to the x-derivative of u at x=1. In the limit, as h->0, the left-hand side converges to u[x](1,t). (With your revision, the right-hand side should now be 8.) These conditions are to be applied at each time step, i.e., for each j. Look at the equations that you use in the interior. These should be discrete approximations to the PDE (heat equation). Something like u[t](x,t)=u[xx](x,t) would be discretized as (u[i,j+1]-u[i,j])/k = (u(i-1,j)-u(i+1,j))/h^2 where k is the size of the timestep and h is the size of the spatial step. You start with an initial condition (u[i,0], i=0..nxpoints). This allows you to compute the solution at the first timestep (u[i,1], i=1..nxpoints-1) and thenyou could use the boundary conditions to solve for u[0,1] and u[nxpoints,1]. Now, repeat for later timesteps. Different approximations for any of these derivatives lead to different approximate solutions. If you approximate the time derivative in the opposite direction, (u[i,j]-u[i,j-1])/k, then you cannot solve directly for u[i,j]. Instead, you have to solve the full set of linear equations simultaneously on ech timestep. All of this can be found in almost any introductory textbook on numerical solutions to PDEs. 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 trick to fill in between two curves is to fill below the upper curve with the desired fill color and then to fill below the lower curve with the background color (white). Of course, this works only when you can clearly identify the upper and lower curves (unless you use max and min, ...). To complicate the matter, the plots have to be received by plots[display] so that they the lower curve is on top. Here's a simple example:

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

It seems that this is something that could be automated. 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 stand by the description in my original post. Here is what I get when I sent Jean-Marc's response to my original post to myself: -------------CUT HERE-------------- MaplePrimes Doug Meade thought you would like to see this page from MaplePrimes Message from Sender: This should be Jean-Marc's message, as posted to MaplePrimes. n/a Click here to read more on our site -------------CUT HERE-------------- Note that I did not receive Jean-Marc's message. I am using Firefox as my browser (and read e-mail with Thunderbird). Could this be a factor in the different experiences? Doug P.S. Jean-Marc, if you will send me your e-mail address, I will send the message to you (from within MaplePrimes) and you can see exactly what I receive. 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 O terms in a series solution are often useful, but ultimately you usually want to get rid of them (for plotting, etc.). The easiest way I know to do this is:

eval( ..., O=0 );  # capital oh = zero

You can also convert the series to a polynomial - even if the powers are not all positive integers:

convert( ..., polynom );

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 O terms in a series solution are often useful, but ultimately you usually want to get rid of them (for plotting, etc.). The easiest way I know to do this is:

eval( ..., O=0 );  # capital oh = zero

You can also convert the series to a polynomial - even if the powers are not all positive integers:

convert( ..., polynom );

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/