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

To have Maple setup but not evaluate the integrals, you can use the intert version of int and limit: Int and Limit, respectively.
f := t;
                                      t
F := Int( f*exp(-s*t), t=a..infinity );
                           /infinity               
                          |                        
                          |          t exp(-s t) dt
                          |                        
                         /a                        
Limit( F, a=0, right );
                            /  /infinity               \
                            | |                        |
                      lim   | |          t exp(-s t) dt|
                    a -> 0+ | |                        |
                            \/a                        /
The eval command can be used to evaluate an inert expression. If you are not interested in the evaluation, then you might be able to use the inert integral to create something that looks more like what you are requesting. For example:
Int( f*exp(-s*t), t=0^`+`..infinity );
                           /infinity               
                          |                        
                          |          t exp(-s t) dt
                          |                        
                         / +                       
                          0                        
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 problem is not completely defined. How do you want the function extended beyond [0,1]? If you want it to be periodic with period 1, then take a look at the worksheet I just uploaded. View 178_FourierExpansion.mw on MapleNet or Download 178_FourierExpansion.mw
View file details Here are the commands, without output, from that worksheet:
restart;
f := x -> 2*x;
a := k -> 2*int( f(t)*cos(2*Pi*k*t), t=0..1 );
b := k -> 2*int( f(t)*sin(2*Pi*k*t), t=0..1 );
F := k -> a(0)/2 + add( a(j)*cos(2*Pi*j*x) + b(j)*sin(2*Pi*j*x), j=1..k );
F(0);
F(1);
F(2);
F(10);
plot( [f(x),F(0),F(1),F(2),F(10)], x=0..1, legend=["y=f(x)","k=0", "k=1", "k=2", "k=10"] );
If you want the function to be periodic with period 2, then is it an even or odd extension? These can be handled similarly, but will require slight adjustments to the computations of the coefficients. 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/
Your problem is not completely defined. How do you want the function extended beyond [0,1]? If you want it to be periodic with period 1, then take a look at the worksheet I just uploaded. View 178_FourierExpansion.mw on MapleNet or Download 178_FourierExpansion.mw
View file details Here are the commands, without output, from that worksheet:
restart;
f := x -> 2*x;
a := k -> 2*int( f(t)*cos(2*Pi*k*t), t=0..1 );
b := k -> 2*int( f(t)*sin(2*Pi*k*t), t=0..1 );
F := k -> a(0)/2 + add( a(j)*cos(2*Pi*j*x) + b(j)*sin(2*Pi*j*x), j=1..k );
F(0);
F(1);
F(2);
F(10);
plot( [f(x),F(0),F(1),F(2),F(10)], x=0..1, legend=["y=f(x)","k=0", "k=1", "k=2", "k=10"] );
If you want the function to be periodic with period 2, then is it an even or odd extension? These can be handled similarly, but will require slight adjustments to the computations of the coefficients. 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/
When I went back to the solution for your first post I could not remember why the Equal proc used sort - so I took it out. But, this is what is needed to do what you want. The only change is to the Equal proc. Here is the new procedure and the three examples:
EquivalenceClass := proc( A )
  local a, b, B, S, n, N, AddA, Equal, PrintA;
  uses combinat;
  AddA  := ()   -> add( A[i], i=args[1]):
  Equal := (a,b)-> evalb( sort(AddA(a))=sort(AddA(b)) ):
  PrintA := ()  -> add( a[i], i=args[1]):

  if nargs>1 then n := args[2] else n := 2 end if;
  N := {$1..nops(A)};
  S := choose(N,n);
  B := {};
  while nops(S)>0 do
    b,S := selectremove( s->Equal(s,S[1]), S );
    B := B union {b};
  end do;
  return seq( [map( PrintA, b),AddA(b[1])], b=B );
end proc:

A := [ [1,0,1,0], [1,1,0,0], [0,0,1,1], [0,1,0,1] ]:

EquivalenceClass( A );
   [{a[1] + a[2], a[1] + a[3], a[2] + a[4], a[3] + a[4]}, [2, 1, 1, 0]], 

     [{a[1] + a[4], a[2] + a[3]}, [1, 1, 1, 1]]

EquivalenceClass( A, 1 );
                  [{a[1], a[2], a[3], a[4]}, [1, 0, 1, 0]]

EquivalenceClass( A, 3 );
       [{a[1] + a[2] + a[3], a[1] + a[2] + a[4], a[1] + a[3] + a[4], 

         a[2] + a[3] + a[4]}, [2, 1, 2, 1]]
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/
When I went back to the solution for your first post I could not remember why the Equal proc used sort - so I took it out. But, this is what is needed to do what you want. The only change is to the Equal proc. Here is the new procedure and the three examples:
EquivalenceClass := proc( A )
  local a, b, B, S, n, N, AddA, Equal, PrintA;
  uses combinat;
  AddA  := ()   -> add( A[i], i=args[1]):
  Equal := (a,b)-> evalb( sort(AddA(a))=sort(AddA(b)) ):
  PrintA := ()  -> add( a[i], i=args[1]):

  if nargs>1 then n := args[2] else n := 2 end if;
  N := {$1..nops(A)};
  S := choose(N,n);
  B := {};
  while nops(S)>0 do
    b,S := selectremove( s->Equal(s,S[1]), S );
    B := B union {b};
  end do;
  return seq( [map( PrintA, b),AddA(b[1])], b=B );
end proc:

A := [ [1,0,1,0], [1,1,0,0], [0,0,1,1], [0,1,0,1] ]:

EquivalenceClass( A );
   [{a[1] + a[2], a[1] + a[3], a[2] + a[4], a[3] + a[4]}, [2, 1, 1, 0]], 

     [{a[1] + a[4], a[2] + a[3]}, [1, 1, 1, 1]]

EquivalenceClass( A, 1 );
                  [{a[1], a[2], a[3], a[4]}, [1, 0, 1, 0]]

EquivalenceClass( A, 3 );
       [{a[1] + a[2] + a[3], a[1] + a[2] + a[4], a[1] + a[3] + a[4], 

         a[2] + a[3] + a[4]}, [2, 1, 2, 1]]
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 took much longer than I expected. First, I must apologize for the erroneous solution posted yesterday. I must have been dreaming when I claimed it worked. The solution I am about to present does not use aliases (because they do not work within a procedure - this is well documented, so not a bug). Instead, I work directly with the indices into the list of vectors, checking for equivalence with the procedure Equal (which uses AddA) and then printing the equivalence classes with the procedure PrintA. Note that these procs can be implemented without using the convert( ..., `+` ) trick. So, now, here is my latest solution:
EquivalenceClass := proc( A )
  local a, b, B, S, n, N, AddA, Equal, PrintA;
  uses combinat;
  AddA  := ()   -> add( A[i], i=args[1]):
  Equal := (a,b)-> evalb( AddA(a)=AddA(b) ):
  PrintA := ()  -> add( a[i], i=args[1]):

  if nargs>1 then n := args[2] else n := 2 end if;
  N := {$1..nops(A)};
  S := choose(N,n);
  B := {};
  while nops(S)>0 do
    b,S := selectremove( s->Equal(s,S[1]), S );
    B := B union {b};
  end do;
  return seq( [map( PrintA, b),AddA(b[1])], b=B );
end proc:
And, in action:
A := [ [1,0,1,0], [1,1,0,0], [0,0,1,1], [0,1,0,1] ]:

EquivalenceClass( A );
 [{a[1] + a[3]}, [1, 0, 2, 1]], [{a[2] + a[3], a[1] + a[4]}, [1, 1, 1, 1]], 

   [{a[1] + a[2]}, [2, 1, 1, 0]], [{a[2] + a[4]}, [1, 2, 0, 1]], 

   [{a[3] + a[4]}, [0, 1, 1, 2]]

EquivalenceClass( A, 1 );
  [{a[1]}, [1, 0, 1, 0]], [{a[2]}, [1, 1, 0, 0]], [{a[4]}, [0, 1, 0, 1]], 

    [{a[3]}, [0, 0, 1, 1]]

EquivalenceClass( A, 3 );
[{a[1] + a[2] + a[3]}, [2, 1, 2, 1]], [{a[1] + a[2] + a[4]}, [2, 2, 1, 1]], 

  [{a[1] + a[3] + a[4]}, [1, 1, 2, 2]], [{a[2] + a[3] + a[4]}, [1, 2, 1, 2]]
Note, in the output, that I also include the common sum of all members of each equivalence class. If you don't want this, remove the AddA(b[1]) and the surrounding []'s in the return command at the end of EquivalenceClass. Hoping this is closer to being 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/
This took much longer than I expected. First, I must apologize for the erroneous solution posted yesterday. I must have been dreaming when I claimed it worked. The solution I am about to present does not use aliases (because they do not work within a procedure - this is well documented, so not a bug). Instead, I work directly with the indices into the list of vectors, checking for equivalence with the procedure Equal (which uses AddA) and then printing the equivalence classes with the procedure PrintA. Note that these procs can be implemented without using the convert( ..., `+` ) trick. So, now, here is my latest solution:
EquivalenceClass := proc( A )
  local a, b, B, S, n, N, AddA, Equal, PrintA;
  uses combinat;
  AddA  := ()   -> add( A[i], i=args[1]):
  Equal := (a,b)-> evalb( AddA(a)=AddA(b) ):
  PrintA := ()  -> add( a[i], i=args[1]):

  if nargs>1 then n := args[2] else n := 2 end if;
  N := {$1..nops(A)};
  S := choose(N,n);
  B := {};
  while nops(S)>0 do
    b,S := selectremove( s->Equal(s,S[1]), S );
    B := B union {b};
  end do;
  return seq( [map( PrintA, b),AddA(b[1])], b=B );
end proc:
And, in action:
A := [ [1,0,1,0], [1,1,0,0], [0,0,1,1], [0,1,0,1] ]:

EquivalenceClass( A );
 [{a[1] + a[3]}, [1, 0, 2, 1]], [{a[2] + a[3], a[1] + a[4]}, [1, 1, 1, 1]], 

   [{a[1] + a[2]}, [2, 1, 1, 0]], [{a[2] + a[4]}, [1, 2, 0, 1]], 

   [{a[3] + a[4]}, [0, 1, 1, 2]]

EquivalenceClass( A, 1 );
  [{a[1]}, [1, 0, 1, 0]], [{a[2]}, [1, 1, 0, 0]], [{a[4]}, [0, 1, 0, 1]], 

    [{a[3]}, [0, 0, 1, 1]]

EquivalenceClass( A, 3 );
[{a[1] + a[2] + a[3]}, [2, 1, 2, 1]], [{a[1] + a[2] + a[4]}, [2, 2, 1, 1]], 

  [{a[1] + a[3] + a[4]}, [1, 1, 2, 2]], [{a[2] + a[3] + a[4]}, [1, 2, 1, 2]]
Note, in the output, that I also include the common sum of all members of each equivalence class. If you don't want this, remove the AddA(b[1]) and the surrounding []'s in the return command at the end of EquivalenceClass. Hoping this is closer to being 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/
You are correct. One of my last edits must have broken something. I'll recheck this and then repost the correct version. My comment about remembering this problem was not a complaint. It was just an observation. I found it interesting that in the original thread I believe I was the person to suggest using combinat[choose], but this idea did not occur to me when you reposted the problem. More to follow shortly, I hope, 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 are correct. One of my last edits must have broken something. I'll recheck this and then repost the correct version. My comment about remembering this problem was not a complaint. It was just an observation. I found it interesting that in the original thread I believe I was the person to suggest using combinat[choose], but this idea did not occur to me when you reposted the problem. More to follow shortly, I hope, 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 knew this problem sounded familiar. With your last post I remembered the exact exchange from almost 3 weeks ago. It's pretty easy to modify that approach to handle your new version of the problem. Here is the update to the procedure to find equivalence classses that I introduced in the original discussion:
EquivalenceClass := proc( A )
  local b, B, s, S, n, Add, Equal, NPrintF;
  uses combinat;
  Add := () ->  convert( [args], `+` ):
  Equal := (a,b)->evalb(sort(a)=sort(b)):
  NPrintF := proc()
    if nargs<2 then return args end if;
    nprintf("%a+%A",args[1],NPrintF(args[2..-1]));
  end proc:
  if nargs>1 then n := args[2] else n := 2 end if;
  S := choose(A,n);
  B := {};
  while S<>{} do
    for b in B do
      S := select( s->not Equal(Add(s[]),Add(b[])), S );
    end do;
    if S<>{} then B := B union {S[1]} end if;
  end do;
  return map( b->NPrintF(b[]), B );
end proc:
Note that this is more specific to this problem in that I no longer have the binary operation (Add) and the equality test (Equal) as arguments to EquivalenceClass. The main functional difference is allowing an extra argument to specify the number of terms to consider when forming the equivalence classes. If omitted, the default is 2. This uses combinat[choose] as suggested earlier in this discussion. The main programming difference was the work to create a way to use nprintf with an arbitrary number of arguments. My solution is the recursive procedure NPrintF (and learning the difference between format a and A). Note that the implementation of Add needed to be modified to accept any number of arguments. The conversion of a list (or set) to a sum (`+`) was useful here. Now, the moment you have all been waiting for, some examples of EquivalenceClass in action:
aa := [ [1,0,1,0], [1,1,0,0], [0,0,1,1], [0,1,0,1] ]:
alias( seq( a[k]=aa[k], k=1..nops(aa) ) );
                           a[1], a[2], a[3], a[4]
A := { a[k] $ k=1..nops(aa) };
                          {a[1], a[2], a[3], a[4]}
EquivalenceClass( A );
     {a[1]+a[2], a[1]+a[3], a[1]+a[4], a[2]+a[3], a[2]+a[4], a[3]+a[4]}
EquivalenceClass( A, 1 );
                          {a[1], a[2], a[3], a[4]}
EquivalenceClass( A, 3 );
      {a[2]+a[3]+a[4], a[1]+a[2]+a[3], a[1]+a[2]+a[4], a[1]+a[3]+a[4]}
I added a fourth vector to the initial set (to make things more interesting). Note also that the construction of the aliases and the set of vectors, A, are now much more general. This could easily be included in the EquivalenceClass procedure. The new procedure would accept aa as the input. I'll leave this for you, or someone else. I hope this is another step forward, 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 knew this problem sounded familiar. With your last post I remembered the exact exchange from almost 3 weeks ago. It's pretty easy to modify that approach to handle your new version of the problem. Here is the update to the procedure to find equivalence classses that I introduced in the original discussion:
EquivalenceClass := proc( A )
  local b, B, s, S, n, Add, Equal, NPrintF;
  uses combinat;
  Add := () ->  convert( [args], `+` ):
  Equal := (a,b)->evalb(sort(a)=sort(b)):
  NPrintF := proc()
    if nargs<2 then return args end if;
    nprintf("%a+%A",args[1],NPrintF(args[2..-1]));
  end proc:
  if nargs>1 then n := args[2] else n := 2 end if;
  S := choose(A,n);
  B := {};
  while S<>{} do
    for b in B do
      S := select( s->not Equal(Add(s[]),Add(b[])), S );
    end do;
    if S<>{} then B := B union {S[1]} end if;
  end do;
  return map( b->NPrintF(b[]), B );
end proc:
Note that this is more specific to this problem in that I no longer have the binary operation (Add) and the equality test (Equal) as arguments to EquivalenceClass. The main functional difference is allowing an extra argument to specify the number of terms to consider when forming the equivalence classes. If omitted, the default is 2. This uses combinat[choose] as suggested earlier in this discussion. The main programming difference was the work to create a way to use nprintf with an arbitrary number of arguments. My solution is the recursive procedure NPrintF (and learning the difference between format a and A). Note that the implementation of Add needed to be modified to accept any number of arguments. The conversion of a list (or set) to a sum (`+`) was useful here. Now, the moment you have all been waiting for, some examples of EquivalenceClass in action:
aa := [ [1,0,1,0], [1,1,0,0], [0,0,1,1], [0,1,0,1] ]:
alias( seq( a[k]=aa[k], k=1..nops(aa) ) );
                           a[1], a[2], a[3], a[4]
A := { a[k] $ k=1..nops(aa) };
                          {a[1], a[2], a[3], a[4]}
EquivalenceClass( A );
     {a[1]+a[2], a[1]+a[3], a[1]+a[4], a[2]+a[3], a[2]+a[4], a[3]+a[4]}
EquivalenceClass( A, 1 );
                          {a[1], a[2], a[3], a[4]}
EquivalenceClass( A, 3 );
      {a[2]+a[3]+a[4], a[1]+a[2]+a[3], a[1]+a[2]+a[4], a[1]+a[3]+a[4]}
I added a fourth vector to the initial set (to make things more interesting). Note also that the construction of the aliases and the set of vectors, A, are now much more general. This could easily be included in the EquivalenceClass procedure. The new procedure would accept aa as the input. I'll leave this for you, or someone else. I hope this is another step forward, 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 recommend using a conditional based on whether there are any "indeterminates" (parameters) in the result from LinearSolve:
MyProc := proc( ... )
  ...
  sol := LinearSolve( A, b );
  if indets(sol)={} then return end if;
  ...
end proc;
The conditional could also be written as nops(indets(sol))=0. 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 recommend using a conditional based on whether there are any "indeterminates" (parameters) in the result from LinearSolve:
MyProc := proc( ... )
  ...
  sol := LinearSolve( A, b );
  if indets(sol)={} then return end if;
  ...
end proc;
The conditional could also be written as nops(indets(sol))=0. 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, the numerical constant know as pi is represented, in Maple, as Pi (the capital P is essential). When you used pi, you got the Greek letter pi, with no numerical value. Here's proof of what I am saying:
a:= 10*(2*pi/360); # 10 degrees
                                    1    
                                    -- pi
                                    18   
evalf( a );
                              0.05555555556 pi
a := 10*(2*Pi/360);
                                    1    
                                    -- Pi
                                    18   
evalf( a );
                                0.1745329252
The output using this value of a is the same as before, but if you pass this output through evalf (EVALuate Floating-point) you will get the result I think you are seeking:
Rx := Matrix([[1,0,0],[0,cos(a),-sin(a)],[0,sin(a),cos(a)]]);
Rx.b;
evalf(Rx.b);
Maple's representation of vectors and matrices does not make it easy to copy the contents of a vector to MaplePrimes, but converting the final result to a list shows how the entries have changed:
convert(Rx.b,list);
     [          /1    \         /1    \        /1    \         /1    \]
     [10, 20 cos|-- Pi| - 30 sin|-- Pi|, 20 sin|-- Pi| + 30 cos|-- Pi|]
     [          \18   /         \18   /        \18   /         \18   /]
convert(evalf(Rx.b),list);
                       [10., 14.48670973, 33.01719614]
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/
First, the numerical constant know as pi is represented, in Maple, as Pi (the capital P is essential). When you used pi, you got the Greek letter pi, with no numerical value. Here's proof of what I am saying:
a:= 10*(2*pi/360); # 10 degrees
                                    1    
                                    -- pi
                                    18   
evalf( a );
                              0.05555555556 pi
a := 10*(2*Pi/360);
                                    1    
                                    -- Pi
                                    18   
evalf( a );
                                0.1745329252
The output using this value of a is the same as before, but if you pass this output through evalf (EVALuate Floating-point) you will get the result I think you are seeking:
Rx := Matrix([[1,0,0],[0,cos(a),-sin(a)],[0,sin(a),cos(a)]]);
Rx.b;
evalf(Rx.b);
Maple's representation of vectors and matrices does not make it easy to copy the contents of a vector to MaplePrimes, but converting the final result to a list shows how the entries have changed:
convert(Rx.b,list);
     [          /1    \         /1    \        /1    \         /1    \]
     [10, 20 cos|-- Pi| - 30 sin|-- Pi|, 20 sin|-- Pi| + 30 cos|-- Pi|]
     [          \18   /         \18   /        \18   /         \18   /]
convert(evalf(Rx.b),list);
                       [10., 14.48670973, 33.01719614]
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/
First 54 55 56 57 58 59 60 Last Page 56 of 76