I assume you have a system of ordinary differential equations with initial (or boundary) conditions and  you are using odeplot to graph their solution.

I know of no way to "subtract one odeplot from another". This needs to be expressed mathematically. I'm guessing that you will need to come up with a new system of ODEs (or new initial conditions) that describes the situation you are trying to study.

If you can tell us more about your specific problem, that will help us come up with a reasonable approach to suggest.

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

Try this:

DEplot(eq, y(x), x = -1 .. 1, y = -1 .. 3, [y(0)=1], colour = red, linecolor=green );

Notice that I changed the viewing window. Without this you won't see the solution curve to the right of x=0. Also, there is no need to solve for the derivative the way you did; just give DEplot the differential equation. (You should have written y as y(x) on the RHS, but Maple only issues a warning about this.)

If you want more solution curves, just put more initial conditions in the list, separated by commas.

The online help for DEplot contains some good examples showing how to create this kind of plot. See ?DEplot .

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 are fitting to a parabola, your general form will be y=a*x^2+b*x+c and the fit finds values of a, b, and c.

If you are fitting to a parabola that passes through the point (x0,y0), the general form will be y=a*(x-x0)^2+b*(x-x0)+y0 and the data can be used to find the parameters a and b.

You also want to have the parabola have a vertex at the point (x0,y0). This requires an additional modification to the general form and removes one more parameter. On the chance that this is a homework assignment, I'm not going to give this form. If you follow my reasoning, and know a little about parabolas, this won't be hard.

I hope this helps,

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

Without seeing exactly what you are doing it is difficult to explain what you are reporting. Here's a simple illustration to show that this idea can work.

with( ArrayTools ):
A := Array( 1..5, 1..5, 1..5, (i,j,k)->i+k ):
IsEqual( A[1][2], A[2][5] );
                                    false
IsEqual( A[2][2], A[2][5] );
                                    true

Did you, by chance, forget to load the ArrayTools package?

Does this help?

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

Change sum to add and all goes as you desire.

f := a*b;
B[1] := a;
B[2] := b;

add((diff(f, B[j]))^2, j = 1 .. 2);
                                    2    2
                                   b  + a 

You should use sum (or Sum) only when you have an indefinite sum that you expect can be rewritten without the summation.

See the help for ?add and ?sum for more explanations of the differences in these commands. You are not the first person to encounter this type of a problem with sum. Search MaplePrimes for more topics on this subject.

I do not know how to control the color of text in a MaplePrimes post. I worry more about changing the format to "formatted" to distinguish Maple input (and output) from what I have typed.

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

Have you looked at a regular 3D plot of your function? Have you seen the singularity at the origin? For example:

plot3d( -z*(1/(x^2+z^2)^(3/2)-1), x=-3..3, z=-3..3, grid=[50,50], axes=boxed);

The default window extends from (about) z=-90 to z=90. When you ask for 10 contours, Maple shows the following level curves:

seq( z, z=-90..90, 180/9 );
                 -90, -70, -50, -30, -10, 10, 30, 50, 70, 90

There are not very many points on these contours, and those that are all have abs(x) and abs(z) less than 0.5. This contributes to the poor resolution in the contourplot. To illustrate, consider this contour plot on an even smaller domain:

contourplot( f, x=-0.5..0.5, z=-0.5..0.5, grid=[50,50], contours=10, numpoints=100, axes=boxed );

While the plot looks essentially the same as the first contour plot, note that the plot is displayed on an even smaller domain.

You aren't seeing smooth curves because you are zoomed in too far on the plot. Increasing the grid does not help because it leads to sampling points even closer to the origin - which increases the vertical range that has to be covered by the 10 contours.

What to do?

A long time ago it was suggested to me that I use contourplot3d instead of contourplot, and my experience has shown that this was wise advise. I would apply it here as well. The contourplot3d command is more efficient and provides more information than contourplot. A fundamental difference is that contourplot3d returns a 3d plot (use orientation=[-90,0] to view the 3d plot from above with the axes in the standard orientation). To get started, look at these 2 plots:

f := -z*(1/(x^2+z^2)^(3/2)-1):
contourplot3d( f, x=-3..3, z=-3..3, grid=[50,50], contours=10, numpoints=100, axes=boxed, orientation=[-90,0] );
contourplot3d( f, x=-3..3, z=-3..3, grid=[50,50], contours=10, numpoints=100, axes=boxed );

The first plot is similar to the contourplot plot, except that it shows the full viewing window with x and z in [-3,3]. The second plot is the same, except the initial orientation is not set to obscure the 3d nature of the plot. Either plot can be rotated.

What is really needed is to convince Maple to show lower level curves on the surface. For this we supply an explicit list of contour values. Say

contourplot3d(f, x=-3..3, z=-3..3, grid=[50,50], contours=[$-10..10], numpoints=100,
                 axes=boxed, orientation=[-90,0]);

Remember that you can rotate this plot to see exactly how the levels stack up to show the surface. (The corresponding plot with contourplot is a total mess.)

I encourage you to consult the online help for contourplot3d (?contourplot3d) to see additional options that might further improve your plot.

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

Your code has so many errors that I hesitate to guess the real problem that you are trying to solve. Remember to include * (or at least a space) for each multiplication and that all parentheses must be matched. If you have done this work in Maple, just copy and paste from your Maple worksheet to MaplePrimes. Also, be careful about mixing upper and lower case; Maple is case sensitive (Dsol and dsol are different objects).

There is still an issue about how to integrate a solution to an ODE. The key is that you need to define the function H so that it returns unevaluated when the argument is not numeric:

HH := proc(y)
  if not type(y,numeric) then 
    return( 'procname'(y) );
  else
    eval( F(x), dsol(y) )
  end if;
end proc:

Then, you could find:

h2 := int( HH(y), y=-10..10 );
evalf( h2 );

There are many other ways to do this, including setting up a new ODE for the integral that you want.

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

You have described your problem perfectly. All you to do is have faith that Maple can handle a recursive call to your procedure - it can! Here is how I would modify your procedure:

addList := proc (a)::integer; local i, sum;
             description "add a list of digits";
#             sum := 0;
#             for i to length(a) do
#               sum := sum+convert(a, base, 10)[i]
#             end do;
             sum := add(i,i=convert(a,base,10));
             if length(sum)>1 then
               addList(sum)
             else
               return(sum)
             end if;
           end proc;

Note also that I replaced your loop to find the sum of the digits by a single call to the add command (see ?add). Here's how this procedure works in practice:

addList(999);
                                      9
debug( addList ):
addList( 999 );
{--> enter addList, args = 999
                                     27
{--> enter addList, args = 27
                                      9
<-- exit addList (now in addList) = 9}
                                      9
<-- exit addList (now at top level) = 9}
                                      9
addList( 99999099999999999999999986665432133399999999999787544556998888888888888888888888888888888888888888899999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999666666666666666666666666666699888888888777777777777888777777777777777777788888888888887777777888888888888888888888888999999777777777777777777777777777777777777777777777777777777777777999 );
     
{--> enter addList, args = 99999099999999999999999986665432133399999999999787544556998888888888888888888888888888888888888888899999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999666666666666666666666666666699888888888777777777777888777777777777777777788888888888887777777888888888888888888888888999999777777777777777777777777777777777777777777777777777777777777999
                                    3088
{--> enter addList, args = 3088
                                     19
{--> enter addList, args = 19
                                     10
{--> enter addList, args = 10
                                      1
<-- exit addList (now in addList) = 1}
                                      1
<-- exit addList (now in addList) = 1}
                                      1
<-- exit addList (now in addList) = 1}
                                      1
<-- exit addList (now at top level) = 1}
                                      1

Recursion can be very powerful, simplifying many algorithms. And, it is sufficiently efficient in Maple so that it is useful in many situations.

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'll show two more ideas that you might be able to apply towards your goal.

But, first, you really should get rid of D as a parameter in your ODE. Remember that D is the name of Maple's differentiation operator. Using D as a parameter will lead to unexpected problems. I note also that jakubi changed your D to DD (a common respelling of D for many Maple users - I will do the same).

Now, odeadvisor is a command in the DEtools package that attempts to classify an ODE in advance of actually solving it. Here:

ode1:=diff(phi(x),x$2)+1/x*diff(phi(x),x)-phi(x)/x^2=-1/DD*(q*x/2+P/(2*Pi*x)):
with( DEtools ):
odeadvisor( ode1 );
              [[_2nd_order, _exact, _linear, _nonhomogeneous]]

This is not as informative as we might have expected. The fact that it is exact is probably the most useful, but this is not a classification that most people expect, or know how to handle, for a second-order ODE.

Next, I'll point you to the interactive ODE analyzer. This is the ODE Analyzer under the list of Assistants (in the Tools menu), or start this maplet directly with the command:

dsolve[interactive](ode1);

In the window that appears, click on Solve Symbolically. Select a "method", decide if you want to see the Maple commands, then click Solve (in the second column). Back in the first window (which you can return to via the Back button) you can add initial and/or boundary conditions and give numerical values to the parameters. (As an idea of the problems caused by using just D as a parameter, try using the ODE Analyzer with the OP's original ODE.)

This does not really answer your question, but this combined with the information from increasing infolevel (as jakubi suggests) are probably the best you are going to be able to do at this time.

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

According to your post, your solve command is multiplied by assume. This does not make sense. I'm surprised this did not create an error.

I doubt you can use assume to find the critical points of your function. Instead, I suggest that you plot the derivative and use fsolve on intervals where you know there is exactly one root (of df):

restart;
f := x*sin(x^2)+1:
df := diff( f, x ):
solve( df=0, x );
                    /  2                                      \
              RootOf\_Z  - RootOf(tan(_Z) + 2 _Z), label = _L3/
plot( df, x=-1..3 );

fsolve( df=0, x=-1..1 );
                                     0.
fsolve( df=0, x=1..2 );
                                 1.355211129
fsolve( df=0, x=2..2.5 );
                                 2.194502750
fsolve( df=0, x=2.5..3 );
                                 2.813725766

Does this help?

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

Why do you believe there will be values of p for which the inequality is satisfied (for all 1<x<2)?

Consider, for x=1:

restart;
ineq :=  (1+x)^p > 2^p*(x^p+x);
                            p / p    \          p
                           2  \x  + x/ < (1 + x) 
ineq1 := eval( ineq, x=1 );
                                     p    p
                                  2 2  < 2 

Clearly, this inequality is never satisfied (for a real value of p).

solve( ineq1, p );

Here Maple's response, NULL, is correct. It is obvious that this inequality is not satisfied for any value of p. (But, you never said this was supposed to apply for x=1, just 1<x<2.)

The case for x=2 is a little more interesting:

ineq2 := eval( ineq, x=2 );
                               p / p    \    p
                              2  \2  + 2/ < 3 

For large values of p the LHS behaves like 4^p and the RHS like 3^p while for small (large negative) values of p the LHS behaves like 2^p and the RHS like 3^p. So, there should be exactly one value of p when the two sides are equal and the inequality will be satisfied for all p smaller than this value (which we estimate using fsolve):

solve( ineq2, p );
Warning, solutions may have been lost
fsolve( (lhs-rhs)(ineq2)=0, p );
                                -537.9214476

So, for p<-537.921 the inequality is satisfied when x=2.

More generally, let's define P(x) to be the value of p for which the inequality is satisfied as an equality:

ineqX := unapply( (lhs-rhs)(ineq), x ):
ineqX(2):
ineqX(3/2):
P := x -> fsolve( ineqX(x), p );
x -> fsolve(ineqX(x), p)
P(2);
                                -537.9214476
P(3/2);
                                -417.5208627

A (point) plot of this function can now be created:

plot( P, 0..2, -1000..0, numpoints=80, adaptive=false, style=point );

The real question is whether limit( P(x), x=1, right ) is finite, or negative infinity. You can zoom in on the function by restricting the domain:

plot( P, 1..1.1, -1000..0, numpoints=80, adaptive=false, style=point );

This plot suggests there is are values of p for which the inequality is true for all 1<x<2, something like p<-700.

This example shows why it can be dangerous to make too many conclusions based solely on graphs. Note that are some values of x for which this function gives abnormally small values for P. For example:

X := [1.099, 1.0999, 1.09999, 1.1, 1.10001, 1.1001, 1.101, 1.11];
         [1.099, 1.0999, 1.09999, 1.1, 1.10001, 1.1001, 1.101, 1.11]
P~(X);
   [-514.1640620, -13875.23882, -29258.97693, -28617.29360, -28345.31565, -24363.75838, -19674.29255, -685.1032675]

I find it rather difficult to believe the function changes this much. I welcome others to tell me why my intuition is wrong and what is really happening with this problem.

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 knew I had read about an event handler in Maple 13, but never had a chance to play with it until today. I don't claim to yet know how to make optimal use of this facility, but hope this post will encourage those with more knowledge to help me (and others) to improve our use of this intriguing new tool.

Here's the setup. I again use the two different versions of the last ODE in the system.

restart;
with( plots ):
G := 6.63*10^(-11)*10^(-33)*(1.49^3)(3600*(365*24))^12:
R := 2:
M := 3.6*10^37:
ic := r(0) = 7, v(0) = 1, omega(0) = 1, theta(0) = 1:
de1 := diff(r(t), t) = v(t):
de2 := diff(theta(t), t) = omega(t):
de3 := diff(omega(t), t) = -2*v(t)*omega(t)/r(t):
de4a := diff(v(t), t) = r(t)*omega(t)^2-G*M/r(t)^2;             # if r(t)<10
de4b := diff(v(t), t) = r(t)*omega(t)^2-2*G*M*r(t)/R^3;    # if r(t)>10

Note that I will create my own transition, at r(t)=10.

Since the IC starts with r(0)=7<10, the system initially follows the system with de4a:

s1 := dsolve([de1, de2, de3, de4a, ic], numeric, events=[[r(t)=10,halt]] ):
Warning, cannot evaluate the solution further right of .88446367, event #1 triggered a halt

Next, since the solution has r(t) increasing, the next stage uses de4b with the solution at this trigger point as the new initial conditions:

t1 := 0.88446367:
S1 := s1(t1):
ic1 := eval( S1[2..-1], S1[1] )[];
Warning, cannot evaluate the solution further right of .88446366, event #1 triggered a halt
            omega(0.884463667843913392) = 0.490000011199886454, 
              r(0.884463667843913392) = 10., 
              theta(0.884463667843913392) = 1.66681392288931306, 
              v(0.884463667843913392) = 5.06347526668841486

s2 := dsolve([de1, de2, de3, de4b, ic1], numeric, events=[[r(t)=10,halt]] ):
Warning, cannot evaluate the solution further right of 1.9665449, event #1 triggered a halt

The next stage in the solution uses de4a with new initial conditions:

t2 := 1.9665449:
S2 := s2(t2):
ic2 := eval( S2[2..-1], S2[1] )[];
Warning, cannot evaluate the solution further right of 1.9665448, event #1 triggered a halt
            omega(1.96654486510690220) = 0.490000097951576841, 
              r(1.96654486510690220) = 10.0000000000000018, 
              theta(1.96654486510690220) = 2.11175299682023621, 
              v(1.96654486510690220) = -5.06347546127608882

s3 := dsolve([de1, de2, de3, de4a, ic2], numeric, events=[[r(t)=10,halt]] ):

I conclude by putting the three pieces together in a single plot:

P1 := odeplot( s1, [[t,r(t)],[t,omega(t)],[t,theta(t)],[t,v(t)]], t= 0..3, legend=[r,omega,theta,v] ):
P2 := odeplot( s2, [[t,r(t)],[t,omega(t)],[t,theta(t)],[t,v(t)]], t=t1..3, legend=[r,omega,theta,v] ):
P3 := odeplot( s3, [[t,r(t)],[t,omega(t)],[t,theta(t)],[t,v(t)]], t=t2..4, legend=[r,omega,theta,v] ):

display( [P1,P2,P3] );

The plot shows all four components of the solution, on each of the three time intervals.

I imagine there are better ways to implement this transition. Maybe even a way to do this automatically within the event handling system (see ?dsolve,numeric,Events ). I look forward to seeing how this could be done.

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

Aman,

Your worksheet does not contain A, except in the output for de4. So, I assume that is the term in which A will appear. Here's what I have in mind. Define 2 odes:

de4a := diff(v(t), t) = r(t)*omega(t)^2-G*M/r(t)^2;
de4b := diff(v(t), t) = r(t)*omega(t)^2-2*G*M*r(t)/R^3;

Since r(0)=7>R, we start with the first version:

s := dsolve([de1, de2, de3, de4b, ic], numeric, method = classical[rk4], stepsize = 0.5e-2,
            maxfun = 500000, output = listprocedure):

To understand what is going on, let's look at the solution:

with( plots ):
odeplot( s, [[t,r(t)],[t,omega(t)],[t,theta(t)],[t,v(t)]], t=-1..10, legend=[r,omega,theta,v] );

Based on this picture, r(t)>6.8>R=2 for all t. If so, then the other ODE never comes into play. So, probably I have made a bad assumption about how A enters into your problems.

If you can tell us more precisely how your system of ODEs depends on A, I wil see if I can't provide a better solution to your problem.

As soon as I get back from the meeting I am supposed to be in right now I am also going to post another approach to your problem. This one will make use of the new "events handler" that is available within the numeric dsolve.

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

First, it seems that you have posted this question twice. Once under the title Working with Matrices and here. I'm going to flag the other posts to ensure that all responses appear here, where they are more appropriate.

I encourage you to use the LinearAlgebra package, not the linalg package. You are already using a Matrix, and the LinearAlgebra routines will be more useable and provide better results in the long run.

As I understand it, you are trying to find "values" of E for which your matrix A has determinant zero. Here's a direct solution to this:

restart;
with( LinearAlgebra ):
A:=Matrix([[alpha_O-E, beta_CO, 0, 0],[beta_CO, alpha_C-E, beta_CC, 0],[0, beta_CC, alpha_C-E, beta_CC],[0, 0, beta_CC, alpha_C-E]]):
detA := Determinant( A );
                3                    2                        2
 alpha_O alpha_C  - 3 alpha_O alpha_C  E + 3 alpha_O alpha_C E 

                               2            3                      2
    - 2 alpha_O alpha_C beta_CC  - alpha_O E  + 2 alpha_O E beta_CC 

               3            2  2              3                      2    4
    - E alpha_C  + 3 alpha_C  E  - 3 alpha_C E  + 2 E alpha_C beta_CC  + E 

         2        2          2        2            2                    2  2
    - 2 E  beta_CC  - beta_CO  alpha_C  + 2 beta_CO  alpha_C E - beta_CO  E 

             2        2
    + beta_CO  beta_CC 
solve( detA=0, E );
        /  4                            3
  RootOf\_Z  + (-alpha_O - 3 alpha_C) _Z 

       /                             2            2          2\   2   /
     + \3 alpha_O alpha_C + 3 alpha_C  - 2 beta_CC  - beta_CO / _Z  + \
                    2                    2          3                    2
  -3 alpha_O alpha_C  + 2 alpha_O beta_CC  - alpha_C  + 2 alpha_C beta_CC 

                2        \                     3          2        2
     + 2 beta_CO  alpha_C/ _Z + alpha_O alpha_C  + beta_CO  beta_CC 

                                2          2        2\
     - 2 alpha_O alpha_C beta_CC  - beta_CO  alpha_C /
allvalues( % );

The output from allvalues is HUGE! This is not surprising - you are looking for a symbolic solution of a quartic polynomial with coefficients that depend on four symbolic parameters.

Let's look a little more closely at your problem. What you are really asking for is the eigenvalues of the matrix

B := A+E*IdentityMatrix(4);
Eigenvalues( B );

The output from the Eigenvalues command is the same huge list of 4 expressions that was obtained above. (To be honest, I expected to get the RootOf form of this result.)

Now, what do you know about the four parameters in your problem. If you know something about the parameters, you can probably use this to simplify the results.

What do you need to know about the roots?

Because of the structure of your matrix (symmetric) we know the root (eigenvalues) are real-valued. It might be possible to find conditions when the matrix is positive definite (making all eigenvalues positive) - or not.

Before going further, I'll ask if you can give us any more information about your parameters.

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, the same approach works in your case:

ode1:= diff(x(t),t) = -4*x(t) + y(t) + 2;
ode2:= diff(y(t),t) = 2*x(t) - 2*y(t) +3;
plotPP := phaseportrait([ode1,ode2], [x(t), y(t)], t=0..5, [[x(0)=1, y(0)=-1], [x(0)=-1, y(0)=1]], arrows=large, x=-1..1, y=-1..1, linecolor = blue):
equil := eval( [ode1,ode2], [x(t)=x,y(t)=y] );
plotEQ := implicitplot( equil, x=-1..1, y=-1..1, color=green, thickness=2 ):
display( [plotPP,plotEQ] );

I think you will want to change the viewing window (in both plots), but I'll leave that for you to decide. If you don't like the way I obtained the equilibrium curves, replace this with your own method (or explicitly, if you so choose).

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