There is dsolve[interactive] (see ?dsolve,interactive ).

For example,

ode := y(x) = 2*x*diff(y(x),x,x);
dsolve[interactive]( ode );

This is not perfect, but might be of some help (or interest) to you or others.

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, for each row just change the diagonal entry to a number that is larger (in absolute value) than the sum of the absolute value of the other entries in the row. For example:

A := LinearAlgebra:-RandomMatrix( 4, 4 ):
for i to LinearAlgebra:-RowDimension(A) do
  A[i,i] := convert( abs~(A[i,..]), `+` );
end do:
A;

But, I doubt this is what you want. You probably want something that preserves some property of the matrix. If you can give us more information about what you are really trying to do, there's a better chance someone will have something constructive to suggest. (This is, of course, assuming this is not part of a homework assignment.)

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's what I get (with Maple 13.02):

with( DEtools ):
infolevel[dsolve] := 5:

ode := y(x) = 2*x*diff(y(x),x,x);
                                    / d  / d      \\
                         y(x) = 2 x |--- |--- y(x)||
                                    \ dx \ dx     //
odeadvisor( ode );
                             [[_Emden, _Fowler]]

dsolve( ode, y(x) );
Methods for second order ODEs:
--- Trying classification methods ---
trying a quadrature
checking if the LODE has constant coefficients
checking if the LODE is of Euler type
trying a symmetry of the form [xi=0, eta=F(x)]
checking if the LODE is missing 'y' 
-> Trying a Liouvillian solution using Kovacic's algorithm
<- No Liouvillian solutions exists
-> Trying a solution in terms of special functions:
   -> Bessel
   <- Bessel successful
<- special function solution successful
                             (1/2)        /    (1/2)  (1/2)\
                 y(x) = _C1 x      BesselI\1, 2      x     /

                           (1/2)        /    (1/2)  (1/2)\
                    + _C2 x      BesselK\1, 2      x     /

Notice that I increased the infolevel for dsolve. You can see a little more information if you replace the "dsolve" with "all" in the infolevel command.

Based on the solution that Maple returns, it seems likely to me that Maple is making a change of variable u=sqrt(2*x). Let's see how that might work:

PDEtools[dchange]( x=u^2/2, ode );
                            /   d          d  / d      \\
                            |  --- y(u)   --- |--- y(u)||
                            |   du         du \ du     /|
                   y(u) = u |- -------- + --------------|
                            |      2            u       |
                            \     u                     /
ode2 := (lhs-rhs=0)(simplify(%)*u);
                         / d      \   / d  / d      \\      
                u y(u) + |--- y(u)| - |--- |--- y(u)|| u = 0
                         \ du     /   \ du \ du     //      
odeadvisor( ode2 );
                   [[_2nd_order, _with_linear_symmetries]]
dsolve( ode2, y(u) );
Methods for second order ODEs:
--- Trying classification methods ---
trying a quadrature
checking if the LODE has constant coefficients
checking if the LODE is of Euler type
trying a symmetry of the form [xi=0, eta=F(x)]
checking if the LODE is missing 'y' 
-> Trying a Liouvillian solution using Kovacic's algorithm
<- No Liouvillian solutions exists
-> Trying a solution in terms of special functions:
   -> Bessel
   <- Bessel successful
<- special function solution successful
              y(u) = _C1 u BesselI(1, u) + _C2 u BesselK(1, u)

This confirms my conjecture. Maple recognized the original ODE as one that could be converted to Bessel's equation.

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

When you use -> you are creating a special type of Maple procedure. You can also the unapply command to convert an expression into a function.

For multiline functions you would use the proc ... end proc syntax (see ?proc ).

The help for -> (see ?-> or ?arrow ) includes the statement that:

A functional operator of the form:
                               vars -> result
is semantically equivalent to:
                               proc( vars) option operator, arrow; result end

In the ?arrow help page we get more details about the full meaning of the operator and arrow options: disables Maple simplification rules that add any non-local names in the operator expression to the parameter list (and is meaningless for modules).

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 think I can explain what Maple is telling you - but it would be best if you could give us the details of your specific situation. Also, what version of Maple are you using?

Maple uses a trailing tilde to indicate that a variable has some assumptions associated with it. In this case I'm guessing that the symbol i1 is assumed to be an integer (or positive integer). The way in which Maple displays assumed variables can be controlled via the GUI. Under the Tools menu, select the Display tab. The "Assumed variables:" can be set to No Annotation, Trailiing Tilde, or Phrase. No Annotation is clear, and you've already seen Trailing Tilde. Here is how Phrase would look:

assume( i, posint );
i;
                                      i
                            With assumptions on i

There is, to my knowledge, no meaning to tilde as a prefix to an equality. "Not equal" is <>.

In Maple 13 the Maple language was extended to allow the use of tildes to indicate the elementwise application of a procedure. For example,

floor( [ 1.1, 1.8, 2.1, -3.2 ] );
    Error, invalid input: floor expects its 1st argument, a1, to be of type algebraic,
    but received [1.1, 1.8, 2.1, -3.2]
floor~( [ 1.1, 1.8, 2.1, -3.2 ] );
                                [1, 1, 2, -4]

I hope this helps to answer your question. (If not, please do give us more information about your exact situation.)

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

Maple, by default, assumes undeclared variables are local.

restart;
f := proc()
  g := 3;
end proc;
Warning, `g` is implicitly declared local to procedure `f`
                      proc()  local g; g:=3;  end proc;
g;
                                      g
f();
                                      3
g;
                                      g

If you are bothered by the warning message, that can be turned off as follows:

interface( warnlevel=0 ):

(The default setting for warnlevel is 3, see also ?interface and ?warning )

Notice Maple's output from the definition of the proc. It includes the local declaration for g. So, if you really want an explicit local statement in your proc, just define the procedure without a local statement, then copy Maple's output back to a new input region (or an external file) and you will have the local statement. (You might have to work to reformat the procedure body.)

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

What you describe is a parametric plot. Here's how I would do it:

x := t->3*t^3;
y := t->4*x(t)+5;

plot( [x(t),y(t),t=-6..6] );

Note that this is correct. Your equations describe the straight line with slope 4 and y-intercept 5, from x=-648 to x=648.

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

No quaternions. The webpage you listed gives the 3d analog of the 2d iterations:

What's the formula of this thing?

There are a few subtle variations, which mostly end up producing the same kind of incredible detail. Listed below is one version. Similar to the original 2D Mandelbrot, the 3D formula is defined by:

z -> z^n + c

...but where 'z' and 'c' are hypercomplex ('triplex') numbers, representing Cartesian x, y, and z coordinates. The exponentiation term can be defined by:

{x,y,z}^n = r^n { sin(theta*n) * cos(phi*n) , sin(theta*n) * sin(phi*n) , cos(theta*n) }
...where:
r = sqrt(x^2 + y^2 + z^2)
theta = atan2( sqrt(x^2+y^2), z )
phi = atan2(y,x)

And the addition term in z -> z^n + c is similar to standard complex addition, and is simply defined by:

{x,y,z}+{a,b,c} = {x+a, y+b, z+c}

The rest of the algorithm is similar to the 2D Mandelbrot!

Here is some pseudo code of the above:

r = sqrt(x*x + y*y + z*z )
theta = atan2(sqrt(x*x + y*y) , z)
phi = atan2(y,x)

newx = r^n * sin(theta*n) * cos(phi*n)
newy = r^n * sin(theta*n) * sin(phi*n)
newz = r^n * cos(theta*n)

...where n is the order of the 3D Mandelbulb. Use n=8 to find the exact object in this article.

The formulae above is the higher power analogue to the first green object shown on the first page here, but also see the formulae on Paul Nylander's page, which corresponds to the second green object shown instead. Both produce amazing (and similar) objects when put to the higher powers.

This would not be difficult to implement in Maple.

Mandelbrot sets have been discussed previously on MaplePrimes. One of the best resources I have seen is a blog post by Samir Khan (http://www.maplesoft.com/blog/view.aspx?sid=6851). I'll bet this could be modified to produce some of the 3D plots of Mandelbulb sets.

See also the Julia Sets and Mandelbrot forum on MaplePrimes (this forum started earlier this month, so the references in it should still be valid).

Good luck!

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 isolve command can also be used to solve this problem.

restart;
isolve({a+b+c=1, a+d+e=2, b+d+f=2, c+e+f=1} );
 {a = _Z2, b = _Z1, c = 1 - _Z1 - _Z2, d = 2 - _Z1 - _Z2, e = _Z1, f = _Z2}

This gives the set of solutions in terms of arbitrary integers _Z1 and _Z2.

The full set of solutions when these are restricted to values 0 and 1 can then be found with

S := { seq( seq( (14), _Z1=[0,1] ), _Z2=[0,1] ) };
                {{a = 0, b = 0, c = 1, d = 2, e = 0, f = 0}, 

                 {a = 0, b = 1, c = 0, d = 1, e = 1, f = 0}, 

                 {a = 1, b = 0, c = 0, d = 1, e = 0, f = 1}, 

                 {a = 1, b = 1, c = -1, d = 0, e = 1, f = 1}}

But, two of these solutions involve values that are neither 0 nor 1, so they must be removed:

select( a->rhs~(a) subset {0,1}, S );
                 {{a = 0, b = 1, c = 0, d = 1, e = 1, f = 0}, 

                  {a = 1, b = 0, c = 0, d = 1, e = 0, f = 1}}

(I really like the new ~ operator for componentwise application of a command.)

This process could be automated if so desired. The trickiest part is probably handling arbitrary number of parameters, and possibly different names used by Maple, in the original solution set. I'll leave this for someone else to implement.

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

Here is my implementation of my understanding of your request.

genRequestDist := proc( N, t )
  local data, pts1, pts2;
  data := sort( RandomTools:-Generate( list( integer(range=1..t), N ) ) );
  pts1 := [seq( [i,data[i]], i=1..N )];
  pts2 := [seq( [i,data[i+1]-data[i]], i=1..N-1 )];
  plot( [pts1,pts2], style=point );
end proc;

This works fine for N=100000 and t=3600000. It does take Maple (13.02) a while to produce this plot of 200000 (-1) points, but it does not complain. However, with the random numbers selected from 1 to 3.6 million, the plot of the differences looks like a plot of zero. If you don't sort the data you will see more differences.

A more appealing plot (at least to me) is obtained by reducing the number of points (N) to, say, 100.

genRequestDist( 100, 3600000 );

You should also consider looking at the visualization commands provided in the Statistics package (see ?Statistics,Visualization ).

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

I'm going to guess that the original problems were caused by the use of 2D input. Note the "." in the definition of z and the "*" after Vector in the definition of R. (Also, there is no need for the VectorField command.) Here is how I would have adapted the original code:

with (VectorCalculus):

x:= (r,t) -> r*sin(t):
y:= (r,t) -> r*cos(t):
z:= (r,t) -> r^2*sin(t)*cos(t):

R := Vector([x(r, t), y(r, t), z(r, t)]);
r1 := diff(R,r);
r2 := diff(R,t);
dS := simplify( CrossProduct(r1,r2) );

Hoping 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

Think about what you are trying to do.

Given your equation:

EQN := V+b = (1/(V*(4*V-b)^3)-c*V-d)/(P*sqrt(T)*V^2);
                                    1                
                               ------------ - c V - d
                                          3          
                               V (4 V - b)           
                       V + b = ----------------------
                                       (1/2)  2      
                                    P T      V       

If we want to try to solve for V by hand, we would move all terms to one side

EQN2 := normal( EQN-rhs(EQN) );
          1            /    7    (1/2)       6    (1/2)         5    (1/2)  2
---------------------- \64 V  P T      + 16 V  P T      b - 36 V  P T      b 
 3          3    (1/2)                                                       
V  (4 V - b)  P T                                                            

         4    (1/2)  3    4  3    (1/2)             5         4           3  2
   + 11 V  P T      b  - b  V  P T      - 1 + 64 c V  - 48 c V  b + 12 c V  b 

        2  3         4         3           2  2        3\    
   - c V  b  + 64 d V  - 48 d V  b + 12 d V  b  - d V b / = 0

The denominator is zero when V=0 or V=b/4:

solve( denom(lhs(EQN2))=0, V );
                                    1    1    1  
                           0, 0, 0, - b, - b, - b
                                    4    4    4  

Maybe these are not physically realistic and can safely be ignored - or maybe not.

The numerator is a 7th degree polynomial in V:

EQN3 := sort( collect( numer( lhs(EQN2) ), V ) );
        (1/2)  7         (1/2)  6     /       (1/2)  2       \  5      3  
  64 P T      V  + 16 P T      V  b + \-36 P T      b  + 64 c/ V  - V b  d

       /      (1/2)  3                \  4
     + \11 P T      b  - 48 b c + 64 d/ V 

       /    (1/2)  4       2           \  3   /  3         2  \  2    
     + \-P T      b  + 12 b  c - 48 b d/ V  + \-b  c + 12 b  d/ V  - 1

There is not going to be a closed form for the solutions of this equation, unless you happen to get lucky. Maybe the origins of your problem can provide ideas about how a symbolic solution can be found, but it's more likely that you'll have to rely upon fsolve to find an approximate solution for any given values of T, P, b, c, and d. This can be automated as follows:

VVall := (TT,PP,bb,cc,dd) -> fsolve( eval(EQN3=0, [T=TT, P=PP, b=bb, c=cc, d=dd]), V, complex ):
VVall(1,1,1,1,1);
     -0.9959275806, -0.2027917024, -0.005040084299 - 0.9998682614 I, 

       -0.005040084299 + 0.9998682614 I, 0.2272988158 - 0.3190790839 I, 

       0.2272988158 + 0.3190790839 I, 0.5042018200

Note that this will return all 7 solutions to the equation (ignoring any possible problems when the denominator is simultaneously zero, as found above):

If it's known that V is positive (and real), we can create a procedure that returns these solutions:

VVpos := (TT,PP,bb,cc,dd) -> select( is, {VVall(TT,PP,bb,cc,dd)}, positive )[]:
VVpos(1,1,1,1,1);
                                0.5042018200

Depending upon your choices for the parameters, there could be more than one positive solution. If this is an issue for you, you will have to modify the selection process to find the correct value to have Maple return. Or, maybe you do know that there is only one positive real solution to this equation.

What do you want to do with your solution? If you want to see how V changes as a function of one of the parameters you will need to write a wrapper that returns unevaluated when the arguments are not all numeric:

VVpos(T,1,1,1,1);
Error, (in fsolve) T is in the equation, and is not solved for
VV := proc(TT,PP,bb,cc,dd)
  if not type( {_passed}, set(numeric) ) then
    return 'procname'( _passed )
  else
    VVpos(_passed)
  end if
end proc:
VV(T,1,1,1,1); VV(1,1,1,1,1);
                              VV(T, 1, 1, 1, 1)
                                0.5042018200

Now, you can create plots as follows:

plot( VV(T,1,1,1,1), T=1..10 );
plot3d( VV(T,P,1,1,1), T=0..10, P=0..10, style=surfacecontour, shading=zhue, axes=boxed );

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

Are you saying you want to connect the end of the curve (red) with the beginning of the (visible) line (black)?

If so, then here is what I suggest. Basically, I just find the last point on the curve (ptA) and the first point on the line (ptB) and connect them with a (blue) line. Note that I define x0 and have modified the creation of plot p2.

S1:=unapply(u*FresnelS(u)+((1/Pi)*cos((1/2)*Pi*u^2))-(1/Pi),u);
C1:=unapply(u*FresnelC(u)-((1/Pi)*sin((1/2)*Pi*u^2)),u);
with(plottools);
with(plots);
u0:=sqrt(0.25);
x0:=evalf(FresnelC(u0)+FresnelS(u0)*cot((3*Pi)/8));
p1:=display(plot([FresnelC(u),FresnelS(u),u=0..u0]));
p2:=display(plot(x-x0,x=x0..1,y=0..1,colour=black));
ptA := [FresnelC(u0),FresnelS(u0)];
ptB := [x0,0];
p3 := plot( [ptA,ptB], color=blue );
display(p1,p2,p3);

If this is not what you are requesting, please restate your request. I do not understand what you mean.

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 first problem is that there is no value for A. You try to give it a value with

A:=fsolve(T=(R0+R1)*(((Pi*C1(B))^2+(Pi*S1(B)+1)^2))^(1/2),B,0..1);

but this equation has no solution. To see this, look at the plot of the rhs and note that its minimum value is 2.2 (when B=0).

If I change the value of T to something larger than 2.2 (and less than about 4.4), then your code concludes by producing two separate plots. Here are a few ways in which you might want to combine these plots:

# Your original plots, given names for future reference
p1:=display(plot([-a0*FresnelC(u),-a0*FresnelS(u),u=u0..0]),c0,c1);
p2:=display(plot([a1*FresnelC(u),a1*FresnelS(u),u=0..u1]),c0,c1);

# The two plots superimposed  (note the use of scaling=constrained to ensure the circles look like circles)
display( [p1,p2], scaling=constrained );

# The two plots as an animation
display( [p1,p2], insequence=true, scaling=constrained );

# All in one command (with the two parametric curves given different colors)
display( plot( [[-a0*FresnelC(u),-a0*FresnelS(u),u=u0..0],
                [ a1*FresnelC(u), a1*FresnelS(u),u=0..u1]] ),
         c0, c1,
         scaling=constrained );

Much more is possible, but I hope this gives you something to start with.

But, first, you have to choose a value of T so that A can be assigned a value.

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

Like this?

delay := proc( t::numeric )
  local st;
  st := time();
  while time()-st<t do end do;
end proc

x := true;
delay(2);
x := false;
                                    true
                                         ... WAIT 2 SECONDS ...
                                    false

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