Advanced Search

acer

32385 Reputation

29 Badges

19 years, 339 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

Contact acer

MaplePrimes Activity


These are answers submitted by acer

something...

acer 32385
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;

-0.689609e-3

.242731

.365908

1

> 

R1 := a*tanh((a^2-mu)/(2*T_c))*ln((2*a^2+2*a*q+q^2-2*mu-(I*2)*Pi*N)
                                   /(2*a^2-2*a*q+q^2-2*mu-(I*2)*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 := q*ln((-q^2-k^2+mu+I*(2*N*Pi*T_c-(2*m+1)*Pi*T_c)+k*q)
           /(-q^2-k^2+mu+I*(2*N*Pi*T_c-(2*m+1)*Pi*T_c)-k*q))
      /(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;

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

.2077817166+0.5380426505e-1*I

> 

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

.2679032606-.3534769797*I

> 

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

.1676174888-.5972959586*I

> 

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

1241.115378-0.2936769042e-1*I

> 

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

1241.120520

> 

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

-0.2936765300e-1

> 

 

Download hard_evalf_int.mw

 

sure...

acer 32385
September 03 2020
0 0 
Threads:-Map(ff,Array([1,2,3]));

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

various...

acer 32385
September 03 2020
3 2

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

> 

G:=proc(s::string,c::string)
  uses ColorTools;
  nprintf("#mn(%a,mathcolor=%a);",s,
          RGB24ToHex(RGBToRGB24([Color(c)[]])));
end proc:
> 

G("Hello, my name is Bob.", "DodgerBlue");

`#mn("Hello, my name is Bob.",mathcolor="#1E90FF");`

> 

G("Hello, my name is Bob.", "xkcd blue");

`#mn("Hello, my name is Bob.",mathcolor="#0343DF");`

> 

G("Hello, my name is Bob.", "Blue");

`#mn("Hello, my name is Bob.",mathcolor="#0000FF");`

> 

G("Hello, my name is Bob.", "SlateBlue");

`#mn("Hello, my name is Bob.",mathcolor="#6A5ACD");`

> 

 

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)" "(&int;)[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(n*Pi*x),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("&sum;", Typesetting:-msemantics = "inert"), mathcolor = "#909090"), Typesetting:-mrow(Typesetting:-mi("n"), Typesetting:-mo("&equals;"), Typesetting:-mn("1")), Typesetting:-mi("N")), Typesetting:-mo("&ApplyFunction;"), Typesetting:-mrow(Typesetting:-mi("sin", fontstyle = "normal"), Typesetting:-mo("&ApplyFunction;"), Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi("n"), Typesetting:-mo("&InvisibleTimes;"), Typesetting:-mi("&pi;"), Typesetting:-mo("&InvisibleTimes;"), Typesetting:-mi("x"))))))

> 

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

{Typesetting:-mfenced, Typesetting:-mi, Typesetting:-mn, Typesetting:-mo, Typesetting:-mrow, Typesetting:-msqrt, Typesetting:-mstyle, Typesetting:-munderover}

> 

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...

acer 32385
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();

4

> 

N := 300:
> 

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

3.369

> 

op(2,points);

1 .. 300, 1 .. 300, 1 .. 3

> 

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;

7.302

> 

op(2,points);

1 .. 300, 1 .. 300, 1 .. 3

> 

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

> 

 

Download MandelLoop_nogrid.mw

exact formula correction...

acer 32385
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

``

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(evalf[10](eval(2*x+exp(-x)/exp(-2), x = j)), j = 2.0 .. 3.0, .1);

5.000000000, 5.104837419, 5.218730754, 5.340818221, 5.470320046, 5.606530660, 5.748811636, 5.896585304, 6.049328964, 6.206569660, 6.367879441

dsol := dsolve({diff(y(x), x) = 2*(1+x)-y(x), y(2) = 5}, y(x), numeric, method = classical[foreuler], stepsize = .1):

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...

acer 32385
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;

.1

> 

Eq := -2*sigma*alpha*beta*mu*u - sigma*alpha*beta^2*u*v - 2.0*mu*alpha*beta^2*u*v - gamma^2*k^2*alpha*beta*u + beta*sigma^2*v + sigma^3 + 2*mu*sigma^2 - alpha*beta*sigma^2*u + sigma*beta^2*u*v + gamma^2*k^2*sigma + 2.0*mu*beta*sigma*v + 2.0*mu*beta^2*u*v;

0.995e-1*sigma+0.4500e-2-0.5e-4*k^2+1.095*sigma^2+sigma^3+0.1e-1*k^2*sigma

> 

K := [seq(2*n*Pi/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/(2*Pi)],'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

fclose, or restart...

acer 32385
September 02 2020
1 4

The files may not be written to disk on the OS until you either restart the session, or fclose them, or close the worksheet.

So, you can add this line inside the end of the loop.
   fclose(pout);

I am presuming that you will also uncomment the line that executes the plot. You may also want it to be,
  print(plot(...));
if your looping is further nested, or if all this is done within a procedure.

For example,

for j from 1 to 3 do
pout := cat("C:/Users/Eli/Documents/Animation/", "file", j, ".bmp"):
print(pout):
plotsetup(bmp, plotoutput = pout):
print(plot(...));
fclose(pout);
od:

 

2D Input...

acer 32385
September 02 2020
3 0

You don't need differential equations for this. You just need to know how to differentiate and substitute.

For brevity, there is this one-liner:

eval( D(2/f)(-1), [D(f)(-1)=2, f(-1)=4] );

                -1
                --
                 4

Or, slightly longer,

eval( D(h=2/f)(-1), [D(f)(-1)=2, f(-1)=4] );

            D(h)(-1) = -1/4

In 2D Input that could be done tersely as follows,

"h(x):=(2)/(f(x)):"

eval((D(h))(-1), [eval(diff(f(x), x), x = -1) = 2, f(-1) = 4])

-1/4

``

Download calc1.mw


I would have liked to do it with the following shorter syntax.

"h(x):=(2)/(f(x)):"

eval((D(h))(-1), [(D(f))(-1) = 2, f(-1) = 4])


But unfortunately f'(-1) produces the D form rather than diff form.
So the forms don't match, so the substitution doesn't succeed.

This happens even if we assigned   h:=2/f  instead.

lprint((D(h))(-1))

-2/f(-1)^2*eval(diff(f(x),x),{x = -1})

lprint((D(f))(-1))

D(f)(-1)
``

``

Download calc2.mw

Here are a few more ways, done using 2D Input. (It's a little awkward, since the 2D Input of  f'(-1) is hooked into the D form rather than the eval(diff(...)...) form.

> 

restart
> 

Typesetting:-Settings(typesetprime = true); Typesetting:-Settings(useprime = true)

true

> 

eq := h(x) = 2/f(x)

h(x) = 2/f(x)

> 

neweq := map(diff, eq, x)

diff(h(x), x) = -2*(diff(f(x), x))/f(x)^2

> 

this := eval(neweq, x = -1)

eval(diff(h(x), x), {x = -1}) = -2*(eval(diff(f(x), x), {x = -1}))/f(-1)^2

> 

eval(this, [eval(diff(f(x), x), {x = -1}) = 2, f(-1) = 4])

eval(diff(h(x), x), {x = -1}) = -1/4

> 

restart
> 

Typesetting:-Settings(typesetprime = true); Typesetting:-Settings(useprime = true)

true

> 

eq := h(x) = 2/f(x)

h(x) = 2/f(x)

> 

neweq := map(diff, eq, x)

diff(h(x), x) = -2*(diff(f(x), x))/f(x)^2

> 

eval(-2*(D(f))(-1)/f(-1)^2, [(D(f))(-1) = 2, f(-1) = 4])

-1/4

> 

restart
> 

eq := h = 2/f

h = 2/f

> 

new := (map(D, eq))(-1)

(D(h))(-1) = -2*(D(f))(-1)/f(-1)^2

> 

eval(new, [(D(f))(-1) = 2, f(-1) = 4])

(D(h))(-1) = -1/4

> 

 

Download calcexample.mw

And, using context-menus,

restart

Typesetting:-Settings(typesetprime = true); Typesetting:-Settings(useprime = true)

h(x) = 2/f(x)

h(x) = 2/f(x)

"(->)"

diff(h(x), x) = -2*(diff(f(x), x))/f(x)^2

"(->)"

eval(diff(h(x), x), {x = -1}) = -2*(eval(diff(f(x), x), {x = -1}))/f(-1)^2

eval(convert(%, D), [(D(f))(-1) = 2, f(-1) = 4])

(D(h))(-1) = -1/4

``

Download calcexample2.mw

symbolic...

acer 32385
September 02 2020
0 2

Even in old Maple 17 the Help pages for radsimp states, "The original radsimp has been removed and replaced by a wrapper to simplify(...,radical,symbolic)."

If you don't want to operate under assumptions about the signum of the squared term under the radical then don't use the symbolic option (or radsimp).

> 

restart;
> 

kernelopts(version);

`Maple 17.01, X86 64 LINUX, Jun 25 2013, Build ID 849430`

> 

m := M-M^2/(2*R)+Q^2/(2*R);

M-(1/2)*M^2/R+(1/2)*Q^2/R

> 

temp := subs(M=4*Pi*R^2*sigma,Q=4*Pi*R^2*sigmae,sigmae=beta*sigma,m-sqrt(m^2-Q^2));

4*Pi*R^2*sigma-8*Pi^2*R^3*sigma^2+8*Pi^2*R^3*beta^2*sigma^2-(-16*Pi^2*R^4*beta^2*sigma^2+(8*Pi^2*R^3*beta^2*sigma^2-8*Pi^2*R^3*sigma^2+4*Pi*R^2*sigma)^2)^(1/2)

> 

new := subs(beta=0,temp);

4*Pi*R^2*sigma-8*Pi^2*R^3*sigma^2-((-8*Pi^2*R^3*sigma^2+4*Pi*R^2*sigma)^2)^(1/2)

> 

radsimp(new);

-16*Pi^2*R^3*sigma^2+8*Pi*R^2*sigma

> 

simplify(new,radical,symbolic);

-16*Pi^2*R^3*sigma^2+8*Pi*R^2*sigma

> 

simplify(new,radical);
factor(%);

4*Pi*R^2*sigma-8*Pi^2*R^3*sigma^2-4*Pi*(R^4*sigma^2*(2*Pi*R*sigma-1)^2)^(1/2)

-4*Pi*(2*Pi*R^3*sigma^2-R^2*sigma+(R^4*sigma^2*(2*Pi*R*sigma-1)^2)^(1/2))

> 

radnormal(new);

-4*Pi*(2*Pi*R^3*sigma^2-R^2*sigma+(R^4*sigma^2*(2*Pi*R*sigma-1)^2)^(1/2))

> 

 

Download radsimp_example.mw

colors...

acer 32385
September 02 2020
1 2

You are using Maple 18 (and not Maple 2018 as your Question's header was incorrectly marked by you). That is why Tom's earlier suggestion did not work for you. I have corrected the Question's header.

Please don't submit duplicates queries for this problem. Post your followups here.

(You should also figure out how to search Mapleprimes -- since I suspect that previous students of your course/supervisor have been tasked with very similar examples and asked here before.)

You can change the name of the color palette (assigned at the top of this attachment). There are a few other choices. (Perhaps see here too.)

 

> 

restart;
> 

with(plots):
> 

palette := "Spring"; # "Dalton" "Niagara"
> 

#
> 

# Define the ODE system
> 

#
> 

  odeSys:= { (diff(F(eta), eta, eta, eta))*(1+epsilon-alpha((diff(F(eta), eta, eta))^2))+F(eta)*(diff(F(eta), eta, eta))+S*(diff(F(eta), eta))-(1/2)*S*eta*(diff(F(eta), eta, eta))-(diff(F(eta), eta))^2-M*(diff(F(eta), eta)), (diff(theta(eta), eta, eta))*(1+R)-delta*(diff(F(eta), eta))^2-Pr((3/2)*S*theta(eta)+(1/2)*S*eta*(diff(theta(eta), eta))-2*(diff(F(eta), eta))*theta(eta)+F*(diff(theta(eta), eta)))}:
> 

#
> 

# Define the first set of boundary conditions
> 

#
> 

  bcs1:= { F(0) = 0, (D(F))(0) = 1, (D(F))(inf) = 0, theta(0) = 1, theta(inf) = 0
> 

         }:
> 

 

> 

  RVals:=[0.1, 0.5, 1]:
> 

  for k from 1 by 1 to numelems(RVals) do
> 

      pList:=[ epsilon = 0.18, M = 0.5, S = 1.5, delta = 0.3, Pr = 1.5, alpha = 0.4, R = RVals[k],inf=1]:
> 

      sol1[k]:= dsolve( eval
> 

                        ( `union`( odeSys, bcs1),
> 

                           pList
> 

                        ),
> 

                        numeric
> 

                      );
> 

    od:
  display
  ( [ seq
      ( odeplot
        ( sol1[i],
          [eta, theta(eta)],
          eta=0..2,
          color=cat(palette," ",i)
        ),
        i=1..numelems(RVals)
      )
    ],
    title = typeset( theta(eta), " versus ", eta),
    titlefont = [times, bold, 20]
  );

"Spring"

Warning, right boundary of range exceeds domain of the problem, it has been adjusted

Warning, right boundary of range exceeds domain of the problem, it has been adjusted

Warning, right boundary of range exceeds domain of the problem, it has been adjusted

> 

 
> 

 p1:= display
      ( [ seq
          ( odeplot
            ( sol1[i],
              [eta, theta(eta)],
              eta=0..1,
              color=cat(palette," ",i)
            ),
            i=1..numelems(RVals)
          )
        ],
        title = typeset( theta(eta), " versus ", eta),
        titlefont = [times, bold, 20]
     ):
> 

 p2:= display
      ( [ seq
          ( odeplot
            ( sol1[i],
              [eta, diff(theta(eta), eta)],
              eta=0..1,
              color=cat(palette," ",i)
            ),
            i=1..numelems(RVals)
          )
        ],
        title = typeset( diff(theta(eta),eta), " versus ", eta),
        titlefont = [times, bold, 20],
        linestyle=solid
     ):
> 

display([p1,p2], title=typeset( theta(eta), " and ", diff(theta(eta),eta), " versus ", eta));

> 

 

NULL

Download lines_ac.mw

 

three ways...

acer 32385
September 01 2020
2 2

You can augment the ODE system with a new differential equation involving g(t) such that g(1) equals the integral of Y0 from 0 to 1. For example, diff(g(t),t)=y0(t) and initial condition g(0)=0.

For example,

restart;
ICs := {y0(0)=3}:
ODESys := {diff(y0(t),t)=y0(t)}:
F := dsolve(ODESys union ICs union {diff(g(t),t)=y0(t), g(0)=0},
            {y0(t),g(t)}, numeric, output=listprocedure):

G := eval(g(t),F):

G(1);
            5.15484531216027

Or you can numerically integrate a Y0 procedure using evalf(Int(....)). For example,

restart;
ICs := {y0(0)=3}:
ODESys := {diff(y0(t),t)=y0(t)}:
F := dsolve(ODESys union ICs,
            {y0(t)}, numeric, output=listprocedure):

Y0 := eval(y0(t),F):

evalf(Int(Y0(t),t=0..1));
            5.15484514547602

I believe that using the augmented ODE system is more often somewhat better (numerically speaking) than numerically integrating the result procedure.

You don't have to use the listprocedure option, however (though I prefer it). There are other valid alternatives for picking off y0(t).

You could also do it with just a slight modification of your original code, using operator-form for the integrand and a floating-point end-point for the range of integration. (I consider this way of picking off the Y0 procedure to be fragile and inferior, and I also prefer my earlier syntax and methodology for the numeric integration.) For example,

restart;
ICs := {y0(0)=3}:
ODESys := {diff(y0(t),t)=y0(t)}:
F := dsolve(ODESys union ICs, {y0(t)}, numeric):

Y0 := t -> rhs(op(2, F(t))):

int(Y0,0..1.0);
            5.15484514547602

a few...

acer 32385
September 01 2020
2 0
> 

restart:
> 

EQ := 1/4*n[1]^4 - 1/2*n[1]^2 + 1/4;

(1/4)*n[1]^4-(1/2)*n[1]^2+1/4

> 

convert(EQ,sqrfree);

(1/4)*(n[1]^2-1)^2

> 

with(Student[Precalculus]):
> 

thaw(CompleteSquare(algsubs(n[1]^2=freeze(n[1]^2),EQ)));

(1/4)*(n[1]^2-1)^2

> 

thaw(simplify(algsubs(n[1]^2=freeze(n[1]^2),EQ),size));

(1/4)*(n[1]^2-1)^2

> 

 

Download sqr_ex.mw

Naturally, the robustness of any technique can depend on the examples provided up front. I don't know what other examples you might have.

another way...

acer 32385
September 01 2020
2 7

I am supposing that you really do want the factor of 1/2 distributed throughout the sum in brackets in the result, otherwise you could take combine(subexp) which returns with the 1/2 factored out front.

> 

restart;
> 

subexp := M__a*sin(omega*t + alpha)*I__a*sin(omega*t + phi);

M__a*sin(omega*t+alpha)*I__a*sin(omega*t+phi)

> 

subexp2 := M__a*I__a*(1/2*(cos(alpha - phi)-cos(2*omega*t + alpha + phi)));

M__a*I__a*((1/2)*cos(alpha-phi)-(1/2)*cos(2*omega*t+alpha+phi))

> 

op(0,subexp)(combine~([selectremove(type,subexp,trig)])[]);

M__a*I__a*((1/2)*cos(alpha-phi)-(1/2)*cos(2*omega*t+alpha+phi))

> 

 

Download combine_trig_ex.mw

note. The order of the additive terms in the bracket depend only on the order in which they appear when subexp2 is first formed.

I could have used `*` instead of op(0,subexp) (since I can see that it is so) but I wanted it more generally programmatically valid and less ad hoc. Ie,

`*`(combine~([selectremove(type,subexp,trig)])[]);

Alternatives -- getting more precise -- might include,

subsindets(subexp,`*`,
           u->`*`(combine~([selectremove(type,subexp,trig)],trig)[]));

subsindets(subexp,And(`*`,satisfies(u->nops(select(type,u,trig))>1)),
           u->`*`(combine~([selectremove(type,subexp,trig)],trig)[]));

solutionmodule...

acer 32385
August 31 2020
0 1
> 

restart;
> 

with(Statistics):
> 

X:=Vector([0,0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.012]):
Y:=Vector([1.103,1.057,1.016,0.978,0.94,0.91,0.88,0.85,0.826,0.8,0.778,0.735]):
> 

res:=Fit((-a*t+b)/(-d+e*t+f*t^3),X,Y,t,output=solutionmodule):-Results():
> 

eval("residualsumofsquares",res);

0.8857314486e-5

> 

Fn:=eval("leastsquaresfunction",res);

(HFloat(5.570317278545915e-8)*t+HFloat(9.299492440664645e-9))/(HFloat(5.6227984603429837e-5)*t^3+HFloat(4.191948978543518e-7)*t+HFloat(8.43004582827281e-9))

> 

plots:-display(plot(X,Y,style=point,color=red,symbol=solidcircle),
               plot(Fn,t=min(X)..max(X),color=blue),
               labels=["X","Y"],size=[400,300]);

> 

eval("parametervalues",res);

[a = HFloat(-5.570317278545915e-8), b = HFloat(9.299492440664645e-9), d = HFloat(-8.43004582827281e-9), e = HFloat(4.191948978543518e-7), f = HFloat(5.6227984603429837e-5)]

> 

res;

["residualmeansquare" = 0.1265330641e-5, "residualsumofsquares" = 0.8857314486e-5, "residualstandarddeviation" = 0.1124869166e-2, "degreesoffreedom" = 7, "parametervalues" = [a = -5.570317278545915*10^(-8), b = 9.299492440664645*10^(-9), d = -8.43004582827281*10^(-9), e = 4.191948978543518*10^(-7), f = 0.56227984603429837e-4], "parametervector" = (Vector(5, {(1) = -0.5570317279e-7, (2) = 0.9299492441e-8, (3) = -0.8430045828e-8, (4) = 0.4191948979e-6, (5) = 0.56227984603429837e-4})), "leastsquaresfunction" = (5.570317278545915*10^(-8)*t+9.299492440664645*10^(-9))/(0.56227984603429837e-4*t^3+4.191948978543518*10^(-7)*t+8.43004582827281*10^(-9)), "residuals" = (Vector[row](12, {(1) = 0.1366412595140698e-3, (2) = 0.16823177724600846e-3, (3) = -0.6784881336476811e-3, (4) = -0.9686555350072457e-3, (5) = 0.1830443681129612e-2, (6) = -0.6708520719832523e-3, (7) = -0.7999853175416627e-3, (8) = 0.11659468614327873e-2, (9) = -0.10090592909720586e-2, (10) = 0.472803591991533e-3, (11) = -0.5626311607461743e-3, (12) = 0.23114746722663337e-3}))]

> 

 

Download fit_solmodule.mw

Paragraph Style...

acer 32385
August 30 2020
2 0

There are various "Paragraph Styles" which you can modify.

There is a "Paragraph Style" names "Maple Output", whose alignment property can be changed from "Center" to "Left".

Here's what you do:

 - On the main menubar, Format -> Styles...

 - In the popup, left combo-box, select "P Maple Output", then click the "Modify" button

 - Change the "Alignment" dropmenu from "Center" to "Left"

Then execute and compute in your Document, and output should appear left-aligned.

You can even save an empty Document with this change as a so-called style sheet, and make it your default for New Documents.

See in the Help system, here and here. 

First 111 112 113 114 115 116 117 Last Page 113 of 336