acer

32385 Reputation

29 Badges

19 years, 339 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

You can use the right-click Unit Formatting facility, or the Units:-UseUnit command.

In the actual Maple GUI the units appear in upright roman without the double-bracing. They appear that way when the worksheet is inlined here on Mapleprimes because that's the "old" way that they used to appear in Maple.

restart

kernelopts(version)

`Maple 2020.1, X86 64 LINUX, Jun 8 2020, Build ID 1474533`

with(Units:-Simple)

Units:-UseUnit(J/(kg*K))

22.5*Unit('kJ')/((.25*Unit('kg')*100)*Unit('K'))

900.0000000*Units:-Unit(J/(kg*K))

12.5*Unit('kJ')/(.17*Unit('kg')*(113*Unit('degC')))

650.7027590*Units:-Unit(J/(kg*K))

simplify(22.5*Unit('kJ')/((.25*Unit('kg')*100)*Unit('K')))

900.0000000*Units:-Unit(J/(kg*K))

simplify(0.9e-3*Unit('km'^2/('s'^2*'K')))

900.0000*Units:-Unit(J/(kg*K))

restart

with(Units:-Simple)

 

For the next one I used the right-click menu and chose
the item:  "Convert Output Units...", then I chose the
"Enter Unit" item and typed in J/(kg*K) .

(I did not type in words like "joule per kilogram times kelvin".)

 

22.5*Unit('kJ')/((.25*Unit('kg')*100)*Unit('K'))

900.0000000*Units:-Unit(m^2/(s^2*K))

 

Download unitpref.mw

Starting with a RealRange, you can convert to relations (inequalities) and use the solve command. Union and intersection can be expressed pretty well using inert And and Or.

You can use Open ends, or not, and get strict or nonstrict inequailities accordingly.

Of course, there is nothing stopping you from beginning with inequalities. You don't have to begin from RealRange structures.

Depending on how you express the variables parameter to solve, the result can be in RealRange or relation form.

And so on.

restart;

A:=RealRange(Open(-infinity),Open(7));
B:=RealRange(Open(-10),Open(25));
C:=RealRange(Open(15),Open(infinity));

RealRange(-infinity, Open(7))

RealRange(Open(-10), Open(25))

RealRange(Open(15), infinity)

convert(x::A,relation);
lprint(%);

And(-infinity <= x, x < 7)

And(-infinity <= x,x < 7)

solve(convert~(And(x::A,x::B),relation));
lprint(%);

RealRange(Open(-10), Open(7))

RealRange(Open(-10),Open(7))

solve(convert~(And(x::A,x::B),relation),{x});
lprint(%);

{-10 < x, x < 7}

{-10 < x, x < 7}

solve(convert~(And(x::B,x::C),relation));

RealRange(Open(15), Open(25))

#for a in {A,B,C} do
#  for b in {A,B,C} do
#    temp := And(x::a,x::b);
#    print(temp = solve(convert~(temp,relation)));
#    temp := And(x::a,x::b);
#    print(temp = solve(convert~(temp,relation)));
#    print();
#  end do;
#end do;

 

Download solveRealRange.mw

You have used both y(t) and y, and it complains.

Perhaps you intended this:

eq1 := diff(y(t), t) = 1/2*cos(t - 2*y(t));
with(DEtools);
DEplot(eq1, y(t), t = -2 .. 2, y(t) = -2 .. 2);

If you want to substitute for diff(x(t), t, t) you don't have to assign to it.

(I find you intention slightly confusing, since it appears that the substitution is merely the equation from isolating `Eq__&lambda;1`=0 for diff(x(t), t, t), hence the result is zero. If you intended something else then please explain.)

For example,

restart;

`Eq__&lambda;1`:=diff(x(t),t,t)-R*diff(psi(t),t,t)*cos(phi(t))
                 +R*diff(psi(t),t)*diff(phi(t),t)*sin(phi(t))  

diff(diff(x(t), t), t)-R*(diff(diff(psi(t), t), t))*cos(phi(t))+R*(diff(psi(t), t))*(diff(phi(t), t))*sin(phi(t))

eq:=diff(x(t),t,t)=-R*diff(psi(t),t)*diff(phi(t),t)*sin(phi(t))
    +R*diff(psi(t),t,t)*cos(phi(t));

diff(diff(x(t), t), t) = -R*(diff(psi(t), t))*(diff(phi(t), t))*sin(phi(t))+R*(diff(diff(psi(t), t), t))*cos(phi(t))

eval(`Eq__&lambda;1`, eq);

0

isolate(`Eq__&lambda;1`=0, diff(x(t),t,t)); # same as `eq`

 

 

diff(diff(x(t), t), t) = -R*(diff(psi(t), t))*(diff(phi(t), t))*sin(phi(t))+R*(diff(diff(psi(t), t), t))*cos(phi(t))

 

Download subst_example.mw

Outside the loop (right after calling dsolve) you can make this assignment,

   Y:=eval(y(t),lambdaE):

Then, replace all your calls lambdaE(tCompute[i]) by Y(tCompute[i]) in the later code.

The above works because you used the output=listprocedure option when you called the dsolve command.

If you had omitted the output=listprocedure option then you could get avoid the assignment to Y above, and instead replace all your calls lambdaE(tCompute[i]) by eval(y(t),lambdaE(tCompute[i])) in the subsequent code.

I prefer the way with listprocedure, but that's a stylistic coding preference.

[edit] I could add that it can be needlessly inefficient to call lambdaE(tCompute[i])  (or whatever you replace it with) multiple times for the same numeric value of tCompute[i]. Instead, inside your loop, you could simply call it once and assign the result to a temporary variable.
    temp := lambdaE(tCompute[i]);  # or whatever
And then use temp in the multiple locations your formulas.

Here are two ways.

The first has Nu[a] in italics (because I used single left-quotes to make it a name).

The second has Nu[a] in upright roman.

I specified the axes font, only so that the labels renders more clearly on this site. Adjust as you wish.
 

plot(axes=box,gridlines,size=[375,300],
     labelfont=["Tahoma",bold,14],
     labeldirections=[horizontal,vertical],
     labels=[eta,`Nu[a]`(eta)]);

plot(axes=box,gridlines,size=[375,300],
     labelfont=["Tahoma",bold,14],
     labeldirections=[horizontal,vertical],
     labels=[eta,Typesetting:-mrow(Typesetting:-mo("Nu",':-fontstyle'="normal"),
                                   Typesetting:-mfenced(Typesetting:-mi("a"),
                                   ':-open'="&lsqb;",':-close'="&rsqb;"))(eta)]);

 

 

Download indexedname_label.mw

[edit] A couple of simpler ways (square brackets matching the bold, which you are free remove)
   labels=[eta,Typesetting:-mtext("Nu[a]")(eta)]
and,
   labels=[eta,`#mtext("Nu[a]");`(eta)]

What matters more is how you process and display the results in the second worksheet, rather than how they are returned from the worksheet which computes them.

Let's suppose that your first worksheet does,

   return [matrix(G), matrix(A), matrix(DR), figplot]

do that it returns a list of the items.

Let's suppose that your are making a call like,

   RunWorksheet(...);

in the second worksheet. You could try making that instead be, say,

   results := RunWorksheet(...):
   map(print, results):

That will print each of the results separately, so that the plot renders with its proper size and the matrices don't interfere with each others lines.

Notice that the command above are termined with full colons.

note: Why use lowercase matrix (instead of modern Matrix) which has been deprecated for 20 years?

note: Plots render small when you print them as part of an expression sequence.

note: A fancier way to do all this would be to use Embedded Components in the second worksheet, or Tabulate. But I guess you don't want to change it all.

There is no point-probe functionality provided by the GUI for 3D plots.

If you want the results to be subsequently accessible as integers then you could cast the result from Statistics:-Sample to a hardware datatype integer Array. (I used integer[8] below because my Maple is 64bit, and that provides for a wide range.)

[edit] On my machine this is just over 4 times faster than Carl's code to apply the trunc command after converting to a list.

It also keeps the Sample output as an rtable (in case the code otherwise expected that). Adjustments would be easy if you wanted it to preserve the rtable subtype, eg. keep Vector[row] as Vector[row], etc.

restart;

Probs := [0.10, 0.30, 0.20, 0.25, 0.15];
Bins  := <50, 60, 70, 80, 90>:

[.10, .30, .20, .25, .15]

fCarl:= (Bins, P::list(realcons), N::posint)->
    (trunc~@[seq]@Statistics:-Sample)(
        Statistics:-RandomVariable(
            'EmpiricalDistribution'(Bins, 'probabilities'= P)
        ),
        N
    )    
:

fCarl(Bins, Probs, 4);

[80, 90, 60, 90]

fac:=proc(Bins, P::list(realcons), N::posint)
  uses Statistics;
  Array(Sample(
    RandomVariable(':-EmpiricalDistribution'(Bins,':-probabilities'=P)),
    N),'datatype'=':-integer[8]');
end proc:

fac(Bins, Probs, 4);

Vector[row](4, {(1) = 80, (2) = 50, (3) = 60, (4) = 70})

n := 10^6:

S3:= CodeTools:-Usage(fCarl(Bins, Probs, n)):

memory used=83.95MiB, alloc change=18.89MiB, cpu time=499.00ms, real time=415.00ms, gc time=173.85ms

S4 := CodeTools:-Usage(fac(Bins, Probs, n)):

memory used=22.90MiB, alloc change=7.63MiB, cpu time=129.00ms, real time=93.00ms, gc time=73.72ms

((u->map(uu->evalf[3](rhs(uu)/n),u))@Statistics:-Tally)~([S||(3..4)]);

[[0.997e-1, .300, .200, .150, .250], [.100, .300, .200, .151, .250]]

 

Download int_sample.mw

(It might be nice if the in-place calling sequence for Sample allowed for such an Array as its second argument. That might save on the extra rtable construction/collection.)

It is not very helpful to show us a mere picture of your code instead of uploading and attaching your actual worksheet to your Question.

Yes, the old (and now deprecated) table-based array and matrix are somewhat more awkward to print in your scenario.

Note that assigning to the names u1, etc, only after creation of the structures that reference them doesn't alter the structures themselves. But you can also access entries, or compute with them, and get the evaluation (resolution).

If you want to get a new structure that both contains as well as prints with the entries fully evaluated then you can map the eval command over the earlier made structure.

Here are some examples to consider, using both matrix and Matrix.

restart;

with(linalg):

ROD:=(cos(phi))*matrix([[1,0,0],[0,1,0],[0,0,1]])
      +((1-cos(phi))*matrix(3,1,[u1,u2,u3])
        &* transpose(matrix(3,1,[u1,u2,u3]))
      +sin(phi)*matrix([[0,-u3,u2],[u3,0,-u1],[-u2,u1,0]]));

ROD := Typesetting[delayDotProduct](cos(phi), Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1}), true)+`&*`(Typesetting[delayDotProduct](1-cos(phi), Matrix(3, 1, {(1, 1) = u1, (2, 1) = u2, (3, 1) = u3}), true), Matrix(1, 3, {(1, 1) = u1, (1, 2) = u2, (1, 3) = u3}))+Typesetting[delayDotProduct](sin(phi), Matrix(3, 3, {(1, 1) = 0, (1, 2) = -u3, (1, 3) = u2, (2, 1) = u3, (2, 2) = 0, (2, 3) = -u1, (3, 1) = -u2, (3, 2) = u1, (3, 3) = 0}), true)

ROD:=evalm(ROD);

ROD := Matrix(3, 3, {(1, 1) = cos(phi)+(1-cos(phi))*u1^2, (1, 2) = (1-cos(phi))*u1*u2-sin(phi)*u3, (1, 3) = (1-cos(phi))*u1*u3+sin(phi)*u2, (2, 1) = (1-cos(phi))*u1*u2+sin(phi)*u3, (2, 2) = cos(phi)+(1-cos(phi))*u2^2, (2, 3) = (1-cos(phi))*u2*u3-sin(phi)*u1, (3, 1) = (1-cos(phi))*u1*u3-sin(phi)*u2, (3, 2) = (1-cos(phi))*u2*u3+sin(phi)*u1, (3, 3) = cos(phi)+(1-cos(phi))*u3^2})

phi:=Pi/4; u1:=1/sqrt(3); u2:=1/sqrt(3); u3:=1/sqrt(3);

(1/4)*Pi

(1/3)*3^(1/2)

(1/3)*3^(1/2)

(1/3)*3^(1/2)

evalm(ROD);

Matrix(3, 3, {(1, 1) = cos(phi)+(1-cos(phi))*u1^2, (1, 2) = (1-cos(phi))*u1*u2-sin(phi)*u3, (1, 3) = (1-cos(phi))*u1*u3+sin(phi)*u2, (2, 1) = (1-cos(phi))*u1*u2+sin(phi)*u3, (2, 2) = cos(phi)+(1-cos(phi))*u2^2, (2, 3) = (1-cos(phi))*u2*u3-sin(phi)*u1, (3, 1) = (1-cos(phi))*u1*u3-sin(phi)*u2, (3, 2) = (1-cos(phi))*u2*u3+sin(phi)*u1, (3, 3) = cos(phi)+(1-cos(phi))*u3^2})

eval(ROD[1,1],1);

cos(phi)+(1-cos(phi))*u1^2

ROD[1,1];

(1/3)*2^(1/2)+1/3

eval~(ROD);

Matrix(3, 3, {(1, 1) = (1/3)*sqrt(2)+1/3, (1, 2) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (1, 3) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 1) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 2) = (1/3)*sqrt(2)+1/3, (2, 3) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 1) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 2) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (3, 3) = (1/3)*sqrt(2)+1/3})

map(eval,ROD);

Matrix(3, 3, {(1, 1) = (1/3)*sqrt(2)+1/3, (1, 2) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (1, 3) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 1) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 2) = (1/3)*sqrt(2)+1/3, (2, 3) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 1) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 2) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (3, 3) = (1/3)*sqrt(2)+1/3})

evalm(ROD &* matrix([[1,0,0],[0,1,0],[0,0,1]]));

Matrix(3, 3, {(1, 1) = (1/3)*sqrt(2)+1/3, (1, 2) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (1, 3) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 1) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 2) = (1/3)*sqrt(2)+1/3, (2, 3) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 1) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 2) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (3, 3) = (1/3)*sqrt(2)+1/3})

restart;

A:=array([u]);

A := Matrix(1, 1, {(1, 1) = u})

T:=table([1=u]);

table( [( 1 ) = u ] )

u:=17;;

17

evalm(A);

Matrix(1, 1, {(1, 1) = u})

map(eval,A);

Matrix(1, 1, {(1, 1) = 17})

eval(T);

table( [( 1 ) = u ] )

map(eval,eval(T));

table( [( 1 ) = 17 ] )

restart;

with(LinearAlgebra):

ROD:=(cos(phi))*Matrix([[1,0,0],[0,1,0],[0,0,1]])
      +((1-cos(phi))*Matrix(3,1,[u1,u2,u3])
        . Transpose(Matrix(3,1,[u1,u2,u3]))
      +sin(phi)*Matrix([[0,-u3,u2],[u3,0,-u1],[-u2,u1,0]]));

Matrix(3, 3, {(1, 1) = cos(phi)+(1-cos(phi))*u1^2, (1, 2) = (1-cos(phi))*u1*u2-sin(phi)*u3, (1, 3) = (1-cos(phi))*u1*u3+sin(phi)*u2, (2, 1) = (1-cos(phi))*u1*u2+sin(phi)*u3, (2, 2) = cos(phi)+(1-cos(phi))*u2^2, (2, 3) = (1-cos(phi))*u2*u3-sin(phi)*u1, (3, 1) = (1-cos(phi))*u1*u3-sin(phi)*u2, (3, 2) = (1-cos(phi))*u2*u3+sin(phi)*u1, (3, 3) = cos(phi)+(1-cos(phi))*u3^2})

phi:=Pi/4; u1:=1/sqrt(3); u2:=1/sqrt(3); u3:=1/sqrt(3);

(1/4)*Pi

(1/3)*3^(1/2)

(1/3)*3^(1/2)

(1/3)*3^(1/2)

ROD;

Matrix(3, 3, {(1, 1) = cos(phi)+(1-cos(phi))*u1^2, (1, 2) = (1-cos(phi))*u1*u2-sin(phi)*u3, (1, 3) = (1-cos(phi))*u1*u3+sin(phi)*u2, (2, 1) = (1-cos(phi))*u1*u2+sin(phi)*u3, (2, 2) = cos(phi)+(1-cos(phi))*u2^2, (2, 3) = (1-cos(phi))*u2*u3-sin(phi)*u1, (3, 1) = (1-cos(phi))*u1*u3-sin(phi)*u2, (3, 2) = (1-cos(phi))*u2*u3+sin(phi)*u1, (3, 3) = cos(phi)+(1-cos(phi))*u3^2})

rtable_eval(ROD);

Matrix(3, 3, {(1, 1) = (1/3)*sqrt(2)+1/3, (1, 2) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (1, 3) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 1) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 2) = (1/3)*sqrt(2)+1/3, (2, 3) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 1) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 2) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (3, 3) = (1/3)*sqrt(2)+1/3})

eval(ROD[1,1],1);

cos(phi)+(1-cos(phi))*u1^2

ROD[1,1];

(1/3)*2^(1/2)+1/3

eval~(ROD);

Matrix(3, 3, {(1, 1) = (1/3)*sqrt(2)+1/3, (1, 2) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (1, 3) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 1) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 2) = (1/3)*sqrt(2)+1/3, (2, 3) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 1) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 2) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (3, 3) = (1/3)*sqrt(2)+1/3})

map(eval,ROD);

Matrix(3, 3, {(1, 1) = (1/3)*sqrt(2)+1/3, (1, 2) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (1, 3) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 1) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 2) = (1/3)*sqrt(2)+1/3, (2, 3) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 1) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 2) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (3, 3) = (1/3)*sqrt(2)+1/3})

copy(ROD);

Matrix(3, 3, {(1, 1) = (1/3)*sqrt(2)+1/3, (1, 2) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (1, 3) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 1) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (2, 2) = (1/3)*sqrt(2)+1/3, (2, 3) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 1) = 1/3-(1/6)*sqrt(2)-(1/6)*sqrt(2)*sqrt(3), (3, 2) = 1/3-(1/6)*sqrt(2)+(1/6)*sqrt(2)*sqrt(3), (3, 3) = (1/3)*sqrt(2)+1/3})

 

Download evalm_misc.mw

ps. Don't set Digits to 5. That is a bad idea, and can incur unnecessary roundoff and numeric error. Instead, either keep Digits=10 the default (or raise it) and utilize evalf[5] as a command to print results more tersely.

cos(Pi/4) is a constant.

You need your plotted expression to contain a dependent variable.

But we cannot know what you want... eg.

  cos(x*Pi/4)

  cos(x+Pi/4)

  cos(x/4)

  cos(x)

and so on.

You should try and find a way to test (at least one d value) of the data points for correctness without the integrand being an operator.

Also, I think that a log plot looks more meaningful for this example.

restart;

Omega:=0.01:no:=1:Delta:=0:tau0:=1:c:=0:gamma1:=0:
j:=1:k1:=0:t:=10:xr:=1:Gamma:=0.1:

k:=1:
for i from -20 to 20 do
d[i]:=i:
b:=-gamma1*tau0+I*tau0*Delta-2*a*(k-1)*xr:
a:=1+I*c;
c1:=sqrt(conjugate(a))/tau0:
c2:=0.5*((conjugate(b)/sqrt(conjugate(a)))):
lambda1[i]:=2*(Gamma-I*d[i])*(t-k1)/j+(Gamma-I*d[i])^2:
lambda2[i]:=t*(1-j)*k1/j+1/sqrt(2)*(Gamma-I*d[i]):
lambda3[i]:=c1*(t-k1)/j+c1*(Gamma-I*d[i])+c2:
J1:=sqrt(Pi)/sqrt(2)*(1-erf(lambda2[i])):
J1mod:=(Re(J1))^2+(Im(J1))^2:

g1:=0.5*sqrt(Pi)*tau0*exp(c2^2)
    *exp(-conjugate(a)*((k-1)*xr)^2)/(sqrt(conjugate(a)));

F2[i]:=(-sqrt(2)*Int(unapply(exp(-x^2)*erf(sqrt(2)*c1*x+lambda3[i]),x),
                     -lambda2[i]..100,epsilon=1e-5));

J2:=g1*J1+g1*F2[i];

J2mod:=abs(J2)^2:
f[d[i]]:=J2mod:
end do:

ptsN1 := CodeTools:-Usage(
           [seq([d[i], evalf(f[d[i]])], i = -30 .. 30)]
                         ):

memory used=1.02GiB, alloc change=4.00MiB, cpu time=7.86s, real time=7.87s, gc time=1.21s

plot(ptsN1, color = black, linestyle = solid, thickness = 2,
     color = black, axes = boxed, labelfont = ["HELVETICA", 14],
     axis[2]=[mode=log],
     labels = ["d", "S_M(d)"]);

 

Download MoJalal_ac.mw

The key thing for the speed is to make use operator form for the integrand (but that disables discontinuity checking, so it should be double-checked at at least one value of d).

(If I didn't mess it up, Carl's code could be treated similarly... Mo_Jalal_ac2.mw But check why the plot is not identical.)

You can use PrintToString, concatenate the missing piece, then write the string to a file, then fclose.

For example, with xml being assigned the representation of the XML in Maple (in function-call form),

str:=cat("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n",
         XMLTools:-PrintToString(xml)):
FileTools:-Text:-WriteString("outfile.mw",str):
fclose("outfile.mw");

ps. You can also use fclose after PrintToFile, to handle (only) the closing aspect.

You can choose a file from a popup with Maplets, using the Maplets[Elements][FileDialogue] element.

 

Are you computing as evalf(Eigenvectors(?)) on a Matrix with exact integers. If so then much better would be to apply evalf to the Matrix first, and then call Eigenvectors.

Are you trying to estimate the numeric stability of the eigenvector computation using the ConditionNumber command? That relates to solving linear systems, rather than eigen-computations. Instead, see the EigenConditionNumbers command (and its output=conditionvectors option).

If I recall correctly, you can utilize the output of that EigenConditionNumbers command to estimate the error bounds. (See here and here for background. For "machine epsilon" you could use evalhf(DBL_EPSILON) in the hardware float case of Digits<=15 and 1.0*10^(-Digits) in the software float case.)

For your 18x18 example with random integer entries in -5..5 I am not seeing cases where the eigenvector condition numbers are being computed inaccurately at default Digits=15 (and so it seems to me that those can be used to compute error bounds). If that situation is rare then it's not clear to me what is the nature of your problem.

Why not show us your methodology explicitly, where you are seeing a need to raise working precision high in order to get a good estimate of conditioning? I find your Question to be quite unclear as to what commands you have used.

 

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