Advanced Search

acer

32954 Reputation

29 Badges

20 years, 150 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

Contact acer

MaplePrimes Activity


These are answers submitted by acer

nn is posint...

acer 32954
September 13 2020
0 1

It seems as if you are trying to find the discrete minimum, for nn a positive integer in 1..5.

If I have interpreted the goal properly, then I suspect that the following might not surprise anyone (where each xx[ii] attaints its lower bound, for each posint nn).

Hence my use of n-n^2 in the first computation below (where even that min and seq call are avoidable if we realize it is decreasing).

Your original formulation is an abuse of the notation and syntax, IMHO. But below is (just) one way to deal with the programming.

> 

restart;
> 

min(seq(n-n^2, n=1..5));

-20

> 

W := proc(nn)
  local ii, rf;
  if not nn::posint then return 'procname'(args); end if;
  rf := add(xx[ii]^2-nn,ii=1..nn);
  return rf;
end proc:
> 

K := proc(nn::posint)
  Optimization:-Minimize(W(nn),seq(xx[ii]=1..10,ii=1..nn));
end proc:
> 

min(seq(K(i)[1], i=1..5));

-20.

> 

K(5);

[-20., [xx[1] = HFloat(1.0), xx[2] = HFloat(1.0), xx[3] = HFloat(1.0), xx[4] = HFloat(1.0), xx[5] = HFloat(1.0)]]

> 

 

Download discr_opt.mw

like this?...

acer 32954
September 13 2020
2 4

You could try something like this,

> 

restart;
> 

replace_all_C:=proc(expr::anything)
  subsindets(expr, ':-suffixed(_C, posint)',
           u->c[parse(convert(u,string)[3..])]);
end proc:
> 

sol:=dsolve(diff(x(t),t$3) = x(t));

x(t) = _C1*exp(t)+_C2*exp(-(1/2)*t)*sin((1/2)*3^(1/2)*t)+_C3*exp(-(1/2)*t)*cos((1/2)*3^(1/2)*t)

> 

tmp:=replace_all_C(sol)

x(t) = c[1]*exp(t)+c[2]*exp(-(1/2)*t)*sin((1/2)*3^(1/2)*t)+c[3]*exp(-(1/2)*t)*cos((1/2)*3^(1/2)*t)

> 

 

Download q_acc.mw

references...

acer 32954
September 13 2020
0 2

There are several interesting ways to plot the logistic map, and it's not completely clear from your wording what kind of plot you are after (though I might guess that you might be looking for a straightforward solution, as a beginner in Maple).

On this forum, see here, here, here, and here.

You can also see the Bifurcation command, in more recent Maple versions.

See also here, for the cobweb plot example near the bottom.

faq...

acer 32954
September 12 2020
0 1

This is a fequently asked question.

f:=x->x^2+x^3:

g:=unapply(diff(f(x),x),x);

   g := x -> 3*x^2+2*x

g:=D(f);
   g := x -> 3*x^2+2*x

fit...

acer 32954
September 11 2020
2 4

Let me know if this is the kind of thing that you wanted to accomplish.

I used the DirectSearch ver.2 package for the nonlinear fitting (since the NonlinearFit command does local optimization and without a tight set of parameter ranges it can converge to a local optimum that is a measurably worse fit). You can install the DirectSearch ver.2 package in the Maple 2020 GUI from the menubar's Maple Cloud login, or obtain it from the Maple Cloud website here, or for much older Maple versions install from the Application Center here.
 

> 

restart;
> 

with(LinearAlgebra):
> 

#data:=ExcelTools:-Import("D:/Maple/donnees_cycles_sin_lejeunes.xlsx"):
data:=ExcelTools:-Import(cat(kernelopts(homedir),"/mapleprimes",
                             "/donnees_cycles_sin_lejeunes.xlsx")):
> 

PK1Expe:=Column(data,1):
> 

fHz:=3: #fréquence
lambdaS:=0.5:
lambdaD:=0.3:
lambdaTrig:= piecewise(t < 1/fHz, 1+lambdaS*t*fHz, t >= 1/fHz, 1+lambdaS+lambdaD*sin((2*Pi*(t-1/fHz)*fHz))):

with(LinearAlgebra):
tensF:=Matrix([[lambdaTrig,0,0],[0,1/sqrt(lambdaTrig),0],[0,0,1/sqrt(lambdaTrig)]]):
tensdF:=diff(tensF,t):
tensFDev:=tensF-(1/3)*Trace(tensF)*Matrix(3,shape=identity):
tensB:=Multiply(tensF,Transpose(tensF)):
tensL:=Multiply(tensdF,MatrixInverse(tensF)):
tensD:=(1/2)*(tensL+Transpose(tensL)):
tensBe:=Matrix([[Be11(t),0,0],[0,1/sqrt(Be11(t)),0],[0,0,1/sqrt(Be11(t))]]):
tensBeDev:=tensBe-(1/3)*Trace(tensBe)*Matrix(3,shape=identity):

ode:=tensdBe[1,1]=-(2/3)*tensBe[1,1]*Trace(tensD)+Multiply(tensL[1,1],tensBe[1,1])+Multiply(tensBe[1,1],Transpose(tensL[1,1]))-(2/eta)*a*Multiply(tensBeDev[1,1],tensBe[1,1]): ivp:=[ode,Be11(0)=1]:
> 

tensdBe:=diff(tensBe,t):
> 

ode_sol:=dsolve(ivp,numeric,parameters=[a,eta]):
> 

Be11Num:=proc(aValue,etaValue,tValue)
   global __old_a, __old_eta;
   local res;
   if not [aValue,etaValue,tValue]::list(numeric) then
      return 'procname'(args);
   end if;
   if __old_a<>A_value and __old_eta<>etaValue then
      (__old_a,__old_eta) := aValue,etaValue;
      ode_sol('parameters'=[a=aValue,eta=etaValue]);
   end if;
   res:=rhs(ode_sol(tValue)[2]);
   #evalf[8](res);
end proc:
> 

Be11Num(a,eta,t); # returns unevaluated, good

Be11Num(a, eta, t)

> 

ode_sol(parameters=[a=1,eta=0.1]);
ode_sol(2.0)[2], Be11Num(1,0.1,2.0); # these should agree

[a = 1., eta = .1]

Be11(t) = HFloat(1.2342315801626595), HFloat(1.2342315801626595)

> 

ode_sol(parameters=[a=1.4,eta=0.3]);
ode_sol(2.0)[2], Be11Num(1.4,0.3,2.0); # these should agree

[a = 1.4, eta = .3]

Be11(t) = HFloat(1.1815627486067706), HFloat(1.1815627486067706)

> 

tensTemp1:=2*C1*tensB+4*C2*(Trace(tensB)-3)*tensB+6*C3*(Trace(tensB)-3)*tensB-2*C4*MatrixInverse(tensB):
tensTemp1Dev:=tensTemp1-(1/3)*Trace(tensTemp1)*Matrix(3,shape=identity):
tensBe:=Matrix([[Be11Num(a,eta,t),0,0],[0,1/sqrt(Be11Num(a,eta,t)),0],[0,0,1/sqrt(Be11Num(a,eta,t))]]):
tensTemp2Dev:=2*a*tensBe-(1/3)*Trace(2*a*tensBe)*~Matrix(3,shape=identity):
sigma:=tensTemp1Dev+tensTemp2Dev-p*~Matrix(3,shape=identity):
p:=sigma[2,2]+p:
sigma:=tensTemp1Dev+tensTemp2Dev-p*~Matrix(3,shape=identity):
PK:=sigma.Transpose(MatrixInverse(tensF)):
PK:=PK[1,1]: #expression  to fit
> 

vectorT:=[seq(i,i=0..13/3,0.01)]:
vectorT:=Vector(vectorT):
> 

DS2:=CodeTools:-Usage(
        DirectSearch:-DataFit(PK,vectorT,PK1Expe,t,
                              [ a=0.1 .. 10.0, eta=0.1 .. 10.2
                                ,C1=0..1, C2=0..10, C3=-6..6, C4=0..1
                              ])
);

memory used=7.16GiB, alloc change=6.00MiB, cpu time=110.76s, real time=105.76s, gc time=10.59s

[HFloat(5.798836840043786e-7), [C1 = HFloat(0.07091614374135646), C2 = HFloat(0.19080181427660822), C3 = HFloat(-0.13066108619768044), C4 = HFloat(2.090135897538898e-5), a = HFloat(0.6139583977913079), eta = HFloat(0.5024357747664052)], 615]

> 

plots:-display(
   plot(<vectorT|PK1Expe>, style=point, color=blue),
   plot(eval(PK,DS2[2]), t=min(vectorT)..max(vectorT),
   size=[600,300]));

> 

 

Download donnees_cycles_sin_DS.mw

edit...

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

acer 32954
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".

> 

restart
> 

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

ro := 10^(-15)*Unit('m')

(1/1000000000000000)*Units:-Unit(m)

> 

h := 0.658e-21*Unit('MeV'*'s')

0.658e-21*Units:-Unit(MeV*s)

> 

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

98.70000000*Units:-Unit(MeV*s)*Units:-Unit(m/s)/(Units:-Unit(m)*c)

> 

combine(pmin, units)

98.70000000*Units:-Unit(MeV)/c

> 

simplify(pmin)

98.70000000*Units:-Unit(MeV)/c

> 

By the way...

convert(Units:-Unit('h_bar'), units, MeV*s)

0.6582119515e-21*Units:-Unit(MeV*s)

> 

By the way...

cc := ScientificConstants:-Constant('c'); ScientificConstants:-GetUnit(cc); ScientificConstants:-GetValue(cc); evalf(cc)

Units:-Unit(m/s)

299792458

299792458.

> 

with(ScientificConstants)
> 

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

149896229000000000000000*Units:-Unit(h_bar)*Units:-Unit(m/s)/(Units:-Unit(m)*c)

> 

simplify(Pmin)

98.66348941*Units:-Unit(MeV)/c

> 

 

Download unitsMeVs.mw

Numeric Formatting of NULL assignment (e...

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

 using ray transfer matrix analysis

 

 #We use Rads and mm as units

 

 

 #Ray transfer matrix for free space

Distance := proc (d) options operator, arrow; Matrix(2, 2, [1, d, 0, 1]) end proc

proc (d) options operator, arrow; Matrix(2, 2, [1, d, 0, 1]) end proc

NULL

Lens := proc (f) options operator, arrow; Matrix(2, 2, [1, 0, -1/f, 1]) end proc

proc (f) options operator, arrow; Matrix(2, 2, [1, 0, -1/f, 1]) end proc

Geometry := 2; if Geometry = 1 then f_obj := 4.5; f_fluo := 125; f_TIE := 45; d_sample := f_obj; d_max := 250; d_cam := f_TIE; d_interlens := d_max-d_cam end if; if Geometry = 2 then f_obj := 4.5; f_fluo := 125; f_TIE := 45; d_sample := 4.6; d_max := 250; d_interlens := d_max-d_cam end if

250-d_cam

TIE := Distance(d_cam).Lens(f_TIE).Distance(d_interlens).Lens(f_obj).Distance(d_sample).Vector(2, [distance, angle])

Vector[column](%id = 18446883956361111542)

NULL``

d_cam_solve := fsolve(coeff(TIE[1], angle, 1), d_cam, d_cam = 20 .. 100)

``

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

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

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

acer 32954
September 07 2020
1 1

Are you trying to do something like this?

> 

restart
> 

QPT := proc (q) local nc, qp, qpp, i, p, pp; nc := LinearAlgebra:-Dimension(q); qp := Vector(nc); qpp := Vector(nc); for i to nc do qp[i] := nprintf("#mrow(mi(%s),mi(%a));", p, q[i]); qpp[i] := nprintf("#mrow(mi(%s),mi(%a));", pp, q[i]) end do; return qp, qpp end proc:
> 

 

> 

q := Vector([alpha, beta])

q := Vector(2, {(1) = alpha, (2) = beta})

> 

QPT(q)

Vector(2, {(1) = p*alpha, (2) = p*beta}), Vector(2, {(1) = pp*alpha, (2) = pp*beta})

> 

``

Download test_cat_ac.mw

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

> 

restart
> 

QPT := proc (q) local nc, qp, qpp, i, p, pp; nc := LinearAlgebra:-Dimension(q); qp := Vector(nc); qpp := Vector(nc); for i to nc do qp[i] := cat(p, q[i]); qpp[i] := cat(pp, q[i]) end do; return qp, qpp end proc:
> 

r:=Vector([`&alpha;`,`&beta;`]);

r := Vector(2, {(1) = alpha, (2) = beta})

> 

foo := QPT(r)

foo := Vector(2, {(1) = `p&alpha;`, (2) = `p&beta;`}), Vector(2, {(1) = `pp&alpha;`, (2) = `pp&beta;`})

> 

foo[1]; lprint(convert(foo[1], list))

Vector(2, {(1) = `p&alpha;`, (2) = `p&beta;`})

[`p&alpha;`, `p&beta;`]

> 

 

Download test_cat_ac2.mw

one way...

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

acer 32954
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(-(Nb*t-n*t__rev)
            *exp((1/2)*(L+sqrt(-4*C*L*R^2+L^2))*t/(L*R*C)
                 -(1/2)*(Nb*t-n*t__rev)^2/(Nb^2*sigma__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...

acer 32954
September 05 2020
0 2

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

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

something...

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

 

First 116 117 118 119 120 121 122 Last Page 118 of 342