acer

32385 Reputation

29 Badges

19 years, 334 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

@Lucas2310 

I used simplex[maximize] instead of Optimization:-Maximize, to get an exact rather than a floating-point result. (I don't understand why you post all your inequaties in calls to PolyhedralSet.)

kernelopts(version);

          Maple 16.02, X86 64 WINDOWS, Nov 18 2012, Build ID 788210

S:=[0 <= (1/25)*b11, 0 <= (1/25)*b22, 0 <= (1/25)*b33, 0 <= (1/25)*b44,
    0 <= -(12/25)*b44+(1/25025)*b99+(4/275)*b2020-(306/5)*b11-(544/35)*b22-(16/5)
    *b33+17050058028932207659/392694612710400000, 0 <= (63/25)*b44-(2/3575)*b99
    -(9/55)*b2020+(12852/25)*b11+(612/5)*b22+(112/5)*b33
    -1859886463352366249/5342783846400000, 0 <= -(182/25)*b44+(1/275)*b99
    +(42/55)*b2020-(9282/5)*b11-(10608/25)*b22-(364/5)
    *b33+10569449395416434161/8630650828800000, 0 <= (273/25)*b44-(4/275)*b99
    -(91/55)*b2020+(15912/5)*b11+(3536/5)*b22
    +(2912/25)*b33-8902222402551338533/4315325414400000, 0 <= (1/25)*b99,
    0 <= -(1001/25)*b44-(2/25)*b99+(273/25)*b2020-(68068/5)*b11-(14586/5)*b22
    -(2288/5)*b33+6773030323238725459/784604620800000,
    0 <= (2574/25)*b44+(3/25)*b99-(182/5)*b2020+(918918/25)*b11
    +(38896/5)*b22+(6006/5)*b33-6049870521396536773/261534873600000,
    0 <= -(3861/25)*b44-(24/175)*b99+(351/5)*b2020-(286416/5)*b11-(2100384/175)*b22
    -(9152/5)*b33+5465968075214691641/152562009600000, 0 <= (4004/25)*b44+(3/25)*b99
    -712028368745424949/18681062400000-(468/5)*b2020+(306306/5)*b11+(63648/5)*b22
    +(48048/25)*b33, 0 <= -(3003/25)*b44-(2/25)*b99+91*b2020-47124*b11-9724*b22
    -1456*b33+4579528093219020391/156920924160000, 0 <= (1638/25)*b44+(1/25)*b99
    -(1638/25)*b2020+(131274/5)*b11+(26928/5)*b22+(4004/5)*b33
    -4234170109213463609/261534873600000, 0 <= -(637/25)*b44-(4/275)*b99
    +(1911/55)*b2020-(259896/25)*b11-(10608/5)*b22-(1568/5)*b33
    +3935028629902940653/616475059200000, 0 <= (168/25)*b44+(1/275)*b99
    -(728/55)*b2020+(13923/5)*b11+(14144/25)*b22+(416/5)*b33
    -3671189442456764071/2157662707200000, 0 <= -(27/25)*b44-(2/3575)*b99
    +(189/55)*b2020-(2268/5)*b11-(459/5)*b22-(336/25)*b33
    +3432194498564787803/12466495641600000, 0 <= (2/25)*b44+(1/25025)*b99
    -(6/11)*b2020+34*b11+(48/7)*b22+b33-3183161946172621261/157077845084160000,
     0 <= (1/25)*b2020]:

sol:=simplex[maximize](-x, S); # if not, then try simplex[maximize](x, S);

           /      3699816803902151           93598739500223  
    sol := { b11 = -----------------, b2020 = --------------, 
            \      12544650979660800          92417930035200  

            74493639125039        53122892719699        2058358715156147           
      b22 = --------------, b33 = --------------, b44 = ----------------, b99 = 0, 
            50235360215040        87574556160000        3314169918259200           

           \ 
      x = 0 }
           / 

map(evalb, eval(S, sol));

    [true, true, true, true, true, true, true, true, true, true, true, true, true, 

      true, true, true, true, true, true, true]

Use Matrix with a capital letter M.

A (lowercase) matrix is of type last_name_eval , which is why you are seeing just that name RES.

Note also that the matrix command is deprecated.

acer

f:=[a,b,c,a,a];

                        [a, b, c, a, a]

[ListTools:-SearchAll(a, f)];

                           [1, 4, 5]

acer

If you have only Maple 16, then could you replace the strict inequalities with nonstrict inequalities to obtain a system for which the simplex[feasible] command could provide a preliminary test?

By "preliminary test" I mean that a false result using the nonstrict constraints would imply an empty feasible set for the tighter, strict constraints. Of course, a true result using the looser nonstrict constraints would not mean that the stricter constraints had any feasible point. (I only mention this weak kind of test because you have Maple 16.)

Let S be the set of (possibly mixed strict and nonstrict) inequalities,

   simplex[feasible](subsindets(S,`<`,u->op(1,u)<=op(2,u)));

The above command quickly returned false in my Maple 16.02 for both your original (smaller) and your first followup (larger) sequence of strict inequalities. 

acer

I'm not exactly sure what output you want. Is it something like this? (If so, then I feel there ought be a more efficient way of generating the frames.)

restart

with(plots):

farey := proc (N::posint) local a, b; sort([op({seq(seq(a/b, a = 0 .. b), b = 1 .. N)})]) end proc:

L := farey(50):

F := []:

for k to nops(L) do F := [op(F), [[L[k], 0], [L[k], 1/denom(L[k])]]] end do:

Ox := 1/2;

1/2

0

M := 1;

1

M1 := M;

1

M2 := 10;

10

N := 10;

10

dM := abs(M2-M1+1)/N;

1

P := [];

[]

N-1

9

for k from 0 to N-1 do print("Frame:", k); M := dM*k+M1; w := 1/(2*M); trule := plottools:-transform(proc (a, b) options operator, arrow; [(1/2)*(a+w-Ox)/w, (1/2)*(b-Oy)/w] end proc); P := [op(P), display(trule(plot(F, Ox-w .. Ox+w, Oy .. Oy+2*w, view = [Ox-w .. Ox+w, Oy .. Oy+2*w], color = black), view = 0 .. 1, axes = none, gridlines = false))] end do:

"Frame:", 0

"Frame:", 1

"Frame:", 2

"Frame:", 3

"Frame:", 4

"Frame:", 5

"Frame:", 6

"Frame:", 7

"Frame:", 8

"Frame:", 9

display(P, insequence = true, axes = none, gridlines = false);

 

NULL

 

Download FareySeq_modif.mw

I can imagine another reasonable interpretation of your question: that you want the first frame to consist of the full shape taking up all the plotted area, and then as it "zooms" the view of the shape gets more and more restricted (in step with the ongoing rescaling) so that for every frame the (shrinking) displayed view takes up precisely the whole plotting area. If this is in fact the case then it is a sort of opposite to what I showed above -- though both zoom in, of course. 

But if you do actually want that other version of "zooming in" then it should not be difficult to use the same basic idea: plottools:-transform and modification of the view.

acer

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.

acer

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.

acer


restart

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

NULL

Pure Component Parameters

 

NULL

REVIEW-1

 

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'):   ``

T__sys/T__c

.3383370941

``

 

CPIGDP-1

 

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

NULL

````NULL

Pure Component Ideal Gas Enthalpy Calculation [THRSWT/7]

 

NULL

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:````

NULL

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')))````

``NULL

Ideal Gas Enthalpy

 

NULL

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

NULL

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

-17798.79622*Units:-Unit(('cal')/('mol'))+33.27573135*Units:-Unit(('m')^2*('kg')/(('s')^2*('mol')))*(T__2-298.1500000)-304.1758224*Units:-Unit(('m')^2*('kg')/(('s')^2*('mol')))*(exp(-.6999496898*(-5963.+10.*T__2)/T__2)-1096.081578)/(exp(4173.800000/T__2)-1.)-106.1572019*Units:-Unit(('m')^2*('kg')/(('s')^2*('mol')))*(exp(1983.920000/T__2)-775.9594666)/(exp(1983.920000/T__2)+1.)

``

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

-17798.79622*Units:-Unit(('cal')/('mol'))-7882.145563*Units:-Unit(('m')^2*('kg')/(('s')^2*('mol')))

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

HIG(T__sys)

-17798.79622*Units:-Unit(('cal')/('mol'))-7882.145562*Units:-Unit(('m')^2*('kg')/(('s')^2*('mol')))

simplify(%);

-82352.30894*Units:-Unit(('m')^2*('kg')/(('s')^2*('mol')))

convert(%, units, cal/mol)

-19682.67422*Units:-Unit(('cal')/('mol'))

``

``NULL


Download Int_Units_modif.mw

acer

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

restart

n := 3:

1/10000

 

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

(1)

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)

2*(3*a*(6*a*(9999/10000+(3/10000)*x^2*a-(1/3750)*x^2*ln(2)-(1/3750)*x^2*ln(5))+(4/3)*(9999/10000+(3/10000)*x^2*a-(1/3750)*x^2*ln(2)-(1/3750)*x^2*ln(5))*ln(9999/10000+(3/10000)*x^2*a-(1/3750)*x^2*ln(2)-(1/3750)*x^2*ln(5)))*x^2+(2/3)*(6*a*(9999/10000+(3/10000)*x^2*a-(1/3750)*x^2*ln(2)-(1/3750)*x^2*ln(5))+(4/3)*(9999/10000+(3/10000)*x^2*a-(1/3750)*x^2*ln(2)-(1/3750)*x^2*ln(5))*ln(9999/10000+(3/10000)*x^2*a-(1/3750)*x^2*ln(2)-(1/3750)*x^2*ln(5)))*x^2*ln((1/2)*(6*a*(9999/10000+(3/10000)*x^2*a-(1/3750)*x^2*ln(2)-(1/3750)*x^2*ln(5))+(4/3)*(9999/10000+(3/10000)*x^2*a-(1/3750)*x^2*ln(2)-(1/3750)*x^2*ln(5))*ln(9999/10000+(3/10000)*x^2*a-(1/3750)*x^2*ln(2)-(1/3750)*x^2*ln(5)))*x^2))*x^2

(2)

with(plots):

 

``

 

Download LE_EQ_m.mw

acer

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:

 

restart

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):

u[1]

(1/4)*P*x[1]*(x[3]^2/R^3-(1-2*nu)/(R*(R+x[3])))/(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)

(1/8)*P*x[1]*(2*x[3]/R^3+(1-2*nu)/(R*(R+x[3])^2))/(Pi*mu)

e(1, 1)

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

 

 

Download Quemodif.mw

acer

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: example_CL.mw



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. example_modif1.mw 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. ex_modif2.mw

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

acer

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.

restart;

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

0.5000000000e-3*Units:-Unit('V')/Units:-Unit('A')

combine(expr, units);

0.5000000000e-3*Units:-Unit('`&Omega;`')

simplify(expr);

0.5000000000e-3*Units:-Unit('`&Omega;`')

 

 

Download ohm.mw

acer

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,

a:="5pa556iL5byP":
a:=StringTools:-Decode(a,':-base64');

acer

restart:

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

acer

Change those instance of

with(LinearAlgebra);

to

uses LinearAlgebra;

inside all your proc definitions.

acer

First 220 221 222 223 224 225 226 Last Page 222 of 336