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

In almost any 2-D plot command you can include the optional argument numpoints=###. This is not really changing the default as it applies only to the command in which it appears. To set default settings for plots, you would use the setoptions command in the plots package. The basic usage is:
with( plots ):
setoptions( numpoints=500 );
I know this will be used in any explicit plot command executed within a worksheet or document. I have not tested this with the plot builder. If it doesn't work, I'd call this a bug (certainly not a feature). 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/
Gentlemen, Read your respective postings carefully. The original posting asks about evalhf( int( ... ) ); and Jacques' response references evalhf( Int( ... ) );. This difference has been the subject of several recent postings on MaplePrimes. To summarize, the inert Int prevents Maple from attempting any symbolic simplifications; it goes directly to numerical evaluation. Using int, Maple first attempts to find a symbolic expression for the definite integral. Numerical integration is used only if Maple cannot find a symbolic form for the definite integral. A very subtle change in the Maple usage but a significant change in what Maple will 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/
Gentlemen, Read your respective postings carefully. The original posting asks about evalhf( int( ... ) ); and Jacques' response references evalhf( Int( ... ) );. This difference has been the subject of several recent postings on MaplePrimes. To summarize, the inert Int prevents Maple from attempting any symbolic simplifications; it goes directly to numerical evaluation. Using int, Maple first attempts to find a symbolic expression for the definite integral. Numerical integration is used only if Maple cannot find a symbolic form for the definite integral. A very subtle change in the Maple usage but a significant change in what Maple will 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/
In my previous response I meant to compare your original solution with the definition of a function. Here is an example that is comparable to what you were trying to do:
f(x) := x^2;
f(2);            # Maple knows nothing about f(2), only f(x)
f(y+z);
f(x);            # all that has been defined is f(x)
Note that the definition of f(x) is made with no meaning assigned to x. This only defines the literal name f(x). Here is one way to define a function in Maple:
F(x) := x -> x^2;
F(2);
F(y+z);
F(x);
Sometimes it is necessary to use the first idea to assign specific values to a function. Here is one implementation of the Fibonacci numbers. Note how a function is used in general, and then the first two values are assigned. Without these initial values, the definition would enter an infinite loop. As it is, don't try to compute a general Fib(n).
Fib := n -> Fib(n-1) + Fib(n-2);
Fib(0) := 1;     # define initial values
Fib(1) := 1;     # define initial values
Fib(5);
There are much better implementations of the Fibonacci numbers; see the online help for remember (?remember). 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/
In my previous response I meant to compare your original solution with the definition of a function. Here is an example that is comparable to what you were trying to do:
f(x) := x^2;
f(2);            # Maple knows nothing about f(2), only f(x)
f(y+z);
f(x);            # all that has been defined is f(x)
Note that the definition of f(x) is made with no meaning assigned to x. This only defines the literal name f(x). Here is one way to define a function in Maple:
F(x) := x -> x^2;
F(2);
F(y+z);
F(x);
Sometimes it is necessary to use the first idea to assign specific values to a function. Here is one implementation of the Fibonacci numbers. Note how a function is used in general, and then the first two values are assigned. Without these initial values, the definition would enter an infinite loop. As it is, don't try to compute a general Fib(n).
Fib := n -> Fib(n-1) + Fib(n-2);
Fib(0) := 1;     # define initial values
Fib(1) := 1;     # define initial values
Fib(5);
There are much better implementations of the Fibonacci numbers; see the online help for remember (?remember). 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 see. The basic problem is that your definition of u[i+1][k] is made when i and k have no meaning. Then, you try to use this definition when i and k have specific values. There are several ways to go about fixing this. The simplest is, I believe, like I suggested in my original post:
> restart;
> with(LinearAlgebra):
> interface( rtablesize=20 ):          # needed to be able to display the results as a matrix in the worksheet (default is 10)
> A:=LinearAlgebra:-RandomMatrix(6,5);
> M:=Transpose(A).A;
> u[0]:=LinearAlgebra:-RandomVector[column](5);
> b:=LinearAlgebra:-RandomVector[column](6);
> p:=LinearAlgebra:-RandomVector[column](5);
> w:=1.3;
> e:=0.2;
> for k from 1 to 5 do
> for i from 0 to 10 do                                 # note that starting value is 0, not 1
>   q1 := (add(M[k,j]*u[i][j]  ,j=1..5));               # change sum to add
>   q2 := e*b[k];
>   q3 := (add(A[k,j]*p[j]     ,j=1..5));               # change sum to add
>   q4 := (add(M[k,j]*u[i+1][j],j=1..k-1));             # change sum to add
>   q5 := (add(M[k,j]*u[i][j]  ,j=1..k-1));             # change sum to add
>   u[i+1][k]:=u[i][k]-w*(1/M[k,k])*(q1+q2+q3-q4+q5);   # compute new value of u
> end do;
> end do;
> Matrix( 5,11, (k,i)->u[i-1][k] );                     # display all computed values
Notes:
  1. It was necessary to change the starting value for i from 1 to 0.
  2. Also, you should have been using add instead of sum. (sum is used only for sums that can be simplified by a summation formula; use add to add a set of specific numbers.)
  3. You might be able to consolidate the loops into a single Matrix command, but you would have to be careful that the matrix is computed columnwise so that all necessary values are available.
  4. A more viable alternative would be to implement the algorithm in terms of matrix-vector operations. Each of the adds is really just a dot product and the inner loop could be avoided by seeing the calculation in terms of a matrix-vector product.
Please check that I have not changed your formula with the copying and pasting that I did. I hope this has been instrsuctive. 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 see. The basic problem is that your definition of u[i+1][k] is made when i and k have no meaning. Then, you try to use this definition when i and k have specific values. There are several ways to go about fixing this. The simplest is, I believe, like I suggested in my original post:
> restart;
> with(LinearAlgebra):
> interface( rtablesize=20 ):          # needed to be able to display the results as a matrix in the worksheet (default is 10)
> A:=LinearAlgebra:-RandomMatrix(6,5);
> M:=Transpose(A).A;
> u[0]:=LinearAlgebra:-RandomVector[column](5);
> b:=LinearAlgebra:-RandomVector[column](6);
> p:=LinearAlgebra:-RandomVector[column](5);
> w:=1.3;
> e:=0.2;
> for k from 1 to 5 do
> for i from 0 to 10 do                                 # note that starting value is 0, not 1
>   q1 := (add(M[k,j]*u[i][j]  ,j=1..5));               # change sum to add
>   q2 := e*b[k];
>   q3 := (add(A[k,j]*p[j]     ,j=1..5));               # change sum to add
>   q4 := (add(M[k,j]*u[i+1][j],j=1..k-1));             # change sum to add
>   q5 := (add(M[k,j]*u[i][j]  ,j=1..k-1));             # change sum to add
>   u[i+1][k]:=u[i][k]-w*(1/M[k,k])*(q1+q2+q3-q4+q5);   # compute new value of u
> end do;
> end do;
> Matrix( 5,11, (k,i)->u[i-1][k] );                     # display all computed values
Notes:
  1. It was necessary to change the starting value for i from 1 to 0.
  2. Also, you should have been using add instead of sum. (sum is used only for sums that can be simplified by a summation formula; use add to add a set of specific numbers.)
  3. You might be able to consolidate the loops into a single Matrix command, but you would have to be careful that the matrix is computed columnwise so that all necessary values are available.
  4. A more viable alternative would be to implement the algorithm in terms of matrix-vector operations. Each of the adds is really just a dot product and the inner loop could be avoided by seeing the calculation in terms of a matrix-vector product.
Please check that I have not changed your formula with the copying and pasting that I did. I hope this has been instrsuctive. 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/
You just type the two character "less than or equal to".
F := x -> piecewise( x<=0, 0, 1 );
F(-0.1); # returns 0
F(0);    # returns 0
F(0.1);  # returns 1
Note that the conditionals in the piecewise command can be as simple or complex as you want.
G := (x,y) -> piecewise( x>0 and x=y, 0, 1 );
G(0,1);  # returns 1
G(0,0);  # returns 1
G(1,1);  # returns 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/
You just type the two character "less than or equal to".
F := x -> piecewise( x<=0, 0, 1 );
F(-0.1); # returns 0
F(0);    # returns 0
F(0.1);  # returns 1
Note that the conditionals in the piecewise command can be as simple or complex as you want.
G := (x,y) -> piecewise( x>0 and x=y, 0, 1 );
G(0,1);  # returns 1
G(0,0);  # returns 1
G(1,1);  # returns 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/
I would use eval, as follows:
eval( MyExpression, c=cos(theta) );
You could use subs, with the syntax:
subs( c=cos(theta), MyExpression );
There are some minor technical differences, but I find the eval fits my usage best. A third option is to create a function, with unapply as follows:
MyFunction := unapply( MyExpression, c );
MyFunction( cos(theta) );
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.edu/~meade/
I would use eval, as follows:
eval( MyExpression, c=cos(theta) );
You could use subs, with the syntax:
subs( c=cos(theta), MyExpression );
There are some minor technical differences, but I find the eval fits my usage best. A third option is to create a function, with unapply as follows:
MyFunction := unapply( MyExpression, c );
MyFunction( cos(theta) );
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.edu/~meade/
NOTE: I assume the previous post had a less than sign that was eaten by the interpreter. piecewise is a Maple command; to see the online documentation, try:
?piecewise
Here are some examples that use piecewise:
> V := piecewise( x<0, infinity, 1 );    # your example?
> P := piecewise( x=0, 1, 0 );           # isolated point
> S := piecewise( x<0, -1, x=0, 0, 1 );  # multiple branches
> A := piecewise( x<0, -x, x );          # absolute value
Note that multiple branches are possible. There are some evaluation rules that can cause problems in some situations; these are explained in the online help. Why are you opposed to using if? There might be some cases where piecewise is not sufficient. You might be able to use the `if` operator (with the quotes). This is detailed in the online help for if (?if). 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/
NOTE: I assume the previous post had a less than sign that was eaten by the interpreter. piecewise is a Maple command; to see the online documentation, try:
?piecewise
Here are some examples that use piecewise:
> V := piecewise( x<0, infinity, 1 );    # your example?
> P := piecewise( x=0, 1, 0 );           # isolated point
> S := piecewise( x<0, -1, x=0, 0, 1 );  # multiple branches
> A := piecewise( x<0, -x, x );          # absolute value
Note that multiple branches are possible. There are some evaluation rules that can cause problems in some situations; these are explained in the online help. Why are you opposed to using if? There might be some cases where piecewise is not sufficient. You might be able to use the `if` operator (with the quotes). This is detailed in the online help for if (?if). 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/
To get started with Runge-Kutta methods, you can look at http://en.wikipedia.org/wiki/Runge-Kutta. You will also find this in many textbooks and other references. You also ask about higher-order ODEs. Most numerical methods are designed only for first-order ODEs. Higher-order ODEs are handled by converting them to a SYSTEM of first-order ODEs. The extra variables are the higher derivatives. For example, if y''' = f(x,y,y',y''), let u = y' and v = y'' and solve the system y' = u u' = v v' = f(x,y,u,v) I hope this get you started. 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/
To get started with Runge-Kutta methods, you can look at http://en.wikipedia.org/wiki/Runge-Kutta. You will also find this in many textbooks and other references. You also ask about higher-order ODEs. Most numerical methods are designed only for first-order ODEs. Higher-order ODEs are handled by converting them to a SYSTEM of first-order ODEs. The extra variables are the higher derivatives. For example, if y''' = f(x,y,y',y''), let u = y' and v = y'' and solve the system y' = u u' = v v' = f(x,y,u,v) I hope this get you started. 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 68 69 70 71 72 73 74 Page 70 of 76