acer

32400 Reputation

29 Badges

19 years, 345 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Your expression is a function call, rather than a multiplication (product) of the two bracketed terms.

Put either a space or an explicit multiplication symbol between the two bracketed terms. 

This is a common mistake (especially with users working in 2D Input mode since implicit multiplication fosters omission of the space).

n := 4:
L := [$1..n^2]:

A := Matrix(n,n,L);

rtable_redim(A,0..n-1,0..n-1);

Note that this changes A, in-place. You could also supply copy(A) as the first argument to rtable_redim. (I have sometimes relied on that in-place aspect, for efficiency.)

Matrix(1,4,(i,j)->j*x);

                       [x    2 x    3 x    4 x]

Matrix(1,4,(i,j)->unapply(j*x,x));

          [x -> x    x -> 2 x    x -> 3 x    x -> 4 x]

Matrix(1,4,(i,j)->subs(_J=j,x->_J*x));

          [x -> x    x -> 2 x    x -> 3 x    x -> 4 x]

I believe that the problem is due to over-agressive resampling of the data at the GUI rendering stage. That is to say, resampling when an Array within a CURVES substructure (of a PLOT structure) has many entries.

In my 64bit Maple 2015.2 for Linux there appears to be a threshold of 9999 for the rows of such CURVES Arrays, above which I see the more severe, reported clipping. By splitting such very tall m-by-2 Arrays into several smaller Arrays (all below the purported threshold) then the rendering behaviour improves considerably.

The following attachment was executed in that Maple 2015.2.

example_clipping_1.mw

There is still some resampling that occurs at the rendering stage, for more modest dimensions of CURVES Array. The above attachment also splits into yet a smaller size, exhibiting less beating patterning.

I expect that at least some level of resampling is needed at rendering, not only because the width in pixels of the target area may be smaller than the number of data points, but also for the goal of better GUI performance.

 

Is this what you want?

with(Statistics):

X:=Distribution(Normal(0,1)):

CDF(X, 0.7, numeric);

                       0.758036347776927

Quantile(X, %, numeric);

                       0.700000000000293

FWIW, using Maple 2016,

is( r^2 - r*s + s^2, positive ) assuming r>0, s>0, r<=s;

                            true

is( r^2 - r*s + s^2, positive ) assuming r>0, s>0, r>s;

                            true

So then I tried using an additional conjunction. Here's it failing even when I use a simple (new) name `T` with multiple properties inside an `Or`.

new := subs(s=T+r, r^2 - r*s + s^2):

is( new, positive ) assuming r>0, T+r>0, T<=0;

                                  true

is( new, positive ) assuming r>0, T+r>0, T>0;
                                  true

is( new, positive ) assuming r>0, T+r>0, Or(T<0, T>0, T=0);

                                  FAIL

is( new, positive ) assuming r>0, T+r>0, Or(T<=0, T>0);

                                  FAIL

If you are working in the Standard Java GUI, then you could try the following (I used Maple 2016.0).
 

restart;

twostep:=proc(expr::uneval)
  uses Typesetting;
  (Typeset(expr) =
   subsindets(expr,
     And(name,
         satisfies(t->type(eval(t),
           {constant,
            '`*`'({constant,
                   'specfunc'(anything,
                              Units:-Unit)})}))),
     z->Typeset(EV(z))))
   = combine(eval(expr),':-units');
end proc:

a := 1:
b := 3:
c := 5:

twostep(a + b*c);

(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi("b"), Typesetting:-mo("&InvisibleTimes;"), Typesetting:-mi("c")), Typesetting:-mo("&plus;"), Typesetting:-mi("a")) = Typesetting:-mn("3")*Typesetting:-mn("5")+Typesetting:-mn("1")) = 16

A := 70 * Unit(V):

B := 40 * Unit(ohm):

twostep(A/B);

(Typesetting:-mfrac(Typesetting:-mi("A"), Typesetting:-mi("B")) = Typesetting:-mrow(Typesetting:-mn("70"), Typesetting:-mo("&InvisibleTimes;"), Typesetting:-mfenced(Typesetting:-mi("V"), open = "&lobrk;", close = "&robrk;", Typesetting:-msemantics = "Unit"))/Typesetting:-mrow(Typesetting:-mn("40"), Typesetting:-mo("&InvisibleTimes;"), Typesetting:-mfenced(Typesetting:-mi("&Omega;"), open = "&lobrk;", close = "&robrk;", Typesetting:-msemantics = "Unit"))) = (7/4)*Units:-Unit('A')

 


 

Download units2step.mw


 

restart;

#
# I use prettyprint=1 below only because the %n labeling
# makes it easier to visualize the (desired) form of
# the result.
#

interface(prettyprint=1):

q := 20*x^3 + 10*x^2 + 4*x + 1;

Y:=[solve(q,x,explicit)]:

A,B:=selectremove(u->is(u,real),%):

#F:=u->simplify(radnormal(combine(u)),size):
#F:=u->simplify(simplify(u),size):
F:=u->simplify(evalc(u),size):

ans := `*`(map2(`-`,x,A)[])
       * collect(expand(`*`(map2(`-`,x,B)[])),x,F);

interface(prettyprint=3):
evalf(ans);

                                 3       2          
                        q := 20 x  + 10 x  + 4 x + 1

       /    1        7     1\ / 2   /1   1        7  \     1     1        49  
ans := |x + -- %1 - ---- + -| |x  + |- - -- %1 + ----| x + -- + --- %2 + -----
       \    30      6 %1   6/ \     \3   30      6 %1/     15   900      36 %2
                                                                              
       7      1    \                                                          
   + ----- - --- %1|                                                          
     36 %1   180   /                                                          
                                                                              
                                                                              
                         (1/3)                                                
      /            (1/2)\                                                     
%1 := \350 + 105 15     /                                                     
                         (2/3)                                                
      /            (1/2)\                                                     
%2 := \350 + 105 15     /                                                     

(x+.3423840949)*(x^2+.1576159051*x+.1460348210)

interface(prettyprint=1):

## Pick an example.
## I leave aside the special degree-4 case where A
## below is empty (in which case you may which to
## separate the two conjugate pairs in B).

#q := x^4 + x^3 + x - 1;
#q := x^4 + x - 1;
q := x^4 + x^2 + x - 1;

Y:=[solve(q,x,explicit)]:

A,B:=selectremove(u->is(u,real),%):

## You might experiment with various simpifiers
## which get applied to the coefficients of the
## expanded non-linear part, below.

#F:=u->simplify(radnormal(combine(u)),size):
#F:=u->simplify(simplify(u),size):
F:=u->simplify(evalc(u),size):

ans := `*`(map2(`-`,x,A)[])
       * collect(expand(`*`(map2(`-`,x,B)[])),x,F);

interface(prettyprint=3):
evalf(ans);

                                  4    2        
                            q := x  + x  + x - 1

               /    1       10    1\ / 2   /  2   1       10 \     2   1    
ans := (x + 1) |x - - %1 + ---- - -| |x  + |- - + - %1 - ----| x + - + -- %2
               \    6      3 %1   3/ \     \  3   6      3 %1/     3   36   
                                                                            
     100     10    1    \                                                   
   + ---- + ---- - -- %1|                                                   
     9 %2   9 %1   18   /                                                   
                                                                            
                                                                            
                       (1/3)                                                
      /          (1/2)\                                                     
%1 := \44 + 12 69     /                                                     
                       (2/3)                                                
      /          (1/2)\                                                     
%2 := \44 + 12 69     /                                                     

(x+1.)*(x-.5698402912)*(x^2-.4301597088*x+1.754877666)

 


 

Download realfactors.mw


 

restart;

d := s*x*(E*K*q-K*r+K*sigma[1]+r*x)*
     (1-x*(E*K*q-K*r+K*sigma[1]+r*x)
     /(K*sigma[2]*L))/(K*sigma[2])+sigma[1]*x
     -x*(E*K*q-K*r+K*sigma[1]+r*x)/K;

s*x*(E*K*q-K*r+K*sigma[1]+r*x)*(1-x*(E*K*q-K*r+K*sigma[1]+r*x)/(K*sigma[2]*L))/(K*sigma[2])+sigma[1]*x-x*(E*K*q-K*r+K*sigma[1]+r*x)/K

L := ldegree(d,x);

1

x^L * collect(d/x^L,x,u->simplify(u,size));

x*(-s*r^2*x^3/(K^2*sigma[2]^2)-2*s*(E*q-r+sigma[1])*r*x^2/(K*sigma[2]^2)+(-s*(E*q-r+sigma[1])^2*K+r*sigma[2]*(s-sigma[2]))*x/(K*sigma[2]^2)+((-E*q+r)*sigma[2]+s*(E*q-r+sigma[1]))/sigma[2])

 


 

Download poly.mw

In general what you ask is not always possible, since the characteristic polynomial of your Matrix might not have its roots expressible in "closed form". That's a mathematical consquence more than a limitation of any computer algebra system.

Let's suppose, for fun, that your Matrix entries are very sparse in x, and the charcteristic polynomial of your Matrix is a degree 50 polynomial in x, and that by some chance its structure is such that the roots can be expressed exactly, each taking say 200 pages of output. What would you expect to do with that kind of expression? You're almost certainly not going to get any insight by trying to look at it. Subsequent algebraic manipulation of it would likely be expensive and awkward at best.

Now let's suppose instead that you created a procedure which, when supplied a numeric value for x, would return a sorted list or Vector of all the eigenvalues. You could plot them (or their real and imaginary parts), as x varied. You could extract the max or min or them, as x varied. You might even compute lower order rational polynomial approximations of such functions of x (perhaps over a range). Indeed you could do a wide set of interesting computations with such a black-box function of x.

My question to you is: what manner of (qualitative or quantitative) insights from subsequent computations did you hope to attain with a closed form representation of the eigenvalues of the Matrix that you thing you could definitely not attain from some numeric black-box procedure?

(ps. I've noticed that you leave many of the Answers to your Questions on this site no reply. One consequence is that we cannot tell whether you've understood, or not.)

See the attachment:

Pause_examples.mw 

(The three ways shown are using readstat, DEBUG, and Embedded Components.)

That's just numeric roundoff error. Your Matrices have entries on the order of 10^19, and that final Norm value is small roundoff noise.

If you want Aa with shape=symmetric indexing function then you could make it so yourself:

  Aa := Matrix(Aa,shape=symmetric,datatype=float[8]):

It's not a big deal but I notice that Mt is datatype=complex[8]. You might make it float[8] to get rid of the 0.0*I fluff.

Another choice is just to increase the working precision (Digits), but of course there's a slowdown when you move from hardware to software floats.

Another choice is to rescale your data (and saving the scaling coefficient for later adjustment as needed, naturally). Eg,

Kk := Kk/max(map(abs, Kk));

 

For symbolic Matrix inversion Maple will use two kinds of simplifying call.

One is Testzero, which is used to guard against inadvertently selecting a "hidden zero" as a pivot. Choosing as a pivot -- and dividing by -- an expression which could actually simplify to zero would make the answer an evaluation time bomb. By default Testzero uses Normalizer.

The other is Normalizer which is used when reducing a row, in order to keep expression swell down. Linearsolve first checks whether this has been changed from being identical to the `normal` procedure. If so it uses it as is. If not then it tries a little to figure out what might be appropriate (e.g. radnormal, evala, simplify (...,trig), or 'simplify' as larger bore gun).

[edit: ...cont'd] Sometimes LinearSolve will pass off to `solve` or `SolveTools`, but it mostly does similar things with Normalizer.

If you effectively "turn off" Normalizer by setting it to, say x->x then sometimes symbolic/exact LinearSolve can become considerably faster. Potentially dangerously faster as then it may divide by a pivot whose expression would be exact zero if simplified, leading to an invalid result and a division by zero error later.

If you change Normalizer then by default consequence you change Testzero. But in general one still needs zero-testing of candidate pivots, so better to change both separately in the case that Normalizer is being disabled.

Only you will know what you might want for normalizing across the row while pivoting, depending on your example. It can help in final result size, and in practice it can help with speed because adequate zero-testing for pivot selection can slow down if all the candidates have grown enormously complicated.

Another aspect of choice is in how to choose between candidate pivots. It's not possible in general to know which might lead to the simplest result (without trying everything, which is combinatorically impractical). In the exact case you could select according to expression `length` (naive), leaf-count, types of subexpressions, etc.

So you might try experimenting with adjusting Normalizer and Testzero, and/or using similar in your own code.

Why not use the ImportMatrix command?

I copied and rearranged elements of your sample line, to get a two-line file. The following handled both as expected, using Maple 2016.

ImportMatrix("foo.txt", source=csv);

Also, using your call to readline, how about following that with,

[parse(StringTools:-Escape(current_line,regexp))];

Since you attempted it using piecewise, here's an alternate way to handle this example. (It's just for fun, since Carl's suggestion to use plots:-inequal is more straightforward.)
 

restart;

f := x -> -x^2+3*x+2:
l := x -> 4*x:
m := x -> 2*x-4:

P1 := min(f(x), l(x)):
P2 := max(m(x), 0):
plots:-shadebetween(piecewise(P1 > P2, P1, undefined),
                    piecewise(P1 > P2, P2, undefined),
                    x = -1 .. 4, view = -5 .. 10);

Download shadebetween.mw

[edit] For those with older versions without the modern plots:-shadebetween or plots:-inequal functionality, an old technique (shown by Alec Mihailovs and Robert Israel if I recall) can also be used:

restart;
f := x -> -x^2+3*x+2:
l := x -> 4*x:
m := x -> 2*x-4:
P1 := min(f(x), l(x)):
P2 := max(m(x), 0):
fcn1 := piecewise(P1 > P2, P1, undefined):
fcn2 := piecewise(P1 > P2, P2, undefined):
plots:-display(
  plot([f,l,m,0],-1 .. 4),
  plottools:-transform(
    unapply([x,y+fcn2],x,y))(plot(fcn1-fcn2,x=-1 .. 4,
                                  filled=true,color="gray")),
  view=[-1..4,-5..10]);

First 198 199 200 201 202 203 204 Last Page 200 of 336