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

Get what curve?

My example does plot the two vertical lines through the points found by fsolve.

I did this by giving a parametric representation of the lines ( x=x0, y=y0, z=t for -infinity < t < infinity ).

They should show up as thick red lines in the plot (but they looked more like discrete point when I plotted them).

Is this not what you wanted?

Might you want the line segment in the x-y plane that joins the two points found by fsolve? If so, let us know. There are several ways this can be found, including with pointplot3d, as follows:

q:=2*10^(-7)*8.99*10^9;
v:=q/sqrt((x-a)^2+(y-b)^2);

g:=-diff(v,x);

charges:=subs(a=-1/2,b=0,v)+subs(a=1/2,b=0,v)-subs(a=0,b=3/4,v);

P1 := plot3d(charges,x=-3..3,y=-4..4,view=-1e4..1e4);

xg:=-diff(charges,x);
yg:=-diff(charges,y);
fsolve({xg=0,yg=0},{x,y},{x=-1..1},{y=-1..0});
L1 := eval( [x,y,t], % );
p1 := eval( [x,y,0], %% );

fsolve({xg=0,yg=0},{x,y},{x=-1..1},{y=0..4});
L2 := eval( [x,y,t], % );
p2 := eval( [x,y,0], %% );

P3 := pointplot3d( [p1,p2], style=line, color=green, thickness=3 ):
P2 := map( spacecurve, [L1,L2], t=-1e4..1e4, thickness=4, color=red )[]:
display( [P1,P2,P3], axes=boxed, view=[-3..3,-4..4,-1e4..1e4] );

You might need to rotate the plot a little to be able to see the green line in the x-y plane. Part is above the potential surface and part is below.

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

This sounds like something that should be addressed by Technical Support.

They will likely want to know exactly what sort of modifications you made to the theme. It's pretty clear to me this is the most likely source of your problems.

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

This sounds like something that should be addressed by Technical Support.

They will likely want to know exactly what sort of modifications you made to the theme. It's pretty clear to me this is the most likely source of your problems.

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?

P1 := plot3d(charges,x=-3..3,y=-4..4,view=-1e4..1e4);

xg:=-diff(charges,x);
yg:=-diff(charges,y);

fsolve({xg=0,yg=0},{x,y},{x=-1..1},{y=-1..0});
L1 := eval( [x,y,t], % );

fsolve({xg=0,yg=0},{x,y},{x=-1..1},{y=0..4});
L2 := eval( [x,y,t], % );

P2 := map( spacecurve, [L1,L2], t=-1e4..1e4, thickness=4, color=red )[]:
display( [P1,P2], axes=boxed, view=[-3..3,-4..4,-1e4..1e4] );

Hopes 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

Like this?

P1 := plot3d(charges,x=-3..3,y=-4..4,view=-1e4..1e4);

xg:=-diff(charges,x);
yg:=-diff(charges,y);

fsolve({xg=0,yg=0},{x,y},{x=-1..1},{y=-1..0});
L1 := eval( [x,y,t], % );

fsolve({xg=0,yg=0},{x,y},{x=-1..1},{y=0..4});
L2 := eval( [x,y,t], % );

P2 := map( spacecurve, [L1,L2], t=-1e4..1e4, thickness=4, color=red )[]:
display( [P1,P2], axes=boxed, view=[-3..3,-4..4,-1e4..1e4] );

Hopes 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

It sounds as though you are wanting to plot a curve, not necessarily a (straight) line.

If you have an explicit formula for the electric field, say E (as a vector of 3 functions of x,y, and z) then I expect you could use the implicitplot3d command.

with( LinearAlgebra ):
with( plots ):
E := < a, b, c >;
e := simplify( Norm(E,2) ) assuming real;
implicitplot3d( e=0, x=-3..3, y=-3..3, z=-3..3, axes=boxed );

But, this is will do better when the vector field vanishes on a surface, not a curve.

It would be best if you could identify parametric equations for the curve where your electric field is zero. Then you would use spacecurve to plot the curve:

spacecurve( [1+t,2+t,3+t] ,t=-1..1, axes=boxed );

Or, are you hoping to use the actual data points in Maple's plot of the electric potential? In this case there would be more work to do. I'll wait for you to confirm exactly which form you are working with before giving the details of this approach.

And, finally, to answer your final question: YES, you would use the display command to combine two (or more) plots. Just be sure you remember to load the plots package (or use the long form of the command name).

P1 := plot3d( ... ):
P2 := spacecurve( ... ):
display( [P1,P2] );

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

It sounds as though you are wanting to plot a curve, not necessarily a (straight) line.

If you have an explicit formula for the electric field, say E (as a vector of 3 functions of x,y, and z) then I expect you could use the implicitplot3d command.

with( LinearAlgebra ):
with( plots ):
E := < a, b, c >;
e := simplify( Norm(E,2) ) assuming real;
implicitplot3d( e=0, x=-3..3, y=-3..3, z=-3..3, axes=boxed );

But, this is will do better when the vector field vanishes on a surface, not a curve.

It would be best if you could identify parametric equations for the curve where your electric field is zero. Then you would use spacecurve to plot the curve:

spacecurve( [1+t,2+t,3+t] ,t=-1..1, axes=boxed );

Or, are you hoping to use the actual data points in Maple's plot of the electric potential? In this case there would be more work to do. I'll wait for you to confirm exactly which form you are working with before giving the details of this approach.

And, finally, to answer your final question: YES, you would use the display command to combine two (or more) plots. Just be sure you remember to load the plots package (or use the long form of the command name).

P1 := plot3d( ... ):
P2 := spacecurve( ... ):
display( [P1,P2] );

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

Axel,

I think I can explain the terminology. It's a mix of Maple and Windows.

The "File tab" is the "File" menu within the Maple GUI. I believe the comment is that trying to Open an existing worksheet (document) or create a New worksheet (document) causes Maple to crash.

The "theme" is a Windows setting. To control the theme, position the cursor on an empty space on your desktop. Hold the right-hand mouse button (for a right-handed mouse) and select the Properties entry from the popup window that appears. The "Display Properties" window will appear. The first tab in this window is called Themes. You can select a theme from the drop-down menu, or find one on the web (maybe for a price).

I have not had the problem reported by the original poster in this thread. I have, however, seen that the theme setting changes the way Maple displays buttons, etc. in maplets.

Even though this does not help the OP, I hope this explanation of terminology is helpful to some.

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

Axel,

I think I can explain the terminology. It's a mix of Maple and Windows.

The "File tab" is the "File" menu within the Maple GUI. I believe the comment is that trying to Open an existing worksheet (document) or create a New worksheet (document) causes Maple to crash.

The "theme" is a Windows setting. To control the theme, position the cursor on an empty space on your desktop. Hold the right-hand mouse button (for a right-handed mouse) and select the Properties entry from the popup window that appears. The "Display Properties" window will appear. The first tab in this window is called Themes. You can select a theme from the drop-down menu, or find one on the web (maybe for a price).

I have not had the problem reported by the original poster in this thread. I have, however, seen that the theme setting changes the way Maple displays buttons, etc. in maplets.

Even though this does not help the OP, I hope this explanation of terminology is helpful to some.

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 if you view your problem as a DAE (Differential Algebraic Equation)?

Here is how that might proceed:

restart;
with( plots ):
G := 6.63*10^(-11):
R := 2*10^6:
M := 3.6*10^30:
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):
de4 := diff(v(t), t) = r(t)*omega(t)^2-G*M*A(t);
               d                      2                 20     
              --- v(t) = r(t) omega(t)  - 2.386800000 10   A(t)
               dt                                              

Note that the last ODE contains the, as yet, unspecified A(t). We define this via an algebraic equation:

eq0 := A(t) = piecewise(r(t) > R, 2*r(t)/R^3, 1/r(t));
                      /                         1                 1  \
      A(t) = piecewise|2000000 < r(t), ------------------- r(t), ----|
                      \                4000000000000000000       r(t)/

Now, specify the initial conditions and ask Maple to provide a numeric solution to this DAE IVP:

ic := r(0) = 4*10^6, v(0) = 1, omega(0) = 1, theta(0) = 1:
s := dsolve([de1, de2, de3, de4, eq0, ic], numeric, maxstep = 0.1e-5, maxfun = 5000000);
proc(x_rkf45_dae)  ...  end;

Note that the argument of this procedure is x_rkf45_dae (not the usual x_rkf45). This indicates that Maple detects a DAE and is using the version of RKF45 that supports DAEs.

The best way to understand this solution is by several different plots:

odeplot( s, [[t,r(t)],[t,R]], t=0..1 );   # graphs of r(t) and R (to show different parts of A's def)
odeplot( s, [r(t)*cos(theta(t)),r(t)*sin(theta(t))], t=0..1 ); # the polar plot you are interested in

The first time r(t) crosses R is around t=0.1377. To investigate what happens at that time consider:

odeplot( s, [[t,A(t)],[t,2*r(t)/R^3],[t,1/r(t)]], t=0.135..0.140, style=[line,point,point] );

Notice that r(t)<R at only 3 points, between t=0.1376 and t=0.1378. This essentially instantaneous drop is not visible in the plot of r(t) on a longer time interval (as above).

If you try to go much beyond t=1, the numerical solution believes it encounters a singularity. Notice that I have reduced the maximum step and increased the number of function calls in the dsolve command. You might find additional tweaking of these enable you to get a solution over a longer time interval.

I still have a hard time believing this problem is correctly formulated. First of all, the function A(t) is not continuous.

You might also benefit from rescaling the problem so that all components of the solution have the same magnitude. This will help with the displaying of results, as well as some of the numerics.

I hope this is helpful and allows you to move closer to a final solution.

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 if you view your problem as a DAE (Differential Algebraic Equation)?

Here is how that might proceed:

restart;
with( plots ):
G := 6.63*10^(-11):
R := 2*10^6:
M := 3.6*10^30:
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):
de4 := diff(v(t), t) = r(t)*omega(t)^2-G*M*A(t);
               d                      2                 20     
              --- v(t) = r(t) omega(t)  - 2.386800000 10   A(t)
               dt                                              

Note that the last ODE contains the, as yet, unspecified A(t). We define this via an algebraic equation:

eq0 := A(t) = piecewise(r(t) > R, 2*r(t)/R^3, 1/r(t));
                      /                         1                 1  \
      A(t) = piecewise|2000000 < r(t), ------------------- r(t), ----|
                      \                4000000000000000000       r(t)/

Now, specify the initial conditions and ask Maple to provide a numeric solution to this DAE IVP:

ic := r(0) = 4*10^6, v(0) = 1, omega(0) = 1, theta(0) = 1:
s := dsolve([de1, de2, de3, de4, eq0, ic], numeric, maxstep = 0.1e-5, maxfun = 5000000);
proc(x_rkf45_dae)  ...  end;

Note that the argument of this procedure is x_rkf45_dae (not the usual x_rkf45). This indicates that Maple detects a DAE and is using the version of RKF45 that supports DAEs.

The best way to understand this solution is by several different plots:

odeplot( s, [[t,r(t)],[t,R]], t=0..1 );   # graphs of r(t) and R (to show different parts of A's def)
odeplot( s, [r(t)*cos(theta(t)),r(t)*sin(theta(t))], t=0..1 ); # the polar plot you are interested in

The first time r(t) crosses R is around t=0.1377. To investigate what happens at that time consider:

odeplot( s, [[t,A(t)],[t,2*r(t)/R^3],[t,1/r(t)]], t=0.135..0.140, style=[line,point,point] );

Notice that r(t)<R at only 3 points, between t=0.1376 and t=0.1378. This essentially instantaneous drop is not visible in the plot of r(t) on a longer time interval (as above).

If you try to go much beyond t=1, the numerical solution believes it encounters a singularity. Notice that I have reduced the maximum step and increased the number of function calls in the dsolve command. You might find additional tweaking of these enable you to get a solution over a longer time interval.

I still have a hard time believing this problem is correctly formulated. First of all, the function A(t) is not continuous.

You might also benefit from rescaling the problem so that all components of the solution have the same magnitude. This will help with the displaying of results, as well as some of the numerics.

I hope this is helpful and allows you to move closer to a final solution.

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 apologies for misreading your inequality. I see now that I was missing one factor of (1+x), and this is critical.

For your actual inequality it appears to hold for all 0 < p < p* where p* is just a little larger than 3 (maybe 3.08).

I do not know how to prove this symbolically.

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 apologies for misreading your inequality. I see now that I was missing one factor of (1+x), and this is critical.

For your actual inequality it appears to hold for all 0 < p < p* where p* is just a little larger than 3 (maybe 3.08).

I do not know how to prove this symbolically.

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 exporting the worksheet to a "Maple input" file. (You can find Export under the File menu.)

This file will,  by default, have .mpl as its extension. It's just a text file that can be edited with Notepad, TextPad, ....

This is the file you should be using as your input, not the .mw file (although it would be nice if the worksheet could be read).

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 exporting the worksheet to a "Maple input" file. (You can find Export under the File menu.)

This file will,  by default, have .mpl as its extension. It's just a text file that can be edited with Notepad, TextPad, ....

This is the file you should be using as your input, not the .mw file (although it would be nice if the worksheet could be read).

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 16 17 18 19 20 21 22 Last Page 18 of 76