I agree with the book's answer. Here's how I would perform this change of variables in Maple. (Note that you have to indicate that the coefficients p and q are known functions.) I also massage the output to make it easy to compare with what you say your book presents.

I know it's a little odd that this requires the use of the PDETools package, but it's good that there is only one command to handle all such transformations. (There used to be a Dchanbevar command in the DEtools package, but it's been replaced by the more powerful dchange in PDETools - see the online help for full description of dchange.)

I hope this has been helpful.

Download DE-ChangeOfVariable.mwDE-ChangeOfVariable.mw

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

Here are a few ways in which this function can be implemented in Maple using recursion.

The first approach is to simply give the general formula, and then specify the initial condition.

P := t -> m[tau]/(1-psi*(beta-1)) - psi*(beta-1)*P(t-1)/(1-psi*(beta-1)):
P(0) := -1:

and, making explicit use of the remember table (and performing some simplification of the value)

Q := proc(t) option remember;
  collect(
    simplify(
      m[tau]/(1-psi*(beta-1)) - psi*(beta-1)*Q(t-1)/(1-psi*(beta-1))
    ), psi
  )
end proc:
Q(0) := -1:

As for making a table of values, here is how I would do this (using a Matrix):

MakeTable := (f,n) -> < <`j`|`p(j)`>, Matrix( n, 2, (i,j) -> `if`(j=1,i,f(i)) ) >:
MakeTable(Q,3);

I hope this has been helpful. Let us know if you have more questions.

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

Does this work for you?

M10 := {[0,[1,4,2,1]],[1,[2,1,4,2]],[2,[3,6,4,3]],[3,[3,6,4,3]],[4,[4,2,1,4]],[5,[4,3,6,4]],[6,[4,3,6,4]],[7,[5,8,6,5]],[8,[5,8,6,5]],[9,[5,8,6,5]],[10,[6,4,3,6]]}:

Categorize := S -> seq( [ [op~(1,select( a->evalb(op(2,a)=b), S ))[]], b], b=op~(2,S) ):

Categorize( M10 );
[[0], [1, 4, 2, 1]], [[1], [2, 1, 4, 2]], [[2, 3], [3, 6, 4, 3]], [[4], [4, 2, 1, 4]], [[5, 6], [4, 3, 6, 4]], [[7, 8, 9], [5, 8, 6, 5]], [[10], [6, 4, 3, 6]]

 To explain, a little, about how this works:

1. op~(2,S)
   returns the set of all second components in the input set.
2. select( a-> evalb(op(2,a)=b, S )
   returns all of the elements of S with second component equal to b
3. op~(1, ... )
   returns the set of all of the first components of the elements of S with second component equal to b
4. [ ...[] ]
   converts the set to a list
5. seq( [...,, b], b=S )
   assembles the list for each second component with the first component being the list of of first components with the same second component

It would not surprise me if someone else comes up with an even more elegant solution, possibly with somethign from the ListTools package. But, I like my one-liner.

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

My approach is a little different - but equivalent: I also got a final "answer" of 247, 500, 271, 500, 1322, 4000.

Maple's numeric dsolve has the ability to work with parameters. This allows the IVP to be "solved" once and then reused with different values of the parameters.

It is still necessary to deal with the multiple uses of the name q. In general you should not try to use the names of the functions in your IVP as names for any quantity - even itself - elsewhere in that Maple worksheet. My approach is to use a capital letter, Q, but others might prefer a more descriptive name like, SOLq.

Here's my revision of your code:

restart;
varepsilon := 1/50; s := (1/10)*varepsilon; N := (2-1)/s; L := 8; omega := 1+cos(varepsilon*t);
F := 0; S := 0; T := 0; F1 := N-F; S1 := N-S; T1 := L*N-T;
sys := diff(q(t), t) = p(t), diff(p(t), t) = -omega*sin(q(t));
ics := q(0) = 1, p(0) = p0:
sol := dsolve({ics, sys}, numeric, output = operator, maxfun = 0, parameters=[p0]):
Q := rhs(sol[3]):
for n from 1 by s to 2 do
  Q(parameters = [p0 = n]);
  if Q(0) < Q(2*Pi/varepsilon) then F := F+1 else F := F end if;
  if Q(2*Pi/varepsilon) < Q(4*Pi/varepsilon) then S := S+1 else S := S end if;
  for k from 0 by 2 to L do
    if Q(k*Pi/varepsilon) < Q((k+2)*Pi/varepsilon) then T := T+1 else T := T end if
  end do
end do;
F, F1, S, S1, T, T1;

                                   [p0 = 1.]
                                   ...
                                   [p0 = 2.]
                                   247, 500, 271, 500, 1322, 4000

Note the use of a parameter, p0, in the IC, the addition of parameters=[p0] in the dsolve, and the definition of Q as the numerical value of the 3rd component of the solution of the IVP, i.e. q. All of this takes place outside the loop. Inside the loop all that has to be done is to specify a value of the parameter, and change all occurrences of q to Q.

In addition to making the loop run faster, the code is cleaner and it's easier (for me) to understand what's going on inside the loops.

My revised worksheet is attached, as count-DBM.mw.

I have to admit that I had a little trouble with your epsilon. It was not "epsilon" but "varepsilon". Anyone familiar with TeX will understand, but others - including most students - will find this pretty odd. Another reason I prefer to stick with 1D input.

I hope this has been 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

I am the author of the shoot procedure that you are attempting to use. This code is pretty dated, but it's still very popular. I hope I can help you to find a way to use shoot for your problem.

Without digging too deep, it appears to me as though your "initial conditions" (IC) include some values at t=0 and others at t=etainf (and your "boundary conditions" (BC) involve values at t=15). I think the error message is trying to tell you that it does not like the IC that give values at more than one time.

You appear to have 13 different functions in your model. Your IC include 16 equations, 11 give information about t=0 and 5 at t=etainf. The 11 conditions at t=0 involve all functions exctp i and r; the also involves 5 parameters: alpha, beta, gamma, delta, and s There are 2 coditions that involve f(0). Taking into account other values specified in IC, it appears to me as though these reduce to

6.2*(0.08...*s*beta) + 0.5*beta = 0

from which you actually have a specific value for s (or beta=0).

Likewise, the two conditions that involve t(0) reduce to

gamma = -0.1 + 0.1*alpha

which allows for the removal of another parameter.

I encourage you to look closely at your problem and simplify the problem as much as possible.

Remember that shoot is a numerical solver. You have to give values to every parameter before it can find a solution.

You might have to get a little creative to massage your problem into one for which shoot can be used. Try to gain experience and confidence by considering simpler problems.

Let us know if you have additional questions.

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

Sure, but I think you want to substitute the first equation into the derivative of the second equation:

eq1 := diff(VL(t),t) = iL(t)/C:
eq2 := diff(iL(t),t) = (V-R*iL(t)-VL(t))/L:
eq3 := eval( diff(eq2,t), eq1 );

If you want to bring everything back to the left-hand side of the equation:

eq4 := eq3-rhs(eq3);

Note that this preserves the object as an equation. If you had done lhs(%)-rhs(%) you would have only an expression, without the = 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

You defined Q as the difference between the solution and cos(x+t), in absolute value.

pds is the module; change your last line to

pds:-plot(x=1, t=0..5);

Some other plots that you might find interesting: 

pds:-plot([u,cos(x+t)], x=0, t=0..5 );
pds:-plot([u,cos(x+t)], x=0.5, t=0..5 );
pds:-plot([u,cos(x+t)], x=1, t=0..5 );

and

pds:-plot3d([u,cos(x+t)], x=0..1, t=0..5 );

As has already been pointed out, the Maple help for pdsolve and pdsolve,numeric pages are 

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

You might want to look at Maple's identify command.

For your problem it returns 40/23, which actually is a better approximation than your 7/4.

It's not clear to me that specifying the number of bits for the numerator is a meaningful means of controlling the output. It seems to me that the number of bits in the denominator is a better control parameter. In this case, maybe you don't want to allow a denominator larger than 10, so that 40/23 is not acceptable.

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

What is i? Maple has a hard time working with arbitrary-order derivatives. If you give i an integer value, you get something more in line with what you expect:

i := 2:
Fdiffi := [diff(subs(timefull, F[1]), t$i), diff(subs(timefull, F[2]), t$i)];

[ / / d / d \\ / d / d \\\
[p[1] p[6] |-|--- |--- x[1](t)|| - |--- |--- x[2](t)|||
[ \ \ dt \ dt // \ dt \ dt ///

/ d / d \\ / / d / d \\
- p[3] |--- |--- x[1](t)||, p[4] p[6] |-|--- |--- x[1](t)||
\ dt \ dt // \ \ dt \ dt //

/ d / d \\\ / d / d \\]
- |--- |--- x[2](t)||| - p[5] |--- |--- x[2](t)||]
\ dt \ dt /// \ dt \ dt //]

Even more, if you look at the Maple help for pochhammer you quickly see that for all reasonable values of i, these values are zero. If you tell Maple that i is a positive integer, then all is well:

Fdiffi := [diff(subs(timefull, F[1]), t$i), diff(subs(timefull, F[2]), t$i)] assuming i::posint;

[p[1] p[6] (-(diff(x[1](t), [t $ i])) - (diff(x[2](t), [t $ i]))) - p[3] (diff(x[1](t), [t $ i])), p[4] p[6] (-(diff(x[1](t), [t $ i])) - (diff(x[2](t), [t $ i]))) - p[5] (diff(x[2](t), [t $ i]))]

You really can't expect more from Maple. There are times when you are interested in a "fractional derivative" so Maple's response allows Maple to be an appropriate tool for that use as well.

You might also be interested in shortening your code to find this derivative, and having the result in a Vector (like the input):

Fdiffi := diff~(subs(timefull, F), t$i) assuming i::posint;

I hope this has been 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

If you add eq1 and eq2 you get an equation involving a x, y, a, and F(a) - but not F'(a).

Solving for F(a) yields:

sol := isolate(eq1+eq2, F(a));

Now, substituting this result back into eq1 (or eq2 or eq3) should allow you to see if this solution makes any sense. When I do this the result is so complicated that I can't really see what's going on.

The appearance of 10^(-53) in your equations is a concern - if you want to do any numerical work. If you can't rescale the equations to ensure all terms have similar orders of magnitude, you may have to work with floating-point numbers with a lot of significant digits. (The default is Digits:=10; you might need Digits:=100;. Even then you would have to do some work to convince yourself - and others - that Maple's results really do make sense.)

A better option could be to introduce a new parameter, c, whose value will be 10^(-53). For example, with this suggestion, solving eq1+eq2 for F(a) produces:

F(a) = (4*y+c*(1+a)^3)/(2*x^2-2*y^2+2)

Your "initial condition" of F(0)=10^(-6) means that:

1/1000000 = (4*y+c)/(2*x^2-2*y^2+2)

Solving this for y gives

Recall the value of c, and this can be plotted:

sol2 := eval(sol, [a = 0, F(a) = 10^(-6)]);
sol3 := isolate(sol2, y);
plot(eval(rhs(sol3), c = 10^(-53)), x = -10 .. 10);

Do not be misled by the fact that this plot looks a lot like a parabola. For x in [-10,10], the corresponding y values are essentially -2*10^6. (Look closely at the values on the tick marks.) Even the plot of log(-y) looks a lot like a parabola, with values between 14.508657738530 and 14.508657738545. (This is done with Digits:=10; I'd have to do some more work before I put too much confidence in these values.)

Now, this has been done without any reference to eq3, just eq1+eq2.

I hope this gives you some ideas about what you can do to refine your equations. If my analysis has been off-target, please let us know what you are hoping to get from your equations.

Here's my version of your worksheet:

differential_eq2.mw

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

I looked at this system of matrix equations:

C.(A+B) = A1                      (1)
C.(B+I) = A2                      (2)
B+C.A = A3                        (3)

From (2),  assuming C is invertible,

B = inv(C).A2 - I                 (4)

Then from (1), and (4),

A = inv(C).(A1-A2) + I         (5)

Substitute (4) and (5) into (3) to obtain a single matrix equation for C:

inv(C).A2 + C = A3+A2-A1   (6)

While I don't have any immediate insights into solving (6) for C, it is clear that solving this equation would then give you explicit values for A and B.

Maybe there are other approaches that allow you to solve for A or B. I'll let you, or others, explore these ideas.

My bottom line is that there's probably a better way to approach this problem than pure brute force.

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

 If the vector in your post is a zero vector, then each of a, b, and c would have to be zero, and so would sqrt( a^2 + b^2 + c^2 ) - so that the denominator would also be zero.

But, I don't see than any of a, b, or c is zero. As a result, I'm not sure what your question is.

One thought I have is that since a, b, and c is complex valued, do you want a^2 or abs(a)^2, etc., in the denominator? (But this won't affect whether your vector is the zero vector or not.)

I would not even attempt to simplify this by hand if I had Maple at my side to do this for me. You don't have to use Maple only as a black box. You can see lots of intermediate steps.

First is the evalc command. This will simplify almost any expression into the form A+B*I. For example, evalc( a ); would rationalize the denominator and simplify the result as much as possible.

I would evaluate the denominator once, manipulating it until I get it in a desired form. Maybe

d2 := abs(a)^2 + abs(b)^2 + abs(c)^2;
d2 := simplify( d2 );

And, to get an idea of the numerical value of this result:

evalf( d2 );

or of the actual vector:

evalf( <a,b,c>/sqrt(d2) );

If I've missed the point of your post, please accept my apologies. Can you try to restate your problem and what you want Maple to do for you? Posting an actual worksheet could be very 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

The use of return is not critical here. What is important is that a comma-separated collection of values is a single Maple object.

LTTS:=proc(ff)
           local ll,r,r1,r2,r3;
           ll:=rhs(ff)-lhs(ff);
           solve01(ll), solve02(ll), solve03(ll), solve04(ll), solve05(ll)
          end:

Depending upon your exact intentions for using this result, you might prefer to have these 5 values as a list or set. If so, simply enclose the comma-separated values in square brackets (list) or squiggly brackets (set):

LTTSlist:=proc(ff)
           local ll,r,r1,r2,r3;
           ll:=rhs(ff)-lhs(ff);
           [ solve01(ll), solve02(ll), solve03(ll), solve04(ll), solve05(ll) ]
          end:

If you do stick with the comma separated values, you might want to know that you can still assign each value to a separate Maple name - all in one command:

a,b,c,d,e := LTTS( f );

If you return a list, you can still do the multiple assignment but you have to strip away the square brackets:

a,b,c,d,e := LTTSlist( f )[];

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

I can replicate your GetRidOfDumbSolutions procedure in one line.

When I first started to learn Maple I approached the programming from a classical point-of-view, with lots of loops, counters, etc. I had to modify my programming style to exploit Maple's ability to work with mathematical objects in mathematical ways. And, when you use loops, you don't have to index everything - you can loop through the elements of a set or list or the terms in an expression.

For example, here is an implementation of your procedure that avoids all counters and indices:

GetRidOfDumbSolutions2 := proc( sols )
  local GoodSols, s;
  GoodSols := NULL;
  for s in sols do
   if not member(0,rhs~(s)) then GoodSols := GoodSols, s end if;
  end do;
  return [GoodSols]
end proc;

(My Maple is complaining about the local statement. If you have the same, just comment it out and let Maple automatically add it. This will be the subject of my next post.)

The keys to this implementation are

1. stepping through each element of the input without the need for an index
2. the use of the tilde operator to construct the set of RHS's for each solution
3. checking for the presence of a 0 on the RHS of any equation in a solution in terms of set operations
4. the construction of the output as a comma separated list, that is converted to a list only at the very end

Another step toward utilizing Maple's power is learning to use Maple's select, remove, and selectremove commands. In the current context, here is a one-line implementation of your procedure:

GetRidOfDumbSolutions3 := S -> remove( s ->member(0,rhs~(s)), S );

The first argument to remove is a Boolean function that is to be applied to each element of the second argument; elements that have a return value of TRUE are removed, all others are kept. (select is similar, but the opposite.)

To show that these are all equivalent to your original procedure, for your input:

A := solve([p[1]*p[2]*p[3] = q[1]*q[2]*q[3], p[1]+p[3] = q[1]+q[3], p[2]^2+p[3]^2 = q[2]^2+q[3]^2]):
S1 := GetRidOfDumbSolutions( [A] ):
S1L := convert(S1,list);
S2 := GetRidOfDumbSolutions2( [A] );
S3 := GetRidOfDumbSolutions3( [A] );

evalb( S1L = S2 );
 true
evalb( S1L = S3 );
 true

But, as I read your original post, does this do what you want? The first two solutions have p[3]=p[3], and the third one has q[1]=q[1]. If you really do want to get rid of things like these, it might be helpful to move everything to one side of the equation and then look to see if the equation is linear with no constant term?

If you need more help, please post more information, possibly including a more detailed example. You might also want to tell us why you think these solutions are "dumb".

I hope this has been 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

I'm not exactly sure what you are expecting, but I'll bet some of what I'll present will be useful.

Let's start by giving your sample output a name. (I note that it is a set of sets. Each inner set appears to contain three equations, giving values for k, l, and o.)

S := {{k = k, l = RootOf(_Z^2+_Z*k+k^2-1), o = -k-RootOf(_Z^2+_Z*k+k^2-1)}, {k = k, l = RootOf(_Z^2+_Z*k+k^2+1), o = -k-RootOf(_Z^2+_Z*k+k^2+1)}, {k = 0, l = 1, o = -1}, {k = 0, l = -1, o = 1}, {k = 1, l = 0, o = -1}, {k = 1, l = -1, o = 0}, {k = -1, l = 0, o = 1}, {k = -1, l = 1, o = 0}, {k = RootOf(_Z^2+1), l = 0, o = -RootOf(_Z^2+1)}, {k = RootOf(_Z^2+1), l = -RootOf(_Z^2+1), o = 0}}:

The allvalues command is useful for expanding Maple's RootOf into specific values. The output is pretty long, so I'll just comment that it's a comma-separated set of sets. But, the union of all of these sets should be the set of all solutinos (to something).

A := allvalues( S ):
A2 := map(op,{A});

A2 contains 14 elements; 4 of them involve imaginary values. The remaining 10 can be found with:

A3 := remove( has, A2, I );

Lastly, the results might be easier to analyze if they were presented in a table:

A4 := map( a -> eval([k,l,o],a), convert(A3,list) ):
Matrix(convert(A4,listlist));

I hope that when you see the results you will find something that interests you.

If not, please provide us with a better description of what you are trying to 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