## 3068 Reputation

19 years, 295 days

Doug

`---------------------------------------------------------------------Douglas B. Meade  <><Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.eduPhone:  (803) 777-6183         URL:    http://www.math.sc.edu`

## plotting equilbrium solutions...

I assume you created the phase portrait using something like the following:

```with( DEtools ): with( plots ):
sys := { diff( x(t),t ) = -4*x(t)+  y(t)+2,
diff( y(t),t ) =  2*x(t)-2*y(t)+3 };
plotPP := DEplot( sys, {x(t),y(t)}, t=0..1, x=-2..5, y=-2..5, arrows=thin ):```

Then, the equilibrium curves can be found, and plotted, as follows:

```equil := eval( sys, [x(t)=x,y(t)=y] );
{0 = -4 x + y + 2, 0 = 2 x - 2 y + 3}
plotEQ := implicitplot( equil, x=-2..5, y=-2..5, color=blue )```

Lastly, put the two plots together using the display command:

```display( [plotPP,plotEQ] );
```

Adjust the colors, ranges, and othe properties as usual.

Does this help?

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## using Remove and SearchAll...

In the OP only the fractional part was being examined. This means the decimal point will be the first character of the string. In this case, you can get rid of the decimal point by simply dropping the first character:

p:=evalf((1+sqrt(5))/2-1,20);
0.61803398874989484820
ns:=convert(p,string)[2..-1];
"61803398874989484820"

More generally, to drop the decimal point from a general number you could use:

```with( StringTools ):
p:=evalf((1+sqrt(5))/2,20);             # this is the real golden ratio
1.6180339887498948482
ns:=Remove(c->evalb(c="."),convert(p,string));
"16180339887498948482"```

Now, the indices of all instances of a particular string can be obtained using:

```SearchAll("618",ns);
2
SearchAll("48",ns);
16, 18
```

To get more digits, just increase the second argument in the evalf command (or increase the Digits setting).

For more information see the ?StringTools , ?StringTools,Remove , ?StringTools,SearchAll , and ?Digits .

Doug

NOTE: Hyperlinks to the StringTool webpages added manually. Onle the one for Digits was generated automatically by MaplePrimes (as discussed here).

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## using Matrix to present table of result...

I achieve this result using Maple's Matrix data structure, as follows:

```interface( rtablesize=30 ):
dsys := {diff(v(t),t) = -x(t), diff(x(t),t)=v(t), x(0)=1, v(0)=0};
data := dsolve(dsys, 'numeric', 'output'=Array([seq(0..2, 0.1)]));
evalf[6]( < data[1,1], data[2,1] > );
```

All this does is create a matrix from the two elements returned by dsolve. Note that it's necessary to increase the setting of rtablesize if the matrix has more than 10 rows or columns; here, 30 is more than enough.

If only there was a way to copy this result into MaplePrimes.

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## diff of a list...

Do you really need the input to be a vector, actually a "Vector"? If you would be happy to express your vector as a list, then you could use just diff:

```v:=[5*t+1, (1/2*(-10))*t^2+50]:
diff( v, t );
[5, -10 t]
```

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## Like this?...

Like this?

```fsolve( hite(t)=11853 );
33.00303938
```

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## ensuring unique random numbers...

The odds of getting the same number from this sample space (1..10^20) is pretty slight. But, to be sure, you might want to use the following:

```UniqueRandom := proc( n, N )
local R, S;
R := rand( N );
S := {};
while nops(S)<n do
S := { seq(R(),i=1..5) }; print( S );
end do;
return( S );
end proc:
UniqueRandom( 5, 10 );
{3, 4, 8}
{2, 3, 5}
{1, 6, 7, 9}
{2, 4, 6}
{1, 3, 5, 7, 9}
{1, 3, 5, 7, 9}```

The idea is to simply repeat the random number selection unitl you do have 5 unique numbers. For real use, delete the print( S ). Then you will see something like:

```UniqueRandom( 5, 10^20 );
{3055679836444599521, 8706066334105081346, 30512373753908017532,
81071746988471090210, 99401568161139836807}
```

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## ensuring unique random numbers...

The odds of getting the same number from this sample space (1..10^20) is pretty slight. But, to be sure, you might want to use the following:

```UniqueRandom := proc( n, N )
local R, S;
R := rand( N );
S := {};
while nops(S)<n do
S := { seq(R(),i=1..5) }; print( S );
end do;
return( S );
end proc:
UniqueRandom( 5, 10 );
{3, 4, 8}
{2, 3, 5}
{1, 6, 7, 9}
{2, 4, 6}
{1, 3, 5, 7, 9}
{1, 3, 5, 7, 9}```

The idea is to simply repeat the random number selection unitl you do have 5 unique numbers. For real use, delete the print( S ). Then you will see something like:

```UniqueRandom( 5, 10^20 );
{3055679836444599521, 8706066334105081346, 30512373753908017532,
81071746988471090210, 99401568161139836807}
```

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## contours=...

Do you know about the contours= option to the contourplot and contourplot3d commands? Here's the information from Maple's online help for contourplot:

• There are eight contour levels by default. You can alter the number and the location of the contours used with the option contours = c where c is either an integer specifying the number of evenly spaced levels or a list of points [NOTE: should be numbers, not points] representing the contour levels.

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## Maplet Viewer...

But you can!

From the worksheet, under the File menu, select Export As ... and change the "Files of type:" to "Maplet (.maplet)", then save. Now, assuming file associations have not been corrupted, double click on the .maplet file just created. It launches using the Maplet Viewer (called the Maplet Launcher on Macs).

For an example, browse to http://maple.math.sc.edu/research/Irreduc.maplet. Since we have MapleNet, this maplet can be made available to anyone (even if they do not have Maple on their local computer) through the URL: http://maple.math.sc.edu/research/Irreduc.html .

It is useful to note that the .maplet file is just a plain text file with the Maple commands from the associated .mw file. Thus, this is an easy way to extract all of the commands from a worksheet.

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

Can you show use what you typed before getting this error message?

Without anything more, you are askimg us to find your needle in a haystack - blindfolded.

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## elementwise operators and solving comple...

The new (in Maple 13) elementwise operators can be used to find the real and imaginary parts of each element of a list (or vector or matrix ...):

```B := [ 3/13 - (2*I)/13, 3/2 - (5*I)/2, -(7/2) - 13*I, 1 - 3*I ];
[3    2     3   5      7                ]
[-- - -- I, - - - I, - - - 13 I, 1 - 3 I]
[13   13    2   2      2                ]
Re~(B);
[3   3  -7   ]
[--, -, --, 1]
[13  2  2    ]
Im~(B);
[-2  -5         ]
[--, --, -13, -3]
[13  2          ]
```

In previous versions of Maple this would have been done with map:

```map( Re, B );
[3   3  -7   ]
[--, -, --, 1]
[13  2  2    ]
```

Note that both Re and Im behave the same way. It's an error to write Re(B) or Im(B).

You have one equation, in x and y. This is not going to have a unique solution.

```solve( (x+1+I*(y-3))/(5+3*I) = 1+I, [x, y]);
[[x = -I y + (1 + 11 I), y = y]]
```

This means that for any value of y there is a value of x such that the equation is satisfied. Of course, most of these solutions are complex valued. The additional information that both x and y are real, pins down a specific solution. It would be nice if you could impose this constraint as follows, but this is not completely satisfactory:

```solve( { (x+1+I*(y-3))/(5+3*I) = 1+I, Im(x)=0 }, [x, y] );
[[x = -I RootOf(-11 + Re(_Z)) + (1 + 11 I), y = RootOf(-11 + Re(_Z))]]
evalf(%);
[[x = 1. + 0. I, y = 11. + 0. I]]
```

Note that adding Im(y)=0 to the list of equations produces the notorious "solutions may be lost" warning.

Here is how I would use Maple to produce a step-by-step solution to this problem:

```q1 := (x+1+I*(y-3))/(5+3*I) = 1+I;
/5    3   \
|-- - -- I| (x + 1 + I (y - 3)) = 1 + I
\34   34  /
q2 := isolate( q1, x );
x = (1 + 8 I) - I (y - 3)
s1 := solve( evalc(Im~(q2)), {y} );
{y = 11}
q3 := eval( q1, s1 );
/5    3   \
|-- - -- I| (x + (1 + 8 I)) = 1 + I
\34   34  /
solve( q3, {x} );
{x = 1}
```

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## convergence rate: visual and analytic...

Here are some things I would consider:

```q1 := n->(1+1/n)^n:
q2 := n->(1+1/n)^(n+1/2):
plot( [q1,q2], 1..infinity );
```

This plot shows that (as expected) both functions have the same limit (at infinity) and the second one (in green) approaches this limit faster than the first one (in red). The limits are:

```L1 := limit( q1(n), n=infinity );
exp(1)
L2 := limit( q2(n), n=infinity );
exp(1)```

A better way to look at the convergence rate is the following:

```plot( [abs(q1(n)-L1),abs(q2(n)-L2)], n=1..100 );
```

To find the actual convergence rate, look at the ratio of successive errors for the first sequence:

```r1 := (q1(n+1)-L1)/(q1(n)-L1);
(n + 1)
/      1  \
|1 + -----|        - exp(1)
\    n + 1/
---------------------------
n
/    1\
|1 + -|  - exp(1)
\    n/
```

and expand this as a series in n at infinity

```series( r1, n=infinity );
1    23      539      /1 \
1 - - + ----- - ------ + O|--|
n       2        3    | 4|
12 n    144 n     \n /

```

If you are not comfortable with an expansion at infinity, make the substitution n=1/h and expand the result about h=0 (then use h=1/n to express in n):

```series( eval(r1,n=1/h), h );
23  2   539  3    / 4\
1 - h + -- h  - --- h  + O\h /
12      144

```

Repeating the same steps for the second sequence:

```r2 := (q2(n+1)-L2)/(q2(n)-L2):
series( r2, n=infinity );
2   4     /1 \
1 - - + -- + O|--|
n    2    | 3|
n     \n /

```

shows that the second sequence converges twice as fast as the first sequence.

`Doug`
```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## bra-ket disassembly...

I don't know much about bra-kets, or the Physics package, but let's see if this is of use to you.

We start by loading the Physics package and defining the quantum operator:

```with( Physics ):
Setup(op=a);
* Partial match of  'op' against keyword 'quantumoperators'
[quantumoperators = {a}]

```

Then, the expression you include in your post in can be defined as

```ans := Bra(X).a[11].Ket(X) + Bra(X).a[12].Ket(Y) + Bra(X).a[13].Ket(Z)
+ Bra(Y).a[21].Ket(X) + Bra(Y).a[22].Ket(Y) + Bra(Y).a[23].Ket(Z)
+ Bra(Z).a[31].Ket(X) + Bra(Z).a[32].Ket(Y) + Bra(Z).a[33].Ket(Z);
Bracket(Bra(X), a[11], Ket(X)) + Bracket(Bra(X), a[12], Ket(Y))

+ Bracket(Bra(X), a[13], Ket(Z)) + Bracket(Bra(Y), a[21], Ket(X))

+ Bracket(Bra(Y), a[22], Ket(Y)) + Bracket(Bra(Y), a[23], Ket(Z))

+ Bracket(Bra(Z), a[31], Ket(X)) + Bracket(Bra(Z), a[32], Ket(Y))

+ Bracket(Bra(Z), a[33], Ket(Z))

```

This display shows a little of Maple's structure for this object. (It is displayed on the screen as you provided in your post.)

Next, I don't want to assume the terms are always in this order, so have to have a way of mapping from the bra and ket back to the appropriate index in the matrix:

```bk_index := proc(T)
local key;
key := table( [X=1, Y=2, Z=3] );
return ( key[op(op(1,T))], key[op(op(3,T))] ) = op(2,T);
end proc:
A := table( map( bk_index, convert(ans,list) ) ):

```

Then, the table A can be used to define the matrix you seek:

```M := Matrix( 3,3, (i,j)->A[i,j] );
```

Does this do what you need? If so, all of this could be put into a single proc. I'm sure others will have other ways to improve upon this, so won't give the combined procedure at this time (besides, I'm already late leaving the office).

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## would you settle for an explicit solutio...

As I looked at your example, I started to think that Maple just might be able to provide an explicit solution to your question. Here's my approach:

```restart;
PARAM := [ x0=1/2 ]:
s := x->x^2:
ode := diff( y(x), x\$2 ) - (c^2+s(x)*d^2)*y(x) = 0:
ic := y(b)=exp(-b), D(y)(b)=1/2:
sol := unapply( rhs(dsolve( {ode,ic}, y(x) )), (x,b,c,d) );
(x, b, c, d) ->
/       /           /   2          \
| (1/2) |           |  c    1   2  |
- |b      |-WhittakerW|- ---, -, b  d| exp(-b) d
\       \           \  4 d  4      /

/   2          \
|  c    1   2  |          2  2
+ 2 WhittakerW|- ---, -, b  d| exp(-b) b  d
\  4 d  4      /

/   2          \                        /   2          \
|  c    1   2  |          2             |  c    1   2  |
+ WhittakerW|- ---, -, b  d| exp(-b) c  - WhittakerW|- ---, -, b  d| d b
\  4 d  4      /                        \  4 d  4      /

/   2               \  \           /   2          \\/
|  c  - 4 d  1   2  |  |           |  c    1   2  ||
- 4 exp(-b) WhittakerW|- --------, -, b  d| d| WhittakerM|- ---, -, x  d||
\    4 d     4      /  /           \  4 d  4      //

//            /   2          \           /   2               \
||            |  c    1   2  |           |  c  - 4 d  1   2  |
||3 WhittakerW|- ---, -, b  d| WhittakerM|- --------, -, b  d| d
\\            \  4 d  4      /           \    4 d     4      /

/   2          \           /   2               \
|  c    1   2  |           |  c  - 4 d  1   2  |  2
- WhittakerW|- ---, -, b  d| WhittakerM|- --------, -, b  d| c
\  4 d  4      /           \    4 d     4      /

/   2          \           /   2               \  \       \
|  c    1   2  |           |  c  - 4 d  1   2  |  |  (1/2)|
+ 4 WhittakerM|- ---, -, b  d| WhittakerW|- --------, -, b  d| d| x     | +
\  4 d  4      /           \    4 d     4      /  /       /

/       /                   /   2          \
| (1/2) |                   |  c    1   2  |
|b      |-exp(-b) WhittakerM|- ---, -, b  d| d
\       \                   \  4 d  4      /

/   2          \
|  c    1   2  |  2  2
+ 2 exp(-b) WhittakerM|- ---, -, b  d| b  d
\  4 d  4      /

/   2          \
|  c    1   2  |  2
+ exp(-b) WhittakerM|- ---, -, b  d| c
\  4 d  4      /

/   2               \
|  c  - 4 d  1   2  |
+ 3 exp(-b) WhittakerM|- --------, -, b  d| d
\    4 d     4      /

/   2               \
|  c  - 4 d  1   2  |  2
- exp(-b) WhittakerM|- --------, -, b  d| c
\    4 d     4      /

/   2          \\           /   2          \\///
|  c    1   2  ||           |  c    1   2  || ||
- d b WhittakerM|- ---, -, b  d|| WhittakerW|- ---, -, x  d|| ||3
\  4 d  4      //           \  4 d  4      // \\

/   2          \           /   2               \
|  c    1   2  |           |  c  - 4 d  1   2  |
WhittakerW|- ---, -, b  d| WhittakerM|- --------, -, b  d| d
\  4 d  4      /           \    4 d     4      /

/   2          \           /   2               \
|  c    1   2  |           |  c  - 4 d  1   2  |  2
- WhittakerW|- ---, -, b  d| WhittakerM|- --------, -, b  d| c
\  4 d  4      /           \    4 d     4      /

/   2          \           /   2               \  \       \
|  c    1   2  |           |  c  - 4 d  1   2  |  |  (1/2)|
+ 4 WhittakerM|- ---, -, b  d| WhittakerW|- --------, -, b  d| d| x     |
\  4 d  4      /           \    4 d     4      /  /       /

dsol := D[2](sol):                                  # derivative of the solution wrt the 2nd variable (b)
eval( dsol( x0/3, x0/3, 10, 1 ), PARAM ):           # exact solution is too complicated to display here, but
eval( dsol( x0/3, x0/3, 10, 1 ), evalf(PARAM) );    # the numeric approximation is reasonable
-1.346481861
```

See Maple's help for information about the Whittaker M function. You might find a series approximation useful at this point, but I won't go there in this post.

If an explicit formula for the derivative was not available, here is how I would approach your problem:

```SOL := (bb, cc, dd) -> dsolve(eval({ic, ode}, [b = bb, c = cc, d = dd]), y(x), numeric, output = operator):
Sol := (a,b,c,d) -> eval( y(a), SOL(b,c,d) );
eval(('Sol'(x0/3, x0/3+db, 10, 1)-'Sol'(x0/3, x0/3, 10, 1))/db, PARAM);
-1.303773000
eval(('Sol'(x0/3, x0/3+db^2, 10, 1)-'Sol'(x0/3, x0/3, 10, 1))/db^2, PARAM);
-1.346400000```

The last calculation provides a more accurate estimate of the deriviative, agreeing with the value obtained from the analytic solution to 4 decimal places.

My approaches have completely avoided the interpolation step. The interpolant might not be a very good approximation to the actual solution surface. Polynomial interpolants have a tendency to exhibit wild oscillations between data points, which could really affect the quality of any derivative that you compute from the interpolant.

I hope something in here has been useful,

Doug

```---------------------------------------------------------------------
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu```

## evaluating a derivative...

The error message is telling you that dx2Data is being called when x2 has a value (1/6). This is a problem because it does not make sense to differentiate a function with respect to a constant. This value is coming from the fact that x0=1/2 and your call to dx2Data has second argument (1/3)*x0, or 1/6.

What to do? I don't have time to completely analyze your problem, but maybe you will have better results with this definition of dx2Data:

```dx2Data := (x1,x2,x3,x4) -> eval( diff(f(x1, x, x3, x4), x), x=x2 );
```

I hope this helps,

Doug

```---------------------------------------------------------------------