acer

13691 Reputation

29 Badges

12 years, 54 days

On google groups. (comp.soft-sys.math.maple and sci.math.symbolic)

On stackoverflow.

On math.stackexchange.com.

MaplePrimes Activity


These are replies submitted by acer

Do you want map(F[n],Q[n]) instead, which in modern Maple could also be done as F[n]~(Q[n]) ?

If that's the case then you'd need to pass Q in as a list of lists, eg, Do( [ x -> x , y -> y^2] , [ [0,1,2,3] , [4,5,6,7] ] ) .

@jnjn0291 In case anyone was wondering whether all the exact eigenvalues are real and non-negative for all x,y real. (This relates to our ignoring the imaginary components in the floating-point approximations.)

restart;

interface(rtablesize=11):

kernelopts(version);

`Maple 18.02, X86 64 LINUX, Oct 20 2014, Build ID 991181`

M := `<|>`(`<,>`(2, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0),
     `<,>`(-1, 4, -1, 0, -1, 0, 0, 0, 0, 0, -exp(-I*x)),
     `<,>`(0, -1, 2, 0, 0, -1, 0, 0, 0, 0, 0),
     `<,>`(-1, 0, 0, 4, -1, 0, -exp(I*y), -1, 0, 0, 0),
     `<,>`(0, -1, 0, -1, 4, -1, 0, 0, -1, 0, 0),
     `<,>`(0, 0, -1, 0, -1, 4, -1, 0, 0, -1, 0),
     `<,>`(0, 0, 0, -exp(-I*y), 0, -1, 2, 0, 0, 0, 0),
     `<,>`(0, 0, 0, -1, 0, 0, 0, 2, -1, 0, 0),
     `<,>`(0, 0, 0, 0, -1, 0, 0, -1, 4, -1, -1),
     `<,>`(0, 0, 0, 0, 0, -1, 0, 0, -1, 2, 0),
     `<,>`(0, -exp(I*x), 0, 0, 0, 0, 0, 0, -1, 0, 2));

M := Matrix(11, 11, {(1, 1) = 2, (1, 2) = -1, (1, 3) = 0, (1, 4) = -1, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (1, 11) = 0, (2, 1) = -1, (2, 2) = 4, (2, 3) = -1, (2, 4) = 0, (2, 5) = -1, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 0, (2, 11) = -exp(I*x), (3, 1) = 0, (3, 2) = -1, (3, 3) = 2, (3, 4) = 0, (3, 5) = 0, (3, 6) = -1, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0, (3, 10) = 0, (3, 11) = 0, (4, 1) = -1, (4, 2) = 0, (4, 3) = 0, (4, 4) = 4, (4, 5) = -1, (4, 6) = 0, (4, 7) = -exp(-I*y), (4, 8) = -1, (4, 9) = 0, (4, 10) = 0, (4, 11) = 0, (5, 1) = 0, (5, 2) = -1, (5, 3) = 0, (5, 4) = -1, (5, 5) = 4, (5, 6) = -1, (5, 7) = 0, (5, 8) = 0, (5, 9) = -1, (5, 10) = 0, (5, 11) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = -1, (6, 4) = 0, (6, 5) = -1, (6, 6) = 4, (6, 7) = -1, (6, 8) = 0, (6, 9) = 0, (6, 10) = -1, (6, 11) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = -exp(I*y), (7, 5) = 0, (7, 6) = -1, (7, 7) = 2, (7, 8) = 0, (7, 9) = 0, (7, 10) = 0, (7, 11) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = -1, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 2, (8, 9) = -1, (8, 10) = 0, (8, 11) = 0, (9, 1) = 0, (9, 2) = 0, (9, 3) = 0, (9, 4) = 0, (9, 5) = -1, (9, 6) = 0, (9, 7) = 0, (9, 8) = -1, (9, 9) = 4, (9, 10) = -1, (9, 11) = -1, (10, 1) = 0, (10, 2) = 0, (10, 3) = 0, (10, 4) = 0, (10, 5) = 0, (10, 6) = -1, (10, 7) = 0, (10, 8) = 0, (10, 9) = -1, (10, 10) = 2, (10, 11) = 0, (11, 1) = 0, (11, 2) = -exp(-I*x), (11, 3) = 0, (11, 4) = 0, (11, 5) = 0, (11, 6) = 0, (11, 7) = 0, (11, 8) = 0, (11, 9) = -1, (11, 10) = 0, (11, 11) = 2})

LinearAlgebra:-Norm(simplify(M - M^%H)) assuming x::real, y::real;

0

ASD:=combine(LinearAlgebra:-IsDefinite(evalc~(M),
                                       'query'='positive_semidefinite'))
     assuming x::real, y::real;

0 <= 2176-64*cos(x-y)-64*cos(x+y)-1024*cos(y)-1024*cos(x) and 0 <= 25552-88*cos(x-y)-88*cos(x+y)-3184*cos(y)-3184*cos(x)

is(ASD) assuming x::real, y::real;

true

 


Download cplx_semidef.mw

 

Just a note: Your second working example, using an anonymous procedure, is parsed OK as 2D Input but not as 1D Maple Notation.

So if you are going to enter it here as unformatted text, then you might want this instead:

s:=0: seq(proc() global s; s:=s+i end proc(), i=1..10);

@mapleatha One of the things I was trying to get at before was that it often makes sense to write code "defensively (in the sense of "defensive driving").

So, for your case, before one tries to access op(n,g) it can be prudent to check that type(g,`+`) returns true.

 

@mapleatha As mentioned, op(n,g) exists for all n from 1 to nops(g).

That's what nops(g) means -- the number of operands of g.

The nth operand exists if n is between 1 and the number of operands.

For a sum of terms, the terms (summands) are the operands.

So yes, a sum of q terms has q operands. So the nth operand exists for all n from 1 to q.

@jnjn0291 

You're most welcome. I suspect that the people who answered enjoyed the programming challenges of your problem.

I noticed that the display of the upper value of one of the eigenvalues is about 1.999999991 and not quite 2, which presents some difficulties in 3-D plotting as a surface so close to the surface for eigenvalues lambda=2. I opted for transparent filled regions, but there are other ways. For example the boundaries could be wireframe instead of plain surface. Or the corners could get vertical darker line segments. If you want it to look different then please explain just how.

And here I add in filled regions for the spectral bands.

eigenspect2.mw

@mapleatha You can also use the Typesetting:-Supress command to veil the bracketed arguments of the function calls.

(This again was Maple 16.02, but it may well work the same in old Maple 13.)

You only need make the calls to the interface and Typesetting commands once (at the start of the worksheet, say), once you've decided on the format you prefer. I have them interleaved with the statements containg diff calls just because I was illustrating various choices of output format.

restart;

diff(y(t),t,t)-4*diff(y(t),t)+5*y(t);

diff(diff(y(t), t), t)-4*(diff(y(t), t))+5*y(t)

interface(typesetting=extended):

Typesetting:-Settings(useprime=true, prime=t):

diff(y(t),t,t)-4*diff(y(t),t)+5*y(t);

diff(diff(y(t), t), t)-4*(diff(y(t), t))+5*y(t)

Typesetting:-Suppress(y(t));

diff(y(t),t,t)-4*diff(y(t),t)+5*y(t);

diff(diff(y(t), t), t)-4*(diff(y(t), t))+5*y(t)

 


Download diff_typesetting_suppress.mw

In Maple 2017, with the new default of extended typesetting level, the commands would be just this:
 

Typesetting:-Settings(typesetprime=true,prime=t):
Typesetting:-Suppress(y(t));

 

diff(y(t),t,t)-4*diff(y(t),t)+5*y(t);

diff(diff(y(t), t), t)-4*(diff(y(t), t))+5*y(t)

 


Download tp_prime_2017.mw

@Kitonum By giving `option remember` to procedure AllRoots there is a speedup by a factor that is close to the number of seperate surfaces generated.

(Using `option remember, system` would allow eventual collection of remembered values, but it might incur gc() overhead between the multiple calls to generate the surfaces in the since plot3d call.)

The end result is just a tad faster than the code in my Answer.

I also added the option `complex` to the call to fsolve, because there is a problematic central region for one of the higher surfaces where the imaginary components come out not quite as zero (and that confuses AllRoots into grabbing from the wrong root). One could also just apply Re, to get the upper surface without a hole (as I did in my code, which I think may be where it runs a tad slower).

lambda_remembered.mw

@John Fredsted 

That is in a similar vein to the code I gave. I realize that John is addressing the followup question about session dependent ordering of the multiplicative terms (so what I write below is not criticism, but notes for Taro the OP).

[deleted: something wrong that I wrote.]

I tried to do it efficiently with just one selectremove, instead of a select and a remove. But that is not very important. I like that you used select directly on the numerator, while my code was busy fiddling with lists of multiplicands.

I note that your code would match a multiplicative term in the numerator such as sin(Omega), while I strived for more relevant targetting of a term which is polynomial in Omega (or mu) with degree at least one. For the current example the only term in the numerator that matches does happen to have that property, so for this particular example this doesn't matter.

I might also note that selecting according to only Omega (or only mu) might not always get to the mark. A change to make one of the factors be (K*Omega+1) instead, say, can illustrate this. That's why I tried to look at all the multiplicands in the denominator (as itself a product of names).

div2.mw

There are lots of ways to make examples which stymie both codes, of course.

I mention this all just for Taro's understanding.

If I knew this at some point in the past then I'd forgotten it.

The plot command does not always apply evalf to the end-points, when not operating in evalhf mode.

Instead, it can try and use the exact end-points for evaluation of the expression/procedure, and only generate floating-point values from such results.

This is why, when Digits is high, the call plot( (1-x)^(10^9), x=0 .. 10^(-15) )  can goes away for a long time, or hit the stack-limit.

So perhaps the simplest successful variation here would be to use greater working precision as well as floats for the end-points, eg,

restart;
Digits:=20:
plot( (1-x)^(10^9), x=0 .. 10.0^(-15) );

That also explains why the first variant I gave in my Answer also works, since its operator turns its argument x into a float before computing the formula.

restart;
plot(x->eval(evalf[20]((1-x)^(10^9))), 0..10^(-15))

While surely it's easier just to use floating-point values in the supplied range, it's also possible to build a coercing operator from a pre-existing expression.

restart;
expr:=(1-x)^(10^9);

f := unapply( expr, [x::coerce(float, (s::realcons)->evalf(s))] ):

Digits:=20:
plot( f, 0 .. 10^(-15) );

I find it interesting that the expression b is smaller than the expression a according to some metrics, but not perhaps by that used by simplify(...,size).

restart;

a := rho^(epsilon-1)*a__01^(-k)*(N__E2/Omega+N__E1*mu)*(K+1)/(mu):
b:=rho^(epsilon-1)*a__01^(-k)*(N__E2/(Omega*mu)+N__E1)*(K+1):

length(a), `simplify/size/size`(a), MmaTranslator:-Mma:-LeafCount(a);
                            99, 60, 26

length(b), `simplify/size/size`(b), MmaTranslator:-Mma:-LeafCount(b);
                            90, 63, 24

c:=simplify(a,size):
length(c), `simplify/size/size`(c), MmaTranslator:-Mma:-LeafCount(c);
                           103, 55, 26

lprint(c);
  rho^(epsilon-1)*a__01^(-k)*(N__E1*Omega*mu+N__E2)*(K+1)/Omega/mu

I don't think that length would be a good metric, but I've wondered about how useful it might be to have simplify(...,size) optionally use LeafCount.

@John Fredsted Your suggestion to use applyop is of course simple and effective, and I apologize if I sounded at all officious. And since we only have the single example to handle then obtaining the target seems just fine.

I was mostly just wondering how "it" might be approached more generally, but in truth all I really figured out was that the single given example/target pair doesn't really specify thoroughly what "the" more general goal might be.

note: the boilerplate conditional checks in my code are somewhat mirrored, in a direct coding using say applyop or op, by the implicit assumption that the original expression is a product of terms. Again, this aspect may not matter at all to the OP with just this single example. 

Another command that you could experiment with is collect. Reading its help page is worthwhile.

@dgh So it's a Bad Idea to have both c and c[1],etc, in the same expression.

I change unindexed c to cc at the start. I inserted a multiplication symbol in the definition of EI.

I also evaluated sqrt(RRn/RRd) at the results from solving for c[0],..,c[3].

Check that it's what you intended.

Try not suppressing all output with full colons, as it's hidden clues from you (like the very presence of the indexed lam[i].)
 

restart

w := c[0]+c[1]*x+c[2]*exp(cc*x)+c[3]*exp(-cc*x)+c[4]*x^n:

dw := simplify(diff(w, x))

c[1]+c[2]*cc*exp(cc*x)-c[3]*cc*exp(-cc*x)+c[4]*x^(n-1)*n

ddw := simplify(diff(w, x, x))

c[2]*cc^2*exp(cc*x)+c[3]*cc^2*exp(-cc*x)+c[4]*x^(n-2)*n^2-c[4]*x^(n-2)*n

dddw := simplify(diff(w, x, x, x))

c[2]*cc^3*exp(cc*x)-c[3]*cc^3*exp(-cc*x)+c[4]*x^(n-3)*n^3-3*c[4]*x^(n-3)*n^2+2*c[4]*x^(n-3)*n

m := -alpha*x+1:

EI := -x^4*beta[4]-x^3*beta[3]-x^2*beta[2]-x*beta[1]+1:

T := lam^2*(int(m*(R+x), x = x .. L))/L^4+F

lam^2*(-(1/3)*alpha*(L^3-x^3)+(1/2)*(-R*alpha+1)*(L^2-x^2)+R*(L-x))/L^4+F

n := 4

4

RRn := `assuming`([simplify(int(EI*ddw^2, x = 0 .. L)+int(T*dw^2, x = 0 .. L))], [L > 0]):

RRd := `assuming`([simplify(int(m*w^2, x = 0 .. L))], [L > 0]):

NULL

alpha := 0:

Eq1 := subs(x = 0, w = 0):

Eq2 := subs(x = 0, dw = 0):

Eq3 := subs(x = L, ddw = 0):

Eq4 := subs(x = L, dddw = 0):

S := solve({Eq1, Eq2, Eq3, Eq4}, {c[0], c[1], c[2], c[3]}):

cc := `assuming`([simplify(subs(x = 0, sqrt(T/EI)))], [lam > 0])

(1/2)*lam*2^(1/2)

new := eval(sqrt(RRn/RRd), S):

RQ := evalf(subs(lam = 1, new))

3.700371506

NULL


Download Uniform_Beam_1.mw

 

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