9 years, 333 days

These are answers submitted by acer


5 hours ago acer 11639

If T has already been assigned that expression then short and simple is,

V := [2,8,14]:
T := 440*(1-exp(-1/tau)):

evalf(seq(T, tau=V));

              173.1265097, 51.7013629, 30.3323769

And if you want it as a list,

evalf([seq(T, tau=V)]);

             [173.1265097, 51.7013629, 30.3323769]

The stated example problem needs neither eval (or subs) nor unapply.


persistence of PLOT qualities

Yesterday at 8:29 PM acer 11639

There are some rendering qualities of plots that get remembered by the Std GUI with respect to a given Execution Group, and are only cleared when the output is removed. One way to remove output is to comment out (or remove) the input, and re-execute. Or execute something which doesn't ouput a plot. Another way involves using the main menubar's Edit->Remove Output.

You may or may not have noticed such persistence before. Here are a few examples to try, in an in Execution Group in a Worksheet:

1) Execute plot(sin) by default you'll get a plot rendered with axes=normal style. So then right-click on the rendered plot in the output, and in the pop-up context menu choose Axes->Boxed.

Now, without deleting the output, edit the command to be plot(cos), and re-execute. The result will have the axes=boxed style. So the axis style has persisted. But removing the output in the ways I described above will clear this aspect.

2) This next example is not about the relatively new size option for 2D plots. Execute a command such as plot(sin) that renders a 2D plot as the output. Now re-size the rendered output with your mouse pointer, using left-click and drag on the plot's border. Now, without deleting the output, edit the 2D plotting command to instead be plot(cos) and re-execute. The new result will be rendered in the same size as the previously dragged output's border. The plot rendering size has persisted. This works for 3D plots too. Removing the output in the ways I described above will clear this aspect.

Ok, so that is the so-called "plot persistence" thing.

I believe a property something like, say, "render all plots of an Execution Group as thumbnails" has accidentally been implemented with similar persistent behaviour. By this I mean that once the Execution Group gets the quality then it can only be cleared by removing the output in the ways I described above. But otherwise it persists for an Execution Group.

You could experiment, to see whether you agree with this characterization. But that is the behaviour that I consistently see in Maple 2015.1. 

I would say that there are two bugs:

A) The thumbnail persistence aspect. This is dire. There ought to be no such persistence affecting all plots in an Exec. Group. This aspect of rendering any plot as a thumbnail ought to have nothing to do with plot persistence.

B) There is no obvious way to obtain the old elided PLOT(...) for the case of a plot as part of an echoed assinment statement or an expression sequence output. This makes A) more serious.

[update] I believe that it's been shown in another Answer that using the size option on the 2D plot command can result in another kind of persistence overriding that the thumbnail apsect. This is just another similar bug, I'd say. Put all three of these in the same Execution Group. The first time it's executed the 2nd and 3rd plots are rendered differently, with each being rendered correctly. But if re-executed without output removal then the 3rd plot is mis-rendered. And the inappropriate sizing persists until all the plot output is removed from that Exec. Group.

test_plot:=plot(sin(x), x=-Pi..Pi);
plots:-display([test_plot], size=[.6,0.2]);
test_plot; # this should always render in the default square manner!

So now there is another bug to report,

C) The size option of 2D plots induces an analogously incorrect persistence of the sizing of plots rendered in the same Execution Group.



August 27 2015 acer 11639


T__sys := 64.47487*Unit('K'): 


Pure Component Parameters





Aspen Plus | Methods | Parameters | Pure Components


enam := "METHANE":

mwgh := 16.04276*Unit('g'/'mol'): 

DHFORM := -17798.79622*Unit('cal'/'mol'): ``

T__c := 190.564*Unit('K'):   ``







Aspen Plus | Methods | Parameters | Pure Components


CPIGDP__1 := 7.953090666*Unit('cal'/('mol'*'K')): ``

CPIGDP__2 := 19.09166906*Unit('cal'/('mol'*'K')):``

CPIGDP__3 := 2086.9*Unit('K'): ``

CPIGDP__4 := 9.936466991*Unit('cal'/('mol'*'K')): 

CPIGDP__5 := 991.96*Unit('K'): ````

CPIGDP__6 := 50*Unit('K'):  ``

CPIGDP__7 := 1500*Unit('K'): ``NULL



Pure Component Ideal Gas Enthalpy Calculation [THRSWT/7]



Ideal Gas Heat Capacity Equation (DIPPR 107)


Aly and Lee 1981


C__p := proc (T) options operator, arrow; C__1a+C__2a*C__3a^2/(T^2*sinh(C__3a/T)^2)+C__4a*C__5a^2/(T^2*cosh(C__5a/T)^2) end proc:````


Required Pure Component Parameters



C__1a := CPIGDP__1:

C__2a := CPIGDP__2:

C__3a := CPIGDP__3:

C__4a := CPIGDP__4:  

C__5a := CPIGDP__5: 

C__6a := CPIGDP__6:

C__7a := CPIGDP__7:``


C__p(T__sys) = 7.953090666*Units:-Unit(('cal')/(('mol')*('K')))````


Ideal Gas Enthalpy



T__ref := 298.15*Unit('K'):


cand := `assuming`([DHFORM+map(int, C__p(T), T = T__ref .. T__2*Unit(K))], [T__2 > convert(T__ref, unit_free)])



eval(cand, T__2 = convert(T__sys, unit_free))


HIG := proc (T__2) options operator, arrow; DHFORM+map(int, C__p(T), T = T__ref .. T__2) end proc

proc (T__2) options operator, arrow; DHFORM+map(int, C__p(T), T = T__ref .. T__2) end proc





convert(%, units, cal/mol)







August 27 2015 acer 11639

It is not appropriate to use unapply in your definiton of R[m]. As it was your R[m] is an operator that itself just generates returns new operators rather than the expressions you want.

You are missing a multiplication sign between a and y[m-1] inside the definition of R[m].

You mistakenly used z(t) inside the int call in the body of Linv, but it should be just z I suspect.

You have lny[m-1] but you probably(?) mean ln(y[m-1]). However your initial condition of y[0]=0 will not work with ln(y[m-1]), as that will throw an exception.

You cannot get a plot unless variable a has some numeric value.

Here is an edited version which I ran in Maple 13.01. Hopefully I found the syntax and basic programming problems. You could review the math as well.


diff(y(x), x, x)+8*(diff(y(x), x))/x+18*ay+4*ylny = 0, y(0) = 1, (D(y))(0) = 0


n := 3:



proc (z) options operator, arrow; int((int(t^2*z, t = 0 .. s))/s^2, s = 0 .. x), x end proc


for m to n do `&varkappa;`[m] := piecewise(m <= 1, 0, m > 1, 1); R[m] := y[m-1]-1+`&varkappa;`[m]-Linv(18*a*y[m-1]+4*y[m-1]*ln(y[m-1])); y[m] := simplify(y[m-1]*`&varkappa;`[m]-R[m]) end do; U := sum(y[n], i = 0 .. n)









There are two kinds of subscripted name.

Some of the keystroke combinations have changed in recent versions. I'm using Maple 2015 below. If you use an older version then check the Help page for topic worksheet,documenting,2DMathShortcutKeys .

They both look the same when prettyprinted (typeset) as 2D Output. But the GUI marks up only one of them always, (ie. actually displays it as subscripted) as 2D Input.

One kind involves indexed names. Internally it looks like x[j] . If you enter it as 2D Input by typing x[j] then the GUI's editor will not automatically display it as subscripted in the 2D Input. The way to get the GUI to display it as subscripted input is to type the keystroke x and then Ctl-Shift-_ and then j . In other words, Control-Shift-Underscore all together. This is one of the keys to your question, I think.

The other kind of subscripted name is a so-called atomic identifier (name). As 2D Input you can enter the special case of an atomic subscripted name by typing x__j . The GUI's editor will immediately render than as the subscripted name.

The subscript j in the indexed name is mutable, by which I mean that the created name differs as j changes. The subscript j in the atomic name is not mutable, and the created name does not use the value of j. This is the other key to your question.

Because you are creating a procedure whose result depends on the parameters i and j you need to use indexed names for your subscripted names. But it seems that you also want your 2D Input to be rendered with those indexed names as subscripted. So use the keyboard shortcut Ctl-Shift-_ as above.

Here's your worksheet, in which I when back and edited the original assignment to u, and did all subscripts in this way:



There is a sequence of expressions

u := P*x[1]*(x[3]^2/R^3-(1-2*nu)/(R*(R+x[3])))/(4*Pi*mu), P*x[2]*(x[3]^2/R^3-(1-2*nu)/(R*(R+x[3])))/(4*Pi*mu), P*(x[3]^2/R^3-(2*(1-nu))/(R*(R+x[3])))/(4*Pi*mu):



e := proc (i, j) options operator, arrow; (1/2)*(diff(u[i], x[j]))+(1/2)*(diff(u[j], x[i])) end proc;

proc (i, j) options operator, arrow; (1/2)*(diff(u[i], x[j]))+(1/2)*(diff(u[j], x[i])) end proc

e(1, 3)


e(1, 1)







August 27 2015 acer 11639

Since you are interested primarily in the plots then numeric approaches should suffice. That is understood.

I'll do all my computations in 64bit Maple 12.02 for MS-Windows, since that appears to be what you are using.

Let's consider Carl's suggestion to use modify your original to use fsolve. For each numeric instance of x and y this generates a 6th degree univariate polynomial in w, which fsolve can solve to get all 6 roots. As Carl points out picking off just one of these roots to obtain the first plot, and then repeaing all the work to pick off the second roots to obtain the second plot, etc, is inefficient. It does 6 times as much work as necessary.

On my Intel i7 it takes about 105 seconds for computing just the first plot with a grid of 39x39 points, and 850 seconds to make all six plots. Together, they look like this:

Now, you could modify that approach so that all 6 roots from fsolve were utilized, and fsolve were only called once for each plotted x and y grid point. I've done this by populating an Array with all the fsolve results (no duplicated effort), and then using plots:-surfdata on each of the 6 layers of that Array.

On the same machine as above this takes 115 seconds, to compute and produce all six plots (which I happen to display together, but you don't have to). Again I used a 39x39 grid of x-y points. I won't inline the plot here, as it looks the same as above. This approach also shows that the imaginary parts of the eigenvalues are zero, at no extra cost.

It's now time to mention that I disgree with the general approach taken so far, of rootfinding the characteristic polynomial instantiated at particular numeric x-y points. It is slower than what is offered by LinearAlgebra:-Eigenvalues for float Matrices instantiated at those numeric x-y points. I really doubt that it is numerically superior.

On the same machine as above it takes 30 seconds to compute all six plots, with one implementation of this idea. The approach is as before, to populate an Array with roots (six at a time), but now the roots come from LinearAlgebra:-Eigenvalues rather than from fsolve.

I'll leave it to you to figure out whether I got all the orientations and ordering of the data correct. You could probably use fsolve to make a sanity check for a couple or x-y points. (I might have accidentally flipped or transposed data, etc, in that last apporach. You might notice that in my last approach above I applied the sort command to each Vector of results from Eigenvalues. Did I do it all just right?)

This is probably a good time to say that it's not 100% clear that any of the six surfaces represent just one particular eigenvalue. By this I mean we know that you want the first surface, say, to represent the "first" eigenvalue as a function of x and y. And perhaps it does, for this example. But in general what if a pair of surfaces touch, and agree for a while, and then split. Which surface now represents which eigenvalue? For your particular problem this difficulty might not arise. But it can happen in practice, and get tricky. (Note too that we may not be able to factor your sixth degree characteristic polynomial in w and thus may not be able to obtain six explicit formulas.)

Lastly, a bit of advice: if you intend on doing simplification then do not introduce floats until you are forced to. I changed your float constants and coefficients to exact rationals. (Hope I did them right. Please check that.) Using the simplify command on an expression mixed with floats and elementary functions and symbolic names does not always work out well (or even properly).


combine units

August 26 2015 acer 11639

If you do not wish to load the Units:-Standard package (which overloads some arithmetic operators and a few common commands) then you can just use simplify or combine(...,units) to convert volts/amperes to the SI unit ohm.

I mention combine(...,units) since simplify might have other effects on your expression, which you may not want.


expr := 1.0*Unit(volt)/(2000.0*Unit(ampere));


combine(expr, units);









August 23 2015 acer 11639

Does some earlier part of your own full application lay down that .txt file, so that you could change how that happens?

If so, then could you make such strings be base64 encoded, and have such .txt maple code instead be something like,




August 22 2015 acer 11639

ee := C1 + C2*(C3/(T*sinh(C3/T)))^2 + C4*(C5/(T*cosh(C5/T)))^2:

ans:=map(int, ee, T=Tref..Tsys, method=ftocms)
     assuming Tref > 0, Tsys > Tref, C3 > 0, C5 > 0:

op(1,ans) + op(2,ans) + convert(op(3,ans),coth) + convert(op(4,ans),tanh);

                                 /    / C3 \       / C3 \\
      -C1 Tref + C1 Tsys + C3 C2 |coth|----| - coth|----||
                                 \    \Tsys/       \Tref//

                 /    / C5 \       / C5 \\
         - C5 C4 |tanh|----| - tanh|----||
                 \    \Tsys/       \Tref//



August 22 2015 acer 11639

Change those instance of



uses LinearAlgebra;

inside all your proc definitions.


square root of a Matrix

August 21 2015 acer 11639

There are quite a few "lesser" problems with the code (loop summation instead of using the add command, inverting the same Matrix many times, repeated computation inside loops of sub-terms which don't depend on the loop-index, etc).

But as far as can see one of the primary efficiency problems is the cost of computing MatrixPower(V, 1/2) for V a float Matrix. Examination indicated to me that was the primary reason why Digits had to be set so high, and also that it was incurring a large port of the time cost.

It turns out that (at present) MatrixPower uses an inaccurate way to compute the square-root of a float Matrix, since (at present) it is just getting redirected to LinearAlgebra:-MatrixFunction to use a general purpose algorithm.

And the only thing that is being done with Vx the computed Matrix-square-root of V is obtaining and utilizing its eigenvalues and eigenvectors. Now, I don't know what your algorithm is supposed to accomplish overall. But for this particular subtask it seems to me that you could do away entirely with that MatrixPower(V, 1/2) call, and compute the eigenvectors and eigenvalues of the Matrix-square-root of V much more directly.

Here's my modification that makes this change. See the comments where the code changed.

This speeds things up enough that the loop-computations have changed from being relatively incidental to now contributing some of the greatest efficiency costs. So I suspect that the overall computation can be made faster still, but  I don't have time to work on it right now.

If there are still places where Digits needs to be raised, for higher N, then certainly that could be done. But such an increase to the working precision should probably only be done for the the duration of the subcomputations where it is necessary. Set Digits back down again, where it doesn't matter.

It's even possible that the eigen-solving might require a higher value for Digits as N gets large. But even then it'll be faster to utilize the simpler Eigenvectors(V) results than to call MatrixPower.

I'm going to log a bug report, that MatrixPower is not doing the best thing itself in this float Matrix example. Note that even if MatrixPower did the right thing here (accurately) it would still be the wrong thing to use: it'd make Maple compute eigenvectors twice, when once is adequate.



August 21 2015 acer 11639

I tried to speed it up, but a smoothish curve for u=0.1 .. 10 still took about an hour on my i7.

The plots look a lot smoother when viewed in the GUI. (Why is the plot rendering so poor for uploaded workeets on this site?)

I changed the first (of the innermost pair) integral to go from 1 to infinity (instead of 100), as that allowed exact integration after a changed of variables and assumptions that matched the outer integration variables' ranges. Let me know it this is not ok. (The integrand is decaying...)

It was a bit awkward to make some of the integrands be procedures, but I found that I needed to avoid the costly expansion of integrands that evalf/Int attempts. So instead I used the expandoff command to prevent such unwanted expansions. (Without disabling the relevant expansions the computing of the outer integral caused memory exhaustion on my computer.)

I found the final computed plot to be strange over the 0.0 .. 0.1 domain. The are also some wiggles over the 0.1 .. 2.0 domain. The latter I tried to resolve by raising Digits (15, and 30, as attempts), and refining the various epsilon tolerances, and replacement of the NAG methods by _Dexp. But that just made things slower, without the plot being smoother over than subdomain.

Note that repeating the exact integration result, after a full run through the worksheet, needs either the right call to expandon or a restart.



c1:=sqrt((s)/Pi/r^3)/2 * (-sin(r*(u-1)^2/4/(s)+Pi/4)+sin(r*(u+1)^2/4/(s))/sqrt(2)):

c2:=Int(c22, y= 1..infinity):

value(student[changevar](Y=sqrt(y-1),c2,Y)) assuming r>1, r<100, u>1, u<10:
alt_c2:=simplify(combine(simplify(%)),size) assuming r>1, r<100, u>1, u<10;


c3_1:=Int(c33, y= 1..2);

c3_2:=simplify(Int(c33, y= 2..infinity, epsilon=1e-6, method=_d01amc),size) assuming r>0, u>0, u<10;

Int((sin((5/32)*r*y*(u+1)^2/(y-1)+(1/4)*Pi)+sin((5/32)*r*y*(u-1)^2/(y-1)+(1/4)*Pi))*((4/5)*exp(-(5/32)*r*y*(u+1)^2)/(Pi*r^2*y*(y-1)^(1/2))-(8/25)*10^(1/2)*(1/(Pi*(y-1)*y^3*r^5))^(1/2)*exp(u+1+(8/5)/(y*r))*erfc((1/8)*(r*y*(u+1)+16/5)*10^(1/2)/(y*r)^(1/2))), y = 1 .. 2)

Int((sin((5/32)*r*y*(u+1)^2/(y-1)+(1/4)*Pi)+sin((5/32)*r*y*(u-1)^2/(y-1)+(1/4)*Pi))*((4/5)*exp(-(5/32)*r*y*(u+1)^2)/(Pi*r^2*y*(y-1)^(1/2))-(8/25)*10^(1/2)*(1/(Pi*(y-1)*y^3*r^5))^(1/2)*exp((1/5)*(8+5*r*y*(u+1))/(y*r))*erfc((1/8)*(r*y*(u+1)+16/5)*10^(1/2)/(y*r)^(1/2))), y = 2 .. infinity, epsilon = 0.1e-5, method = _d01amc)

eval(c3_2, [u=1.2, r=3.1]):


alt_c3_1 := simplify(combine(student[changevar](Y=1/(y-1), c3_1, Y)),size) assuming r>0, u>0, u<10:
alt_c3_1 := Int(op(1,alt_c3_1), Y=1..100, epsilon=1e-6, method=_d01akc);

Int(-(8/25)*(sin((5/32)*r*(Y+1)*(u+1)^2+(1/4)*Pi)+sin((5/32)*r*(Y+1)*(u-1)^2+(1/4)*Pi))*(Y^2*exp((1/5)*(5*(u+1)*(Y+1)*r+8*Y)/(r*(Y+1)))*erfc((1/8)*((16/5+(u+1)*r)*Y+(u+1)*r)*10^(1/2)/(((Y+1)*r)^(1/2)*Y^(1/2)))*Pi*10^(1/2)-(5/2)*exp(-(5/32)*r*(Y+1)*(u+1)^2/Y)*Y^(3/2)*((Y+1)*Pi*r)^(1/2))/(r^2*(Pi*Y*r+Pi*r)^(1/2)*Pi*Y^2*(Y+1)), Y = 1 .. 100, epsilon = 0.1e-5, method = _d01akc)

evalf(eval(alt_c3_1, [u=1.2, r=3.1]));


plot(eval(op(1,alt_c3_1), [u=1.2, r=3.1]), Y=1..100, size=[500,200]);


evalf(eval(g1, [u=1.2, r=3.1]));


evalf(eval(g1, [u=1.2, r=10.6]));


CodeTools:-Usage( plot( eval(g1, [u=1.2]), r=1..100, adaptive=false, numpoints=300, size=[500,200] ) );

memory used=52.17MiB, alloc change=8.00MiB, cpu time=2.54s, real time=2.71s, gc time=0ns

expandoff(sin, cos, exp, erf, erfc);


#g2_pre := Int(unapply(eval(g1, [u=1.2]), r), 1..100, epsilon=1e-5, method=_Dexp):

g2_pre := Int(g1, r=1..100, epsilon=1e-3, method=_Dexp):
CodeTools:-Usage( evalf(eval(g2_pre, [u=4.1])) );

memory used=189.97MiB, alloc change=32.00MiB, cpu time=11.44s, real time=11.19s, gc time=624.00ms


g2:= unapply(g2_pre,u):



CodeTools:-Usage( evalf(g2(4.0)) );

memory used=190.19MiB, alloc change=0 bytes, cpu time=12.56s, real time=12.22s, gc time=733.20ms


CodeTools:-Usage( plot(g2, 0.1..10, adaptive=false, numpoints=100) );

memory used=43.22GiB, alloc change=384.00MiB, cpu time=60.07m, real time=58.65m, gc time=3.32m

CodeTools:-Usage( plot(g2, 0.1..2.0, adaptive=false, numpoints=50) );

memory used=2.44GiB, alloc change=0 bytes, cpu time=2.49m, real time=2.41m, gc time=9.13s






August 20 2015 acer 11639

Previous answers are fine, though this particular plotting example can also be handled by rescaling the expression. And below is another way to avoid "hardware floats".




# Use software floats by setting Digits>evalhf(Digits)=15 ,
# since UseHardwareFloats=deduced by default.

P1 := plot(ee, t=2000..2001, numpoints=30,
           color=blue, style=point, symbol=cross, symbolsize=20):
Digits := 10:

# Use software floats by setting the value of UseHardwareFloats.

P2 := plot(ee, t=2000..2001, numpoints=20, color=green, style=point):

# Rescale the expression.

ff := expand(eval(ee, t=2000+tbar));
P3 := plottools:-transform((x,y)->[x+2000,y])(plot(ff, tbar=0..1, color=red)):


plots:-display( P1, P2, P3, gridlines=false );






August 19 2015 acer 11639

This is what I have so far...

It produces a much simpler form of the result than what you previously obtained from int, and it demonstrates that this is equal to your expected form (involving coth).

It also shows how your short form involving coth can be turned into this result (much simpler than previously obtained) programatically.

But so far I haven't found a way to turn this new (simpler) form of the result into your equivalent form involving coth.



ee := C1 + C2*(C3/(T*sinh(C3/T)))^2 + C4*(C5/(T*cosh(C5/T)))^2;


desired := C1*(Tsys-Tref) + C2*C3*(coth(C3/Tsys)
           - coth(C3/Tref)) - C4*C5*(tanh(C5/Tsys)-tanh(C5/Tref));


# By mapping int over ee (which is a sum of terms) we can obtain a terser answer.
# That's nice, but ideally we'd like to obtain `desired` above (programmatically).

raw := map(int,ee, T=Tref..Tsys)
       assuming Tref > 0, Tsys > Tref, C3 > 0, C5 > 0;


# The difference produces 0, after conversion to exp form.
# It's awkward to need simplify(..,size) followed by plain simplify. But, OK.

simplify(simplify(convert(raw - desired, exp),size));


# Let's consider the 2nd term of `ee`, and see whether we can convert
# that into the 2nd term of `desired`. If we can it might give a hint
# as to how to go the other way...

raw1 := int(C2*(C3/(T*sinh(C3/T)))^2, T=Tref..Tsys)
        assuming Tref > 0, Tsys > Tref, C3 > 0, C5 > 0;


des1 := convert(C2*C3*(coth(C3/Tsys) - coth(C3/Tref)),exp);


# The corresponding pieces simplify to zero, after conversion to exp form. Good.

simplify(raw1 - des1);


# Note that we can successfully convert `des1` back to the
# middle term of `desired`. Good.


# This converts `des1` to `raw1`. OK.
# Hopefully this will inspire me for a way to go the other way...


normal(% - raw1);








August 18 2015 acer 11639



F:=proc(x, {degree::posint:=3})
    if not type(x,numeric) then return 'procname'(args); end if;
    CurveFitting:-ArrayInterpolation(X, Y, x,
                     ':-method'=':-spline', ':-degree'=degree):
end proc:

evalf(Int(F, 0..Pi));

evalf(Int(x->F(x,degree=3), 0..Pi));

evalf(Int(x->F(x,degree=5), 0..Pi));


1 2 3 4 5 6 7 Last Page 1 of 115