This might of interest (esp. if you like to watch your machine grind), even with limited scope.

acer

Are you sure about eq2?

acer

Alec's code appeared earlier in a worksheet by John Oprea (appeared in the Application Center in 2001, but may date from around 1996 and appears to have been created in MapleV R4). They both used a wrapping call to evalf when initializing local `c` in the scalar iterating procedure. Alec used plots[densityplot] (more attractive) while John used plot3d (faster).

And John's code is very close to that which appeared in the Maple V -- Programming Guide (p.231 of my 1996 edition). That manual's example used exponent 4 rather than 2, producing a kind of "Mandelbroid" if I may use the term. Alexander Walz cites [1, 2] John's revision as early as 1996.

One drawback to using hue is that its shading cycles back, with both lower and upper extremes shaded red. But with a high maximal iteration limit most inputs end up relatively near the extremes. So a lot of the plot is red, and it's difficult to distinguish visually which points escaped early and which did not escape at all.

The code variants mentioned use different accounting schemes for escape, in terms of initial and maximal attained iteration counts. I actually like a scheme where the count starts at 1 but for which unattained escape is recorded as value 0, since that information can be used to say mask out pixel intensity.

It's interesting to me that when using (plot3d or) plots[densityplot] and its restricttoranges=true option the purely real axis is more closely utilized if the grid has an odd number of points in the imaginary (y) direction. This effect can show more clearly if the maximal iteration limit is not high (30, say).

The output=layer1 option to Maple 18's Fractals:-EscapeTime:-Mandelbrot produces output where nonescape is recorded as value 0. This leads to a similar hue shading as where nonescape is recorded using the maximal iteration value (or 1 higher), because of the mentioned hue shading wrap-around with both extremes showing as red.

There is a mistake in that Library routine, where module local procedure Fractals:-EscapeTime:-sourceMandelbrot fails to check and record a possibly early escape at iteration 1. It starts checking at 2.

One thing that routine does have is additional easy checks against the input point being in either the largest or secondary bulbs of nonescaping points. So here is a density plot variant that does those checks and runs a bit faster for the interesting ranges like in the discussed code variants. I haven't tried to optimize the code's subexpressions, but for this code below on my machine the cross-over for speed is about 75 as the maximal iteration limit. More iterations than that and the easy checks pay off, while for fewer iterations their overhead is a penalty.

restart:

Mandelfun2 := proc(a, b, N)
local z1, z2, z1s, z2s, i, t, temp;
 z1 := a; z2 := b;
 z1s := z1^2; z2s := z2^2;
 if (z1+1.0)^2+z2s < .625e-1 then
    return N+1;
 end if;
 t := (z1-.25)^2+z2s;
 if (t+z1-.25)*t < .25*z2s then
    return N+1;
 end if;
 for i to N while z1s+z2s < 4 do
   temp := z1s - z2s + a;
   z2 := 2*z1*z2 + b; z1 := temp;
   z1s := z1^2; z2s := z2^2;
 end do;
 i;
end proc:

Mandelfun := proc(x,y,N)
local n,c,z;
     c:= x+I*y;
     z:= 0;
     for n to N while abs(z) < 2 do
          z:= z^2+c
     end do;
     n
end proc:

n := 250:

CodeTools:-Usage(
   plots:-densityplot(
        'Mandelfun2'(x,y,75), x=-2.0..0.7, y=-1.35..1.35, grid=[n,n+1],
        colorstyle= HUE, restricttoranges=true,
        axes= boxed, scaling= constrained, 
        style= patchnogrid,
        labels= [Re,Im]) ):
memory used=1.39MiB, alloc change=0 bytes, cpu time=499.00ms, real time=497.00ms, gc time=0ns

CodeTools:-Usage(
   plots:-densityplot(
        'Mandelfun'(x,y,75), x=-2.0..0.7, y=-1.35..1.35, grid=[n,n+1],
        colorstyle= HUE, restricttoranges=true,
        axes= boxed, scaling= constrained, 
        style= patchnogrid,
        labels= [Re,Im]) ):
memory used=1.02MiB, alloc change=0 bytes, cpu time=500.00ms, real time=499.00ms, gc time=0ns

In Maple 18 the revised plots:-surfdata command can be used to easily turn a single layer black & white image Array into a hue shaded density plot. Here's a simple procedure that calls the Fractals:-EscapeTime:-Mandelbrot command as then follows that up with such a conversion. The timings seemed to improve in subsequent runs in the same GUI session, even after restart.

restart:

# Initialize the Compiler and ModuleLoad,
# incurring a one-time hit after each restart which
# may reasonably be excluded from timing results
# if we expect to convert more than once.
with(Fractals:-EscapeTime):

MandelDens:=proc(N::posint,
                 cll::complex(numeric),
                 cur::complex(numeric),
                 {cutoff::posint:=4,
                  iterationlimit::posint:=100})
   local M;
   M := Fractals:-EscapeTime:-Mandelbrot(
          N, cll, cur,
          ':-cutoff'=cutoff,
          ':-iterationlimit'=iterationlimit,
          ':-output'=':-layer1');
   rtable_options(M,':-subtype'=Matrix);
   LinearAlgebra:-Transpose(M,':-inplace');
   rtable_options(M,':-subtype'=Array);
   plots:-surfdata(':-color'=':-COLOR'(':-HUE',M),
                   Re(cll)..Re(cur),Im(cll)..Im(cur),
                   ':-dimension'=2,
                   ':-style'=':-patchnogrid',
                   ':-scaling'=':-constrained',
                   ':-axes'=':-box', _rest);
end proc:

n := 251:

CodeTools:-Usage(
  MandelDens(n, -2.0-1.35*I, 0.7+1.35*I,
             iterationlimit=75, labels=[Re,Im]) ):
memory used=3.24MiB, alloc change=17.73MiB, cpu time=63.00ms, real time=30.00ms, gc time=0ns

A more general purpose image to densityplot conversion utility would be useful. It might handle HUE, HSV, or RGB shadings. I'm thinking of Classic GUI users who cannot use ImageTools:-Embed, and for whom densityplots behave better. And another utility for densityplot to image conversion would also be useful. There I'm thinking of Standard GUI users, where images may stress the interface less. Perhaps more on that later.

@Christopher2222 It does not work for me on Maple 18.01 with symbol=point on 64bit Maple on Windows 7 Pro (ATI Radeon HD 4300/4500).

Increasing the symbolsize or expanding the `view` does not remedy it.

But it works fine without explicitly specifying symbol=point in the command and so without the accompanying substructure in the PLOT structure (ie. default symbol as per GUI, and default symbolsize).

@Carl Love Could `curves,rest` be split out simply just with selectremove of that specfunc type (using a list as originally, not a Vector)?

@Carl Love I suspect that one of the key aspects is that you and I switched it from Re(Int(...)) to Int(Re(...)).

And that can be made faster (even still allowing adaptive plotting). The first version I gave took about 8 seconds. The one below takes under 2 seconds on my Maple 18 (Intel i7, 64bit Windows) even without relaxing the requested accuracy of the numeric integration.

restart:

FF := Q-1+(1/5)*K*dp^3*h^5+(1/3)*dp*h^3+h+h1*h:
DDP:=[solve(FF,dp)]:
h:=1+phi*cos(2*Pi*x):
h1:=2*Pi*alpha*beta*phi*cos(2*Pi*x):
beta:=1:alpha:=0: phi:=1/2:
dpdx:=evalf[15](radnormal(DDP[1])):

#dpp:=q->Int(unapply(Re(eval(dpdx,[K=-0.1,Q=q])),x),0..1,epsilon=1e-5):
dpp:=q->Int(unapply(Re(eval(dpdx,[K=-0.1,Q=q])),x),0..1):

plot(dpp,-1..1,axes=box,color=[blue]);

And a little more can be shaved off by instead using the commented line above which imposes a laxer tolerance on the numeric integration, getting it down to about 1.5 seconds on my machine.

This works for me in both Maple 17.02 and 18.00.

restart:                                                 

M := 2: alpha := 1/2:                                    

U := sum(sum(binomial(n-1, i)*x^(n-i-alpha)*(-a*n)^i*c[n]
             *GAMMA(n-i+1)/GAMMA(n-i-alpha+1),           
             i = 0 .. n-1), n = ceil(alpha) .. M):       

radnormal(convert(U,elementary));                        

                                                        1/2
                        (6 a c[2] - 4 x c[2] - 3 c[1]) x
                   -2/3 -----------------------------------
                                         1/2
                                       Pi

@itsme This might produce the same edited version of Grid:-Map (in 17.02 or 18.01 perhaps) that you suggested, using anames('user').

HackGridMap :=
   FromInert(subsop([5,1,1,2,1]
             =subs(_Inert_LOCAL(1)=_Inert_LOCAL(7),
                   op([5,1],
                      ToInert(proc() local var;
                                var:={anames('user'),
                                      '`grid/mapcmd`'};
                              end proc))),
             ToInert(eval(Grid:-Map)))):

@Michael The second argument to `subs`, which you have as m mod(4) , is evaluated to m even before `subs` gets to do the substitution.

@sarra Your line,

   alpha:=(n,m)->beta[1];

makes the same mistake made here (and in at least once other repost of your this issue that you made). It should be,

   alpha:=unapply(beta[1],[n,m]);

instead.

@Alejandro Jakubi I get simplify-with-siderels as you used it to be about twice as fast as applyrule as I gave it, for many repetitions of the given example if using `forget` between each iteration.

It may also be worth noting that the purely syntactic replacement done by `applyrule` won't automatically match anything in, say, (d+x)*(d-x) to be d^2.

I think that whether one wants just syntactic substitution, or targeted simplification, or to handle just the numerator, etc, can be key to these kinds of tasks even though the user may not always realize what all the implications are. How should the various related tasks be documented by example, I wonder?

@Carl Love The time() command reports the total for all threads, and that can greatly exceed the wall-clock time reported by time[real]().

The output that you ought to have seen is (actually),

> n:=5;

                                    n := 5

> error "bad argument: %1", n; 

Error, bad argument: 5

Note the `5` appearing in the emitted error message. The %1 is a placeholder for the first parameter used in the formatting of the error message.

Here's another example,

> n:=5:

> error "%-1 argument does not match %-3 argument: %2 <> x", 13, n, 7;

Error, 13th argument does not match 7th argument: 5 <> x

acer

@Joe Riel If c23 is not zero then isn't Matrix([[0,c12,c13], [0,1,c23], [0,0,1]]) defective?