Carl Love

Carl Love

28095 Reputation

25 Badges

13 years, 100 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

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
);

@Chia

If I understand your Reply "Is it possible to realize the translation?" correctly, you want to translate (also known as "pan") and zoom without changing the viewing angle. (It's the same as doing a "fly-over" view with decreasing altitude.) The key to this is that the look be the projection into the xy plane of the location. (That is, if the location is [x,y,z], then the corresponding look is [x,y,0].) That way you're always looking straight down even if x and y are changing due to the translation. Here's an example:

(zr,zi):= (Re,Im)(Zeta(1/2+I*im)):

plots:-display(
     [
          plots:-spacecurve([im,zr,0], im= 0..60, numpoints= 1000, color= red),
          plots:-spacecurve([im,zi,0], im= 0..60, numpoints= 1000, color= blue)
     ],
     viewpoint= [
          location= [seq([k/5,0,610-2*k], k= 150..300)],
          look= [seq([k/5,0,0], k= 150..300)],
          upvector= [0,1,0]
     ],
     axes= normal, axis[1]= [tickmarks= 60],
     title= typeset("Zeros of ",zeta(1/2+i*y)), titlefont= [HELVETICA,14]             
);

Is that what you're talking about?

You can take an operation such as subtraction an append a tilde to the operator to make it into an "elementwise" operator, which means that the operation will be performed on each element in a container, such as a vector.

theta:= 90 -~ M(..,1);

A "naked" tilde can be used to apply any procedure to each element in a container:

phi:= (x-> `if`(x < 0, 360+x, x)) ~ (M(..,2));

All of the parentheses are required in the above. An alternate syntax is

phi:= map(x-> `if`(x < 0, 360+x, x), M(..,2));

See ?map and ?elementwise .

r=0 is a regular singular point of your differential equation, so you can't give initial conditions at 0.

The error arises because there is no way to compute the recursive function _C when the input is symbolic n. I recommend first using rsolve to get a nonrecursive expression for the coefficients, like this:


restart:

R1:= rsolve({_C(n) = -_C(n-2)/n/(n-1), _C(0)= C[0], _C(1)= C[1]}, {_C(n)});

((1/2)*C[0]+((1/2)*I)*C[1])*(-I)^n/GAMMA(n+1)-((1/2)*I)*(I*C[0]+C[1])*I^n/GAMMA(n+1)

convert(%, factorial);

((1/2)*C[0]+((1/2)*I)*C[1])*(-I)^n/factorial(n)-((1/2)*I)*(I*C[0]+C[1])*I^n/factorial(n)

simplify(%) assuming n::nonnegint;

-(1/2)*(I*exp(((1/2)*I)*n*Pi)*C[1]-I*exp(-((1/2)*I)*n*Pi)*C[1]-exp(((1/2)*I)*n*Pi)*C[0]-exp(-((1/2)*I)*n*Pi)*C[0])/factorial(n)

R2:= evalc(%) assuming n::nonnegint;

(cos((1/2)*n*Pi)*C[0]+sin((1/2)*n*Pi)*C[1])/factorial(n)

sum(R1*x^n, n= 0..infinity);

(1/2)*C[0]*exp(I*x)-((1/2)*I)*C[1]*exp(I*x)+(1/2)*C[0]/exp(I*x)+((1/2)*I)*C[1]/exp(I*x)

sum(R2*x^n, n= 0..infinity);

C[0]*cos(x)+C[1]*sin(x)

 


Download rsolve_sum.mw

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