acer

edit...

Answered: acer 32333

September 11 2020

1 2

You are encountering an (old, eg Maple 11) evaluation issue in animate.

(The problem was fixed many years/releases ago.)

There are several possible workarounds, for your example in your version.

In your Maple 11, replace the call,

animate(pointplot3d, [ [b(t),c(t),d(t)], symbol=box, color=blue],t=15..100,
background = odeplot(sol, [x(t),y(t),z(t)],t=15..100,numpoints=7000),
frames=20);

with,

animate(t->pointplot3d([b(t),c(t),d(t)], symbol=box, color=blue),[t],t=15..100,
background = odeplot(sol, [x(t),y(t),z(t)],t=15..100,numpoints=7000),
frames=20);

units...

Answered: acer 32333

September 09 2020

2 1

I'm not sure what you want to do about "c". Do you want it as some new unit, and if so could you specify completely? I'm not sure I understand what you want with "c=h=1".

>

Units:-UseUnit(MeV); Units:-UseUnit(MeV*s)

>

pmin := 0.3e9*h*Unit('m'/'s')/((2*ro)*c)

>

>	By the way...

>	By the way...

>

Pmin := Unit('h_bar')*GetValue(cc)*GetUnit(cc)/((2*ro)*c)

>

Download unitsMeVs.mw

Numeric Formatting of NULL assignment (e...

Answered: acer 32333

September 09 2020

1 1

The problem seems to lie in how the GUI is attempting to typeset the echoed assignment, where the result from fsolve is NULL (no real solution found).

In particular the problem seems to be that the problematic paragraph in question has had its "Numeric Formatting" set to "fixed" with 3 decimal places. It seems that when the result of the computation is NULL this goes wrong and it mishandles trying to typeset that in the specified numeric format (which is a GUI bug).

I had difficulty fixing it by, say, changing the d_sample value to 4.5 so the an numeric result was attained, toggling off Numeric Formatting (I think successfully), re-executing OK, then changing the value of d_sample back to 4.6 so that no real root would be found. The problem remained. But I was able to reproduce the problem (see at end) in a way that convinces me that Numeric Formatting of the assignment of NULL is the cause.

One way to fix the example is to place a full colon as terminator of the problematic 2D Input line. Then d_cam_solved gets properly assigned NULL, but the problematic echoing of the assignment is avoided.

Another way to avoid the problem is to wrap the call to fsolve in square brackets, eg, [fsolve(...)] so that in the problematic case the result is the empty list [].

Another way to fix the example is to rewrite the whole line (manually), but without setting its Numeric Formatting. But not cut&paste of the input line, as that can inherit the prior Numeric Formatting. See attached.

Download Ray_transfer_matrix_ac.mw

I can reproduce the problem elsewhere. I create a Document with the 2D Input,
foo:=1.234567
and then execute, and then set the Numeric Formatting to say "fixed" with 3 decimal places. Then I change the 2D Input to instead be,
foo:=NULL
and re-execute, and this problem appears. I will submit a bug report.

The problem also occurs if the input is in 1D plaintext Maple Notation, if the 2D Output Numeric formatting is so-adjusted, etc.

like this...

Answered: acer 32333

September 08 2020

1 5

plots:-animate(plot3d,[[Re,Im](u^4*BesselJ(alpha+4,u)*BesselJ(alpha,v)),
                       u=-12..12, v=-12..12, grid=[51,51],
                       color=[red,blue]],
               alpha=-1..4, frames=50);

Below I make no extra attempt at gaining efficiency (eg. by memoization).

plots:-animate(plot3d,[abs(u^4*BesselJ(alpha+4,u)*BesselJ(alpha,v)),
                       u=-12..12, v=-12..12, grid=[51,51],
                       color=[argument(u^4*BesselJ(alpha+4,u)*BesselJ(alpha,v)),
                              1,1,colortype=HSV]],
               alpha=-1..4, frames=50);

various ways...

Answered: acer 32333

September 08 2020

2 0

You can use plot(...,style=point) or you can use Statistics:-ScatterPlot.

You can also add options axes=box and a color, and have plot(...,style=point) produce something similar to the ScatterPlot result here.

Sieving out the non-numeric data has an added benefit that the plot determines the horizontal range of the valid data more nicely.

>

restart;

>

N := 200:

>	V1 := LinearAlgebra:-RandomVector(N,generator=-1.0..1.0):

>	# # V2 contains some Float(undefined) entries. # V2 := map(x->RealDomain:-sqrt(0.4-x)-RealDomain:-sqrt(0.3+x),V1):

>	M:=<V1\|V2>:

>	M[1..6,..];

Matrix(6, 2, {(1, 1) = .589662833766906092, (1, 2) = Float(undefined), (2, 1) = .635255416644524118, (2, 2) = Float(undefined), (3, 1) = 0.215431283442193422e-1, (3, 2) = 0.481407545e-1, (4, 1) = 0.170173107622539899e-1, (4, 2) = 0.558130487e-1, (5, 1) = -.387301055966885244, (5, 2) = Float(undefined), (6, 1) = -.106432501140387492, (6, 2) = .2716776453})

>	plot(M,style=point,size=[300,300]);

>	A:=`<,>`(select(type,[seq(M[i,..],i=1..N)],'Vector'(numeric))[]):

>	A[1..6,..];

Matrix(6, 2, {(1, 1) = 0.215431283442193422e-1, (1, 2) = 0.481407545e-1, (2, 1) = 0.170173107622539899e-1, (2, 2) = 0.558130487e-1, (3, 1) = -.106432501140387492, (3, 2) = .2716776453, (4, 1) = -.128282822838161836, (4, 2) = .3124429564, (5, 1) = -0.264167351936552830e-1, (5, 2) = .1299540470, (6, 1) = .251237121459380708, (6, 2) = -.3567555365})

>	# # If you want the horizontal domain to be restricted # automatically, then you could select only the numeric rows. # plot(A,style=point,size=[300,300]);

>	plot(A,style=point,color="SteelBlue",axes=box,size=[300,300]);

>	# # ScatterPlot is here (mostly) acting like the plot(..,style=point). # Statistics:-ScatterPlot(A,size=[300,300]);

>	# # ScatterPlot can handle Float(undefined). # # It also doesn't determine the nicer horizontal range automatically. # Statistics:-ScatterPlot(<V1\|V2>,size=[300,300]);

>	# # Or, using a list of lists (of the pairs of values). # L1:=[seq([V1[i],V2[i]],i=1..N)]: #L1[1..6,..]; plot(L1,style=point,size=[300,300]);

>	L2:=select(hastype,L1,list(realcons)): #L2[1..6,..]; plot(L2,style=point,size=[300,300]);

>	# # Or, generate the shorter list of lists directly, # without wasting space on the unwanted points. # L3:=[seq(`if`(V2[i]::realcons,[V1[i],V2[i]],NULL),i=1..N)]: #L3[1..6,..]; plot(L3,style=point,size=[300,300]);

>	Statistics:-ScatterPlot(L3,size=[300,300]);

>

Download point_plot_undef.mw

two ways...

Answered: acer 32333

September 07 2020

1 1

Are you trying to do something like this?

Download test_cat_ac.mw

Here is another way to get that effect, using the cat command,

Download test_cat_ac2.mw

one way...

Answered: acer 32333

September 07 2020

2 3

I suppose that you might utilize PDEtools:-dcoeff, except that it doesn't show the zero coefficients.

Here is one simple approach (which you could make stronger, adjust, etc).

>

restart;

>	H := proc(ee::algebraic, f::function) local i,n,v; if nops(f)<>1 or (not op(1,f)::name) then error "second argument f is not a function of a name"; else v := op(1,f); end if; n := PDEtools:-difforder(ee,v); < seq( frontend(coeff,[ee,diff(f,v$i)]), i=1..n) >; end proc:

>

>	z := f0(x) + f1(x)diff(y(x),x) + f2(x)diff(y(x),x,x) + f4(x)*diff(y(x),x,x,x,x):

>	H(z, y(x));

Vector(4, {(1) = f1(x), (2) = f2(x), (3) = 0, (4) = f4(x)})

>	z := f0(x) + f2(x)diff(y(x),x,x) + f4(x)diff(y(x),x,x,x,x):

>	H(z, y(x));

Vector(4, {(1) = 0, (2) = f2(x), (3) = 0, (4) = f4(x)})

>

Download simple_dcoeff.mw

one way...

Answered: acer 32333

September 06 2020

1 0

You can write a procedure to handle the task.

This is for the specific form Int(Sum(...),...). If instead you have a form like Int(A+B*Sum(...),...) then you'd have to enhance the procedure, or first massage into the specified form (using IntegrationTools:-Expand, etc).

Adjust and strengthen or weaken as needed.

[edit] I made an alteration to allow additional arguments of the Int call, see third example.

>

restart;

>

F:=e->subsindets(e,And(specfunc(anything,Int),
                       satisfies(u->nops(u)>1 and
                                    type(op(1,u),
                                         specfunc(anything,Sum))
                                    and
                                    type(op(2,u),
                                         {name,name=range}))),
                 u->Sum(Int(op([1,1],u),op(2..,u)),op([1,2..],u))):

>	ee:=Int(Sum(-(Nbt-nt__rev) exp((1/2)(L+sqrt(-4CLR^2+L^2))t/(LRC) -(1/2)(Nbt-nt__rev)^2/(Nb^2sigma__b^2)), n = 0 .. Nb-1), t);

Int(Sum(-(Nb*t-n*t__rev)*exp((1/2)*(L+(-4*C*L*R^2+L^2)^(1/2))*t/(L*R*C)-(1/2)*(Nb*t-n*t__rev)^2/(Nb^2*sigma__b^2)), n = 0 .. Nb-1), t)

>

F(ee);

Sum(Int(-(Nb*t-n*t__rev)*exp((1/2)*(L+(-4*C*L*R^2+L^2)^(1/2))*t/(L*R*C)-(1/2)*(Nb*t-n*t__rev)^2/(Nb^2*sigma__b^2)), t), n = 0 .. Nb-1)

>

>	foo := bar + ee + Int(K(t),t) + Int(Sum(P(t,n),n=1..N),t=a..b): F(foo);

bar+Sum(Int(-(Nb*t-n*t__rev)*exp((1/2)*(L+(-4*C*L*R^2+L^2)^(1/2))*t/(L*R*C)-(1/2)*(Nb*t-n*t__rev)^2/(Nb^2*sigma__b^2)), t), n = 0 .. Nb-1)+Int(K(t), t)+Sum(Int(P(t, n), t = a .. b), n = 1 .. N)

>

>	blah := Int(Sum(K(t,n),n=0..N),t,allsolutions) + Int(Sum(P(t,n),n=1..N),t=a..b,method=_Dexp,epsilon=1e-5): F(blah);

Sum(Int(K(t, n), t, allsolutions), n = 0 .. N)+Sum(Int(P(t, n), t = a .. b, method = _Dexp, epsilon = 0.1e-4), n = 1 .. N)

>

Download IntofSum.mw

another way...

Answered: acer 32333

September 05 2020

0 2

Suppose that L is a list of your lists. Then,

plots:-display(map(plots:-listplot,L));

something...

Answered: acer 32333

September 04 2020

1 0

I didn't look at trying to set up a double integral (of R1 wrt a and q, as a function of N).

With the iterated single integrations the key seemed to be to get the inner integration (of R1 wrt a, as function of N and q) to be as fast as possible. But requesting more than modest accuracy of that takes longer. And the outer integration's possible accuracy depends on the accuracy of the inner (though to what extent I don't know).

You fiddle with this, and try to corroborate the results. There are some choices for the single integrations using _Dexp for a complex integrand and _d01ajc for split real and imaginary parts.

I used Maple 2017 because I suspect that is what you're using.

>

restart;

>	tt := -0.689609e-3; T_c := .242731; mu := .365908; k := 1;

>	R1 := atanh((a^2-mu)/(2T_c))ln((2a^2+2aq+q^2-2mu-(I2)PiN) /(2a^2-2aq+q^2-2mu-(I2)Pi*N))/q-2:

>

R2 := proc(q,N)
        local res;
        Digits:=15;
        if not [q,N]::list(numeric) then
          return 'procname'(q,N);
        end if;
        res:=evalf(Int(unapply(eval(R1,[:-q=q,:-N=N]),a), 0 .. 10000, epsilon=1e-4, method=_Dexp));
        #res:=evalf(Int(unapply(eval(Re(R1),[:-q=q,:-N=N]),a), 0 .. 10000, epsilon=1e-3, method=_d01ajc)
        #           +I*Int(unapply(eval(Im(R1),[:-q=q,:-N=N]),a), 0 .. 10000, epsilon=1e-3, method=_d01ajc));
        if not res::complex(numeric) then
          res:=evalf(Int(unapply(eval(R1,[:-q=q,:-N=N]),a), 0 .. 10000, epsilon=1e-3, method=_Dexp));
        end if;
        if not res::complex(numeric) then
          error q,N,eval(R1,[:-q=q,:-N=N]);
        end if;
        res;
end proc:

>	R3 := qln((-q^2-k^2+mu+I(2NPiT_c-(2m+1)PiT_c)+kq) /(-q^2-k^2+mu+I(2NPiT_c-(2m+1)PiT_c)-kq)) /(k(tt+R2(q,N))):

>	R4 := proc(q,a,b) if not [q,a,b]::list(numeric) then return 'procname'(q,a,b); end if; add(eval(R3,[:-q=q]), N = a .. b); end proc:

>

m := 1;

>	#plot(eval([Re,Im](op([2,2,1],R2)),[q=1,N=-1])); forget(evalf); forget(`evalf/int`); CodeTools:-Usage( evalf(eval(R3,[q=1,N=1])) ); # d01ajc .2077805664+0.5380483364e-1*I

memory used=22.86MiB, alloc change=34.00MiB, cpu time=151.00ms, real time=151.00ms, gc time=16.38ms

>	forget(evalf); forget(`evalf/int`); CodeTools:-Usage( evalf(eval(R4(q,-100,100),[q=1])) );

memory used=2.70GiB, alloc change=36.00MiB, cpu time=17.10s, real time=16.83s, gc time=969.43ms

>	forget(evalf); forget(`evalf/int`); CodeTools:-Usage( evalf(eval(R4(q,-100,100),[q=0.1])) );

memory used=2.80GiB, alloc change=-2.00MiB, cpu time=17.74s, real time=17.21s, gc time=1.19s

>	CodeTools:-Usage( evalf(Int(q->R4(q,-100,100), 100 .. 10000, epsilon=1e-2, method=_Dexp)) );

memory used=19.48GiB, alloc change=0 bytes, cpu time=2.07m, real time=2.00m, gc time=8.79s

>	CodeTools:-Usage(plot([Re,Im](R4(q,-100,100)), q=10 .. 10000, adaptive=false, numpoints=23, size=[300,200]));

memory used=34.93GiB, alloc change=0 bytes, cpu time=3.77m, real time=3.64m, gc time=17.08s

>	CodeTools:-Usage( evalf(Int(q->Re(R4(q,-100,100)), 100 .. 10000, epsilon=1e-2, method=_d01ajc)) );

memory used=96.88GiB, alloc change=0 bytes, cpu time=10.50m, real time=10.03m, gc time=63.43s

>	CodeTools:-Usage( evalf(Int(q->Im(R4(q,-100,100)), 100 .. 10000, epsilon=1e-2, method=_d01ajc)) );

memory used=137.28GiB, alloc change=0 bytes, cpu time=15.39m, real time=14.71m, gc time=92.31s

>

Download hard_evalf_int.mw

sure...

Answered: acer 32333

September 03 2020

0 0

Threads:-Map(ff,Array([1,2,3]));

            [ff(1), ff(2), ff(3)]

various...

Answered: acer 32333

September 03 2020

3 2

For the simple task of colorizing a string, you could start with this:

Download string_color.mw

More complicated things, like handling 2D Math, are possible. But I'd have to dig because even though it's been asked several times this site's Search doesn't index on code in verbatim/pre tags.

I also found this one (an edit from here). Adjust to taste. I've done similar things with both more or less flexibility, and some that also supported font size and bold/italic, etc. Let me know if you have specific needs.

>

restart;

>

P := proc()
  local x, lastSpace := "    ", mainColor := black,
        lastColor := black, i, q := [], new;
  uses T=Typesetting, C=ColorTools;
  for i from 1 to _npassed do
    x := _passed[i];
    if type(x, string) then
      # Parse strings and handle special modifiers (@color, " #n", etc)
      if x[1] = " " and x[2] = "#" then
        x := cat(" "$parse(x[3..-1])); lastSpace := x;
      end if;
      if numelems(x) = 0 then
        x := lastSpace;
      end if;
      x := x[1..-1-SearchText("@", StringTools:-Reverse(x))];
      lastColor := `if`(numelems(x) < numelems(_passed[i]),
                        _passed[i][numelems(x)+2..-1], lastColor);
      if type(parse(lastColor), integer) then
        lastColor := C:-GetColorNames(C:-GetPalette("HTML"))[parse(lastColor)];
      end:
      if numelems(x) = 0 then
        mainColor := lastColor; continue;
      end if;
      q := [op(q), T:-mtext(x, ':-color'=lastColor)];
    else
      new := T:-Typeset(T:-EV(x));
      new := subsindets(new, identical(:-mathcolor)=anything,
                        ()->':-mathcolor'=lastColor);
      new := subsindets(new, specfunc({T:-mn,T:-mi,T:-mo,T:-ms,T:-mfrac,
                                       T:-mover,T:-msub,T:-msup,T:-msubsup,
                                       T:-mfenced,T:-msqrt,T:-mroot}),
                        u->op(0,u)(op(u),':-mathcolor'=lastColor));
      #new := subsindets(new, specfunc({T:-mo,T:-mi,T:-mn}),
      #                  u->op(0,u)(op(u),':-fontweight'="bold"));
      q := [op(q), new];
    end if;
  end do;
  Typesetting:-mrow(op(q));
end proc:

>	P("test@Red", " #4", [1,3,4], "\n@Purple", (3/4+surd(x,5))^2, "", conjugate(u), "\n@Green", sqrt(y), "", Int(f(x)*sin(x),x=a..b));

$""test"" "[1,3,4]"\n"(3/4+x^(1/5))^2" "u"\n"sqrt(y)" "(∫)[a]^bf(x) sin(x) &DifferentialD;x"$

>

Download Color_Typesetting.mw

And here's another, simpler example. It's not in a procedure, but illustrates the basic idea.

>

restart;

>	expr := sqrt(Sum(sin(nPix),n=1..N));

(Sum(sin(n*Pi*x), n = 1 .. N))^(1/2)

>	foo:=Typesetting:-Typeset(Typesetting:-EV(expr)):

>	lprint(foo);

Typesetting:-msqrt(Typesetting:-mrow(Typesetting:-munderover(Typesetting:-mstyle(Typesetting:-mo("∑", Typesetting:-msemantics = "inert"), mathcolor = "#909090"), Typesetting:-mrow(Typesetting:-mi("n"), Typesetting:-mo("="), Typesetting:-mn("1")), Typesetting:-mi("N")), Typesetting:-mo("⁡"), Typesetting:-mrow(Typesetting:-mi("sin", fontstyle = "normal"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi("n"), Typesetting:-mo("⁢"), Typesetting:-mi("π"), Typesetting:-mo("⁢"), Typesetting:-mi("x"))))))

>	map[2](op,0,indets(foo,function));

>	stuff1 := fontfamily="Helvetica", size = 16, fontweight = "bold", mathcolor = "purple":

>

subsindets(foo, 'specfunc({Typesetting:-msqrt,Typesetting:-mrow,
                           Typesetting:-mi,Typesetting:-mo,Typesetting:-mn,
                           Typesetting:-mfenced,Typesetting:-munderover})',
           u->op(0,u)(op(remove(type,[op(u)],
                         identical(fontfamily,size,
                                   fontweight,mathcolor)=anything)),stuff1));

sqrt(Sum(sin(n*Pi*x), n = 1 .. N))

>

# including mstyle, to replace gray of inert Sum

subsindets(foo, 'specfunc({Typesetting:-msqrt,Typesetting:-mrow,
                           Typesetting:-mi,Typesetting:-mo,Typesetting:-mn,
                           Typesetting:-mfenced,Typesetting:-munderover,
                           Typesetting:-mstyle})',
           u->op(0,u)(op(remove(type,[op(u)],
                         identical(fontfamily,size,
                                   fontweight,mathcolor)=anything)),stuff1));

sqrt(Sum(sin(n*Pi*x), n = 1 .. N))

>

Download Typesetting_misc.mw

modification...

Answered: acer 32333

September 03 2020

2 1

I suppose that you already know that this computation is better suited to Threads than Grid, since these particular in-place computations are not only "embarassingly paralellizable" but are also thread-safe.

Indeed, the Mandelbrot command does it using Threads (and Compile or evalhf which are faster then option hfloat). (See also here and perhaps here for a related work on a Newton fractal.)

But it may still be an interesting example for data passing under the Grid package.

I made a couple of adjustments and got the Grid version to work with a speedup factor of about 2 (3.4sec versus 7.3sec at size N=300) on my 4-node Linux box. (I retained the option hfloat and didn't try to Compile either, as my current desktop has gcc issues...)

By the way, the parsing error about break that you saw seems to be due to the option hfloat.

Here is the revision:

>

restart;

>

Mandelbrot := module()
    local MandelLoop,
            MandelGrid,
            ModuleApply;

    MandelLoop := proc( X, Y, imageArray, i_low, i_high, j_low, j_high, iter, bailout )
        local str, flag, i, j, Xc, Yc, Xtemp, Ytemp, Xold, Yold, k, t;
        option hfloat;
        str:=time[real]();
        for i from i_low to i_high do
           for j from j_low to j_high do
               Xtemp := X[i];
               Ytemp := Y[j];
               Xc := Xtemp;
               Yc := Ytemp;
               k := 0;
               flag := true:
               while k < iter and flag do
                   Xold := Xtemp;
                   Yold := Ytemp;
                   Xtemp := Xold^2-Yold^2+Xc;
                   Ytemp := 2*Xold*Yold+Yc;
                   t := Xtemp^2+Ytemp^2;
                   if t >= bailout then
                       imageArray[i, j, 1] := k - ln( ln( t ) )/ln(2.);
                       imageArray[i, j, 2] := imageArray[i, j, 1];
                       imageArray[i, j, 3] := imageArray[i, j, 1];
                       #break;
                       flag := false:
                  end if;
                  k := k+1;
               end do
           end do;
        end do;
        NULL:
    end proc:

    MandelGrid := proc( X, Y, iter, bailout )
        local i, n, step, imageData, start, endp;

        n := Grid:-NumNodes(); #n:=1;
        i := Grid:-MyNode();
        step := floor( numelems( X )/n );

        if ( i = 0 ) then
            start := 1;
            endp := step;
        elif ( i = n-1 ) then
            start := step*(n-1)+1;
            endp := numelems(X);
        else
            start := step*i+1;
            endp := step*(i+1);
        end;
        imageData := Array( start..endp, 1..numelems(Y), 1..3, datatype=float[8] );
        :-MandelLoop( X, Y, imageData, start, endp, 1, numelems(Y), iter, bailout );

        if ( i > 0 ) then
            Grid:-Send(0,imageData);
        else
            [ imageData, seq( Grid:-Receive(i), i = 1..n-1 ) ];
        end;
    end proc:

    ModuleApply := proc ( ptsX, ptsY, iter, X1, X2, Y1, Y2, bailout )

        local X, Y, imageData, ret, i, l, u:

        X := Vector(ptsX, i->X1+(X2-X1)*(i-1)/(ptsX-1) , datatype = float[8]);
        Y := Vector(ptsY, i->Y1+(Y2-Y1)*(i-1)/(ptsY-1) , datatype = float[8]);
        ret := Grid:-Launch( MandelGrid, X, Y, iter, bailout,
                             imports=[ ':-MandelLoop'=eval(MandelLoop) ] );
        imageData := Array( 1..ptsX, 1..ptsY, 1..3, datatype=float[8] );

        for i in ret
        do
            l := lowerbound( i );
            u := upperbound( i );
            imageData[l[1]..u[1], l[2]..u[2], 1..3] := i;
        end;

        imageData;
    end proc:
end:

>

>	Grid:-NumNodes();

>

N := 300:

>	s := time[real](): points := Mandelbrot( N, N, 100, -2.0, .7, -1.35, 1.35, 10.0 ): time[real]()-s;

>	op(2,points);

>	plots:-display(ImageTools:-Preview(points), size=[300,300], axes=none);

>

Download MandelLoop_try0.mw

And here is the serial version.

>

restart;

>

Mandelbrot := module()
    local MandelLoop,
            ModuleApply;

    MandelLoop := proc( X, Y, imageArray, i_low, i_high, j_low, j_high, iter, bailout )
        local flag, i, j, Xc, Yc, Xtemp, Ytemp, Xold, Yold, k, t;
        option hfloat;

        for i from i_low to i_high do
           for j from j_low to j_high do
               Xtemp := X[i];
               Ytemp := Y[j];
               Xc := Xtemp;
               Yc := Ytemp;
               k := 0;
               flag := true;
               while k < iter and flag do
                   Xold := Xtemp;
                   Yold := Ytemp;
                   Xtemp := Xold^2-Yold^2+Xc;
                   Ytemp := 2*Xold*Yold+Yc;
                   t := Xtemp^2+Ytemp^2;
                   if Xtemp^2+Ytemp^2 >= bailout then
                        imageArray[i, j, 1] := k - ln( ln( t ) )/ln(2.);
                        imageArray[i, j, 2] := imageArray[i, j, 1];
                        imageArray[i, j, 3] := imageArray[i, j, 1];
                        #break;
                        flag := false;
                  end if;
                  k := k+1;
               end do
           end do;
        end do;
    end proc:

    ModuleApply := proc ( ptsY, ptsX, iter, X1, X2, Y1, Y2, bailout )
        local X, Y, imageArray, i:

        X := Vector(ptsX, i->X1+(X2-X1)*(i-1)/(ptsX-1) , datatype = float[8]);
        Y := Vector(ptsY, i->Y1+(Y2-Y1)*(i-1)/(ptsY-1) , datatype = float[8]);
        imageArray := Array(1 .. ptsY, 1 .. ptsX, 1 .. 3, datatype = float[8]);

        MandelLoop( X, Y, imageArray, 1, ptsX, 1, ptsY, iter, bailout );

        return imageArray;
    end proc:
end:

>	N := 300: s := time[real](): points := Mandelbrot( N, N, 100, -2.0, .7, -1.35, 1.35, 10.0 ): time[real]()-s;

>	op(2,points);

>	plots:-display(ImageTools:-Preview(points), size=[300,300], axes=none);

>

Download MandelLoop_nogrid.mw

exact formula correction...

Answered: acer 32333

September 03 2020

1 1

In your code you used,

eval(exp(2-c*h)+(2*c*h)):

as the formula for the exact values. But shouldn't it instead be,

eval(exp(2-(2.0+c*h))+(2*(2.0+c*h))):

where the 2.0 is the starting value (originally assighned to inx, but that gets reassigned).

By the way, you can check your forward Euler scheme's results against dsolve,numeric with its classical[foreuler] method. See below.

restart:

e1:=y[n+1]=y[n]+h*f(n):

inx:=2:

    h               Num.y                Ex.y             Error y
0.10        5.1000000000        5.1048374180             0.00484
0.20        5.2100000000        5.2187307531             0.00873
0.30        5.3290000000        5.3408182207              0.0118
0.40        5.4561000000        5.4703200460              0.0142
0.50        5.5904900000        5.6065306597               0.016
0.60        5.7314410000        5.7488116361              0.0174
0.70        5.8782969000        5.8965853038              0.0183
0.80        6.0304672100        6.0493289641              0.0189
0.90        6.1874204890        6.2065696597              0.0191
1.00        6.3486784401        6.3678794412              0.0192

$seq(printf("x =%5.2f %a\n", x, dsol(x)[2]), x = 2 .. 3, .1);$

x = 2.00 y(x) = 5.

x = 2.10 y(x) = 5.1
x = 2.20 y(x) = 5.21
x = 2.30 y(x) = 5.329
x = 2.40 y(x) = 5.4561
x = 2.50 y(x) = 5.59049
x = 2.60 y(x) = 5.731441
x = 2.70 y(x) = 5.8782969
x = 2.80 y(x) = 6.03046721
x = 2.90 y(x) = 6.187420489
x = 3.00 y(x) = 6.3486784401

Download EULER_Correction_ac2.mw

gamma...

Answered: acer 32333

September 02 2020

1 0

It seems as if you want to utilize the value gamma=0.1 , but haven't. Using gamma=0.1 will produce a different curve than what you had originally (possibly mistakenly).

You could declare local gamma and assign to it (or use a different name throughout).

Also, here are some ideas for the plotting. You can retain or remove whichever features you prefer. It may give you some ideas.

>

restart;

>	with(plots):

>	alpha := 0.1: beta := 0.1: mu := 0.5: u := 0.5: v := 1: local gamma := 0.1;

>	Eq := -2sigmaalphabetamuu - sigmaalphabeta^2uv - 2.0mualphabeta^2uv - gamma^2k^2alphabetau + betasigma^2v + sigma^3 + 2musigma^2 - alphabetasigma^2u + sigmabeta^2uv + gamma^2k^2sigma + 2.0mubetasigmav + 2.0mubeta^2uv;

>	K := [seq(2nPi/10,n=0..15)]:

>	S1 := map(kk->seq(`if`(kkk>=-0.1 and kkk<=0.1,[kk,kkk],NULL), kkk=[fsolve(eval(Eq,k=kk))]), K):

>

display(
plot(S1,color=blue,style=point,symbol=solidcircle,symbolsize=10),
seq(textplot([S1[i][],S1[i][1]],'align'={'right','below'},
              font=["Times",8]),
     i=1..nops(S1)),
implicitplot(Eq, k = 0 .. 10, sigma = -0.1 .. 0.1)
, axes=boxed
, axis[1]=[gridlines=S1[..,1]
            , tickmarks=map(u->u=u,S1[..,1])
           ]
, axesfont=["Times",8]
);

>	display( plot(S1,color=blue,style=point,symbol=solidcircle,symbolsize=10), seq(textplot([S1[i][],S1[i][1]],'align'={'right','below'}, font=["Times",8]), i=1..nops(S1)), implicitplot(Eq, k = 0 .. 10, sigma = -0.1 .. 0.1) , axes=boxed );

>	display( plot(S1,color=blue,style=point,symbol=solidcircle,symbolsize=10), seq(textplot([S1[i][],S1[i][1]10/(2Pi)],'align'={'right','below'}, font=["Times",8]), i=1..nops(S1)), implicitplot(Eq, k = 0 .. 10, sigma = -0.1 .. 0.1) , axes=boxed );

>

S2 := map(kk->seq([kk,kkk],kkk=[fsolve(eval(Eq,k=kk))]),K):
display(
plot(S2,color=blue,style=point,symbol=solidcircle,symbolsize=10),
seq(textplot([S2[i][],S2[i][1]*10/(2*Pi)],'align'={'right','below'},
font=["Times",8]),
i=1..nops(S2)),
implicitplot(Eq, k = 0 .. 10, sigma = -0.1 .. 0.1)
#, view=[min(S1[..,1]) .. 10, -0.1 .. max(S2[..,2])]
, axes=boxed
);

>

Download implicitplot_points.mw

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

32333 Reputation

29 Badges

Social Networks and Content at Maplesoft.com

MaplePrimes Activity

These are answers submitted by acer

edit...

units...

Numeric Formatting of NULL assignment (e...

like this...

various ways...

two ways...

one way...

one way...

another way...

something...

sure...

various...

modification...

exact formula correction...

gamma...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again