Carl Love

Carl Love

28025 Reputation

25 Badges

12 years, 307 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

If your data is still accessible via the %, %%, or %%% operators (i.e., if it is one the three most-recent results), then it can't be garbage collected. Observe:

 

restart:

base:= kernelopts(bytesalloc);

8781824

A:= LinearAlgebra:-RandomMatrix(4000,4000):

kernelopts(bytesalloc) - base;

128004096

unassign('A'): gc();

kernelopts(bytesalloc) - base;

128004096

x: y: #Clear the %-%%-%%% stack.

gc():

kernelopts(bytesalloc) - base;

0

 

 

Download bytesalloc.mw

 

@tomet The norm of a vector is nonnegative. So the plot will have all the "mountains" proceeding along the positive z-axis. Maple is making no mistake in this. What you want is either the Divergence of the vector field or perhaps the Norm times the signum of the Divergence.

plot3d(VectorCalculus:-Divergence(E), x= -10..10, y= -10..10, view= -2*10^10..2*10^10);

I got quick results (a few seconds) using method= rosenbrock, and reasonable results (20-30 seconds) with method= mebdfi. The default method= rkf45 takes a very long time; I didn't let it finish. All other methods, including all the classical methods, produce garbage results, oddly enough.

Sol:= dsolve(
     #All initial conditions set to 0.
     {eq||(1..5), seq(q[k](0)=0, k= 1..5), seq(D(q[k])(0)=0, k= 1..3)},
     numeric, maxfun= 0, method= rosenbrock
):

plots:-odeplot(Sol, [seq([T,q[k](T)], k= 1..5)], T= 0..10);

plots:-odeplot(Sol, [seq([T,q[k](T)], k= 2..5)], T= 0..0.1);

A:= Matrix([[a,b,c], [d,e,f], [g,h,i]]);

#90-degree clockwise rotation:
ArrayTools:-FlipDimension(A,1)^%T;

#Flip over horizontal axis:
ArrayTools:-FlipDimension(A,1);

All eight can be built from those two.

add(nops(select(`=`, convert(n, base, 10), 5)), n= 1..1000);

     300

This can be gotten around by using delay-evaluation quotes:

seq('`if`'(a[i] < b[i], c[i], d[i]), i= 1..10);

seq has special evaluation rules, and `if` has what I call very special evaluation rules. I don't think it's a bug; it turns out to be useful sometimes. Here's an example of the rule:

restart:
A:= a:  B:= 2:
F:= _a-> eval(`if`(A,B,C), a= _a):
F(true);
                               B
eval(B);
                               2

There was a discussion of this in the MaplePrimes archives, but all the Answers and Replies disappeared!!!

There is a discontinuity plainly visible on the south polar contour of the sample plot that you posted. This extends from the pole along a meridian. So, this is not the fault of my code.

Areas of "no data" are bright red. These include a large, nearly circular area at the north pole and the evenly spaced streaks roughly perpendicular to the equator. It should be possible to use some interpolation to get rid of the streaks.

Here are precisely the steps that I used to import your data:

  1. Using FireFox, I clicked on the hyperlink to your data file that you posted.
  2. Using FireFox, I clicked Edit -> Select All and used Control-C to copy.
  3. I opened a new Notepad document on my desktop.
  4. I clicked in the blank document and used Control-V to paste.
  5. I positioned my cursor at the beginning of the first line of actual data.
  6. I hit backspace to delete the first (blank) line of the document.
  7. I saved the document to my desktop as globedata.txt.

Then I executed the following worksheet.


restart:

G:= ImportMatrix(
     "C:/Users/Carl/desktop/globedata.txt",
     source= delimited, delimiter= " ", datatype= float[8]
):

G:= G[.., 3]: #We only need column 3.

(M,m):= (max,min)(G):

#Set "no data" values to the max.
map[inplace](x-> `if`(x=m, M, x), G):

#Linear rescaling to 0..1.

m:= min(G):

map[inplace](x-> (x-m)/(M-m), G):

#Make into a matrix.

G:= ArrayTools:-Alias(G, [360, 180]):

#Plot it.

subsop(

     [1,2]= COLOR(HUE, ListTools:-Flatten(convert(G, listlist))[]),

     plot3d(

          1, theta= 0..2*Pi, phi= 0..Pi,

          grid= [360,180], coords= spherical, color= theta,

          style= patchnogrid, lightmodel= none, axes= none,
          viewpoint= [circleleft, frames= 50]

     )

);

 


Download globeplot.mw

The following procedure uses foldl to build the appropriate nested add command, and then evaluates it. Thus, this procedure maintains all of the efficiency associated with add.

CartProdAdd:= proc(f::{name,appliable}, L::seq(`..`))
local A, _i, V:= _i||(2..nargs);
     eval(foldl(A, f(V), V=~ L), A= add)
end proc:

Example of its use:

CartProdAdd(f, 1..2, 3..4);

The contour plot on a sphere (or on any plottable surface) can be achieved if the data can be arranged in a grid. That is, you need to have an evenly spaced set T of values of theta and an evenly spaced set P of values of phi such that there is a value of D for every element of the cartesian product T x P.  And T P must cover the globe. It may require a significant amount of interpolation for you to achieve this grid.

Let's suppose that the grid is stored in a matrix G:= Matrix(nops(T),nops(P)) each entry of which is a value of D. We do a linear rescaling of G to the range 0..1. Then we take a basic plot of the unit sphere in spherical coordinates with a specific coloring and replace the color matrix by G, like this:


# Linear rescaling to 0..1:
map[inplace](evalf, G):  M:= max(G):   m:= min(G):
map[inplace](x->
(x-m)/(M-m), G):

subsop(
     [1,2]= COLOR(HUE, ListTools:-Flatten(convert(G, listlist))[]),
     plot3d(
          1, theta= 0..2*Pi, phi= 0..Pi,
          grid= [nops(T),nops(P)], coords= spherical, color= theta,
          style= patchnogrid, lightmodel= none
     )
);

If you have trouble with the interpolation, then post your data.

Edit: Text and code updated with the linear rescaling idea due to a bug discovered by Markiyan Hirnyk in the Reply "Yet output" below.

Change the definitions of A, B, C, d, and E to

A:= unapply(int(1/y*exp(-C1*y^2), y= N..M), C1,N,M) assuming N > 0, M > 0;
B:= unapply(int(exp(-C1*y^2), y= N..M), C1,N,M);
C:= unapply(int(y*exp(-C1*y^2), y= N..M), C1,N,M);
d:= unapply(int(y^2*exp(-C1*y^2), y= N..M), C1,N,M);
E:= unapply(int(y^3*exp(-C1*y^2), y= N..M), C1,N,M);

Doing this, the Minimize returns a value instantly.

Maple's 2D input allows you to define functions using the syntax

F(x,y):= ...;

This is not robust, and it makes it difficult to share your code. So please break that habit and define functions like this:

F:= (x,y)-> ...;

The reason that I used unapply to define the functions is because all of the integrals can easily be done symbolically by Maple. Using evalf(int(...)) is a waste.

 

Use an assuming clause to make a temporary assumption during the solving phase, like this:

solve(...) assuming w::real;

If you use an extra space at the beginning of the name, then it will never be interpretted as an infix operator.

` &Delta;&epsilon;`(T,z);

While it's true that the extra space does show up in the typesetting, it usually doesn't look that bad.  For example,

x*` &Delta;&epsilon;`(T,z);

Of course, rather than repeatedly typing the weird name, I would first do

de:= ` &Delta;&epsilon;`:

and then exclusively use de for input.

The Direct Search optimization package does constrained global optimization, among other things. It is available for free download here.

Suppose that M is your matrix, a..b is your range of horizontal values, and c..d is your range of vertical values. Then do

plots:-surfdata(
     M, a..b, c..d,
     shading= zgrayscale, orientation= [-90,0], lightmodel= NONE
);

First 297 298 299 300 301 302 303 Last Page 299 of 395