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

I think Axel was a little premature in his comments about "mw" files. I agree that it is frustrating to not see the full input. The real culprit here is 2d math input. If you use Maple notation for input, then there is no hiding of characters. Well, it would be even better if - like Word - there was a way to truly see every character in a worksheet. Just my two cents, 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/
If you look closely at the output you report, you will see that Maple did exactly as you asked. ALL occurrences of x(t) have been replaced by -0.5. This includes the occurrence of x(t) inside diff(x(t),t). You no longer have diff(x(t),t) so Maple does not see that the last substitution can be made. This problem can be avoided because the derivative dy/dt does not depend on x. So, you can omit the value for x(t) in the subs. Even though the subs will work, you will notice that the value for y(t) is inserted into dy/dt, thereby ruining this equation. Here is my suggested modification to your approach. Create a new expression that gives new names to the derivatives, say:
> dF2 := subs( diff(x(t),t)=DxDt, diff(y(t),t)=DyDt, dF1 );
Then you can substitute values using:
> subs( x(t)=-0.5, y(t)=(5/2)^(1/3), DxDt=-0.01, dF2 );
Note that you have mixed floating point and exact values. If you want all floating point numbers, change 5/2 into 5./2. If you want all exact values, change 0.5 into 1.2 and 0.01 into 1/100 - don't forget to do this in the definition of F1 also. This should get you closer to what you need.
---------------------------------------------------------------------
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/
If you look closely at the output you report, you will see that Maple did exactly as you asked. ALL occurrences of x(t) have been replaced by -0.5. This includes the occurrence of x(t) inside diff(x(t),t). You no longer have diff(x(t),t) so Maple does not see that the last substitution can be made. This problem can be avoided because the derivative dy/dt does not depend on x. So, you can omit the value for x(t) in the subs. Even though the subs will work, you will notice that the value for y(t) is inserted into dy/dt, thereby ruining this equation. Here is my suggested modification to your approach. Create a new expression that gives new names to the derivatives, say:
> dF2 := subs( diff(x(t),t)=DxDt, diff(y(t),t)=DyDt, dF1 );
Then you can substitute values using:
> subs( x(t)=-0.5, y(t)=(5/2)^(1/3), DxDt=-0.01, dF2 );
Note that you have mixed floating point and exact values. If you want all floating point numbers, change 5/2 into 5./2. If you want all exact values, change 0.5 into 1.2 and 0.01 into 1/100 - don't forget to do this in the definition of F1 also. This should get you closer to what you need.
---------------------------------------------------------------------
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/
Another useful tool is the showstat command. showstat( P ); produces a structured printout of the procedure P. This output includes line numbers. These are particularly useful when you want to focus a discussion on a specific portion of Maple code. For example,
> showstat( `int/xsincos`, 30..37 );

`int/xsincos` := proc(dx, ds, dc)
local a;
       ...
  30   if dx < 0 then
  31     if dc = 0 and ds = 1 then
  32       return -sin(_X)/(-dx-1)/(_X^(-dx-1))+1/(-dx-1)*`int/xsincos`(dx+1,0,1)
         end if;
  33     if ds = 0 and dc = 1 then
  34       return -cos(_X)/(-dx-1)/(_X^(-dx-1))-1/(-dx-1)*`int/xsincos`(dx+1,1,0)
         end if;
  35     if 0 <= dc and 0 <= ds then
  36       return `int/indef1`(expand(`trig/tfourier`(ds,dc,_X)*_X^dx,sin,cos))
         end if;
  37     return FAIL
       end if;
       ...
end proc

showstat( `trig/tfourier` );

`trig/tfourier` := proc(n, m)
local f, g, h;
   1   if m = 0 then
   2     `tfourier/sin`(n)/(2^n)
       elif n = 0 then
   3     `tfourier/cos`(m)/(2^m)
       else
   4     f := `tfourier/sin`(n);
   5     g := `tfourier/cos`(m);
   6     h := frontend(expand,[f*g]);
   7     `trig/mul`(h)/(2^(m+n))
       end if
end proc

showstat( `tfourier/sin` );

`tfourier/sin` := proc(n)
local k, x;
   1   x := _X;
   2   if n = 0 then
   3     1
       elif irem(n,2) = 0 then
   4     (-1)^modp(1/2*n,2)*((-1)^(1/2*n)*binomial(n,1/2*n)+2*sum('(-1)^k*binomial(n,k)*cos((n-2*k)*_X)',k = 0 .. 1/2*n-1))
       else
   5     2*(-1)^modp(1/2*n-1/2,2)*sum('(-1)^k*binomial(n,k)*sin((n-2*k)*_X)',k = 0 .. 1/2*n-1/2)
       end if
end proc

showstat( `tfourier/cos` );

`tfourier/cos` := proc(n)
local k, x;
   1   x := _X;
   2   if n = 0 then
   3     1
       elif irem(n,2) = 0 then
   4     2*sum('binomial(n,k)*cos((n-2*k)*_X)',k = 0 .. 1/2*n-1)+binomial(n,1/2*n)
       else
   5     2*sum('binomial(n,k)*cos((n-2*k)*_X)',k = 0 .. 1/2*n-1/2)
       end if
end proc
I hope some find this useful. Doug
---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/
Glenn, It sounds as though you might be trying to implement a shooting method. This has been done before in Maple. I wrote a general procedure (Shoot) to implement the shooting method. Others have done similar (often more specialized) implementations. A Google search for "Maple shooting method" turns up a nice list of options. My collection of shooting method tools can be found on my website at http://www.math.sc.edu/~meade/publ.html. This includes implementations for various versions of Maple. While my code is pretty complicated, you can see how I handle the Maple output from dsolve/numeric and the order in which results are reported. This experience was a major factor in the response I posted yesterday in this forum (under a different topic name). I am available to answer your questions if needed. 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 hope this thread (or, at least this comment) is sufficient to get this idea considered for a future release of Maple. I believe it is consistent with other parts of Maple and would be helpful. It is not sufficient to give dsolve a list instead of a set with the unknown functions. This changes the output but does not guarantee that the functions are in the order specified. A common way to re-order the output from dsolve is to use eval. For example, with yn as defined in the original post,
yn2 := eval( [y1(x),y2(x),y3(x)], yn );
The only drawback to this approach is that the list in yn2 contains only the right-hand sides of the corresponding solutions. (The output in yn is a list of equations.) But, in many cases the list of expressions is more useful than the list of equations. Note also that this approach can be used to create more interesting functions based on the solution. For example, the sum of the squares of the three components of the solution is
yn3 := eval( y1(x)^2+y2(x)^2+y3(x)^2, yn );
I hope some of 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 hope this thread (or, at least this comment) is sufficient to get this idea considered for a future release of Maple. I believe it is consistent with other parts of Maple and would be helpful. It is not sufficient to give dsolve a list instead of a set with the unknown functions. This changes the output but does not guarantee that the functions are in the order specified. A common way to re-order the output from dsolve is to use eval. For example, with yn as defined in the original post,
yn2 := eval( [y1(x),y2(x),y3(x)], yn );
The only drawback to this approach is that the list in yn2 contains only the right-hand sides of the corresponding solutions. (The output in yn is a list of equations.) But, in many cases the list of expressions is more useful than the list of equations. Note also that this approach can be used to create more interesting functions based on the solution. For example, the sum of the squares of the three components of the solution is
yn3 := eval( y1(x)^2+y2(x)^2+y3(x)^2, yn );
I hope some of 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/
Bikash, Your system is nonlinear, with quite a few unspecified parameters. It is highly unlikely that Maple will be able to find a general solution. My recommendation, particularly for a first step, is to look for a numerical solution. For this, you will need to specify numeric values of the parameters and initial conditions. Here is how that might look:
sys := { ... }; # as in your worksheet
param := { a=1, b=2, ... };
ic := { x1(0)=1, x2(0)=0, ... };
fns := { x1(t), x2(t), ... }
sol := dsolve( {sys,ic}, fns, type=numeric );
sol(1);  # approximate solution when t=1
plots[odeplot]( sol, [[t,x1(t)],[t,x2(t)]], 0..5 );
You should consult the plots[odeplot] help document for additional options regarding the types of plots that can be created. Note, in particular, that it is possible to plot in the phase-space, e.g., [x1(t),x2(t)], and to plot other functions of the solution, e.g., [t,x1(t)^2+x2(t)^2]. This is a pretty powerful command. The help document for dsolve,numeric will give more information about the types of numerical routines available for use. 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/
<\pre>
Bikash, Your system is nonlinear, with quite a few unspecified parameters. It is highly unlikely that Maple will be able to find a general solution. My recommendation, particularly for a first step, is to look for a numerical solution. For this, you will need to specify numeric values of the parameters and initial conditions. Here is how that might look:
sys := { ... }; # as in your worksheet
param := { a=1, b=2, ... };
ic := { x1(0)=1, x2(0)=0, ... };
fns := { x1(t), x2(t), ... }
sol := dsolve( {sys,ic}, fns, type=numeric );
sol(1);  # approximate solution when t=1
plots[odeplot]( sol, [[t,x1(t)],[t,x2(t)]], 0..5 );
You should consult the plots[odeplot] help document for additional options regarding the types of plots that can be created. Note, in particular, that it is possible to plot in the phase-space, e.g., [x1(t),x2(t)], and to plot other functions of the solution, e.g., [t,x1(t)^2+x2(t)^2]. This is a pretty powerful command. The help document for dsolve,numeric will give more information about the types of numerical routines available for use. 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/
<\pre>
How could a user find information about this use of the axes option of the plot command? I tried
?plot,options
<\pre>
and
?plot,axes
<\pre>

Via the Help Browser, I found the plot,structure document that gives information about the _AXIS[n] option.

---------------------------------------------------------------------
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/
<\pre>
How could a user find information about this use of the axes option of the plot command? I tried
?plot,options
<\pre>
and
?plot,axes
<\pre>

Via the Help Browser, I found the plot,structure document that gives information about the _AXIS[n] option.

---------------------------------------------------------------------
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/
<\pre>
The procedure in my previous post has one typo. It the P that appears in the indices command should be T. Here is the correct version:
DisplayTable := ( T::table ) ->
  plots[display]( [ seq(T[j], j=sort(map(op,[indices(T)]))) ], args[2..-1] );
Note that this procedure can be used even if the table is indexed by non-numeric indices. For example,
Q[red] := implicitplot( x^2+y^2=1, x=-1..1, y=-1..1, color=red ):
Q[blue] := implicitplot( x^2+y^2=1, x=-1..1, y=-1..1, color=blue ):
Q[green] := implicitplot( x^2+y^2=1, x=-1..1, y=-1..1, color=green ):
Q[pink] := implicitplot( x^2+y^2=1, x=-1..1, y=-1..1, color=pink ):
Q[cyan] := implicitplot( x^2+y^2=1, x=-1..1, y=-1..1, color=cyan ):
Q[magenta] := implicitplot( x^2+y^2=1, x=-1..1, y=-1..1, color=magenta ):
Q[turquoise] := implicitplot( x^2+y^2=1, x=-1..1, y=-1..1, color=turquoise ):

DisplayTable( Q, insequence=true );
If you want to use a different ordering that the one produced by sort, you will have to modify DisplayTable.
---------------------------------------------------------------------
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/
<\pre>
Dan, Here are some more ideas for you to consider. While Will's comments created the plots within the display command, my comments deal with the situation where you have created the plots separately, and want to combine them into a single display. First, create plots with two different naming conventions.
restart;
with( plots ):

for i from 1 to 10 do
  p||i := implicitplot( (x-i)^2+y^2=1, x=i-1..i+1, y=-1..1 );
  P[i] := implicitplot( x^2+(y-i)^2=1, x=-1..1, y=i-1..i+1 );
end do:
i; # note: i has a value (11) at the end of this loop
To work with the plots with concatenated names, we can use either $ or seq. With $ it is necessary to be aware that the loop variable in the for .. do .. end do loop retains its final value (11) at the end of the loop.
display( [p||i$i=1..10] );  # error because i has a value from before
display( [p||j$j=1..10] );  # works

display( [seq(p||i,i=1..10)] ); # uses a local variable named i
display( [seq(p||i,i=1..10)], scaling=constrained );
display( [seq(p||i,i=1..10)], scaling=constrained, insequence=true ); # animation
The last example shows one way to create animations. The table of plots can be used in the same way.
display( [ P[j] $ j=1..10 ] );
display( [ seq(P[j],j=1..10) ] );
As Will suggests in his post, it is possible to plot a table without explicit reference to the dimension of the table. It's possible to go even further than Will did with his example.
display( [ seq(P[i], i=sort(map(op,[indices(P)]))) ], scaling=constrained, insequence=true );

DisplayTable := ( T::table ) ->
  plots[display]( [ seq(T[j], j=sort(map(op,[indices(P)]))) ], args[2..-1] );

DisplayTable( P, scaling=constrained, insequence=true, thickness=2 );
I hope you found this useful, Doug
---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/
<\pre>
The issue, as I see it, is that the symbolic names of the elements of the matrix are not carried into the procedure. To overcome this it is necessary to assign values to the symbolic names in the original matrix. There are many ways to do this. All are essentially equivalent to the following:
f3 := proc( M::Matrix(2,2) )
  local i, j;
  global u;
  for i from 1 to RowDimension(M) do
    for j from 1 to ColumnDimension(M) do
      u[i,j]:=M[i,j]
    end do;
  end do;
  r;
end proc:
The problem with this implementation is that the elements, u[1,1], etc., retain their values outside this procedure. This means that r does as well. These problems can be overcome with some additional trickery:
f4 := proc( M::Matrix(2,2) )
  local i, j, _ReturnValue;
  global u;
  for i from 1 to RowDimension(M) do
    for j from 1 to ColumnDimension(M) do
      u[i,j]:=M[i,j]
    end do;
  end do;
  _ReturnValue := r;
  unassign(`u`);
  return _ReturnValue
end proc:
As I wrote the above, an even simpler approach came to mind:
f5 := proc( M::Matrix(2,2) )
  eval( r,
        [ seq( seq( u[i,j]=M[i,j],
                    j=1..ColumnDimension(M) ),
               i=1..RowDimension(M) ) ]
      );
end proc;
This implemenation avoids all intermediate names and all concerns about residual values of assignments. It could be simplified by creating a general procedure for creating the list of substitutions to be applied. I will leave it to others to decide if it is worth creating such a procedure, and to find an efficient implementation. I hope this has been helpful. Doug
----------------------------------------------------------------------
Prof. Douglas B. Meade     Phone:  (803) 777-6183
Department of Mathematics  URL:    http://www.math.sc.edu/~meade/
USC, Columbia, SC 29208    E-mail: mailto:meade@math.sc.edu
The issue, as I see it, is that the symbolic names of the elements of the matrix are not carried into the procedure. To overcome this it is necessary to assign values to the symbolic names in the original matrix. There are many ways to do this. All are essentially equivalent to the following:
f3 := proc( M::Matrix(2,2) )
  local i, j;
  global u;
  for i from 1 to RowDimension(M) do
    for j from 1 to ColumnDimension(M) do
      u[i,j]:=M[i,j]
    end do;
  end do;
  r;
end proc:
The problem with this implementation is that the elements, u[1,1], etc., retain their values outside this procedure. This means that r does as well. These problems can be overcome with some additional trickery:
f4 := proc( M::Matrix(2,2) )
  local i, j, _ReturnValue;
  global u;
  for i from 1 to RowDimension(M) do
    for j from 1 to ColumnDimension(M) do
      u[i,j]:=M[i,j]
    end do;
  end do;
  _ReturnValue := r;
  unassign(`u`);
  return _ReturnValue
end proc:
As I wrote the above, an even simpler approach came to mind:
f5 := proc( M::Matrix(2,2) )
  eval( r,
        [ seq( seq( u[i,j]=M[i,j],
                    j=1..ColumnDimension(M) ),
               i=1..RowDimension(M) ) ]
      );
end proc;
This implemenation avoids all intermediate names and all concerns about residual values of assignments. It could be simplified by creating a general procedure for creating the list of substitutions to be applied. I will leave it to others to decide if it is worth creating such a procedure, and to find an efficient implementation. I hope this has been helpful. Doug
----------------------------------------------------------------------
Prof. Douglas B. Meade     Phone:  (803) 777-6183
Department of Mathematics  URL:    http://www.math.sc.edu/~meade/
USC, Columbia, SC 29208    E-mail: mailto:meade@math.sc.edu
First 70 71 72 73 74 75 76 Page 72 of 76